Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   1:  #region Translated by Jose Antonio De Santiago-Castillo.
   2:   
   3:  //Translated by Jose Antonio De Santiago-Castillo. 
   4:  //E-mail:JAntonioDeSantiago@gmail.com
   5:  //Web: www.DotNumerics.com
   6:  //
   7:  //Fortran to C# Translation.
   8:  //Translated by:
   9:  //F2CSharp Version 0.71 (November 10, 2009)
  10:  //Code Optimizations: None
  11:  //
  12:  #endregion
  13:   
  14:  using System;
  15:  using DotNumerics.FortranLibrary;
  16:   
  17:  namespace DotNumerics.CSLapack
  18:  {
  19:      /// <summary>
  20:      /// -- LAPACK routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLAED8 merges the two sets of eigenvalues together into a single
  27:      /// sorted set.  Then it tries to deflate the size of the problem.
  28:      /// There are two ways in which deflation can occur:  when two or more
  29:      /// eigenvalues are close together or if there is a tiny element in the
  30:      /// Z vector.  For each such occurrence the order of the related secular
  31:      /// equation problem is reduced by one.
  32:      /// 
  33:      ///</summary>
  34:      public class DLAED8
  35:      {
  36:      
  37:   
  38:          #region Dependencies
  39:          
  40:          IDAMAX _idamax; DLAMCH _dlamch; DLAPY2 _dlapy2; DCOPY _dcopy; DLACPY _dlacpy; DLAMRG _dlamrg; DROT _drot; DSCAL _dscal; 
  41:          XERBLA _xerbla;
  42:   
  43:          #endregion
  44:   
  45:   
  46:          #region Fields
  47:          
  48:          const double MONE =  - 1.0E0; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; 
  49:          const double EIGHT = 8.0E0;int I = 0; int IMAX = 0; int J = 0; int JLAM = 0; int JMAX = 0; int JP = 0; int K2 = 0; 
  50:          int N1 = 0;int N1P1 = 0; int N2 = 0; double C = 0; double EPS = 0; double S = 0; double T = 0; double TAU = 0; 
  51:          double TOL = 0;
  52:   
  53:          #endregion
  54:   
  55:          public DLAED8(IDAMAX idamax, DLAMCH dlamch, DLAPY2 dlapy2, DCOPY dcopy, DLACPY dlacpy, DLAMRG dlamrg, DROT drot, DSCAL dscal, XERBLA xerbla)
  56:          {
  57:      
  58:   
  59:              #region Set Dependencies
  60:              
  61:              this._idamax = idamax; this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dcopy = dcopy; this._dlacpy = dlacpy; 
  62:              this._dlamrg = dlamrg;this._drot = drot; this._dscal = dscal; this._xerbla = xerbla; 
  63:   
  64:              #endregion
  65:   
  66:          }
  67:      
  68:          public DLAED8()
  69:          {
  70:      
  71:   
  72:              #region Dependencies (Initialization)
  73:              
  74:              IDAMAX idamax = new IDAMAX();
  75:              LSAME lsame = new LSAME();
  76:              DLAMC3 dlamc3 = new DLAMC3();
  77:              DLAPY2 dlapy2 = new DLAPY2();
  78:              DCOPY dcopy = new DCOPY();
  79:              DLAMRG dlamrg = new DLAMRG();
  80:              DROT drot = new DROT();
  81:              DSCAL dscal = new DSCAL();
  82:              XERBLA xerbla = new XERBLA();
  83:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  84:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  85:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  86:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  87:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  88:              DLACPY dlacpy = new DLACPY(lsame);
  89:   
  90:              #endregion
  91:   
  92:   
  93:              #region Set Dependencies
  94:              
  95:              this._idamax = idamax; this._dlamch = dlamch; this._dlapy2 = dlapy2; this._dcopy = dcopy; this._dlacpy = dlacpy; 
  96:              this._dlamrg = dlamrg;this._drot = drot; this._dscal = dscal; this._xerbla = xerbla; 
  97:   
  98:              #endregion
  99:   
 100:          }
 101:          /// <summary>
 102:          /// Purpose
 103:          /// =======
 104:          /// 
 105:          /// DLAED8 merges the two sets of eigenvalues together into a single
 106:          /// sorted set.  Then it tries to deflate the size of the problem.
 107:          /// There are two ways in which deflation can occur:  when two or more
 108:          /// eigenvalues are close together or if there is a tiny element in the
 109:          /// Z vector.  For each such occurrence the order of the related secular
 110:          /// equation problem is reduced by one.
 111:          /// 
 112:          ///</summary>
 113:          /// <param name="ICOMPQ">
 114:          /// (input) INTEGER
 115:          /// = 0:  Compute eigenvalues only.
 116:          /// = 1:  Compute eigenvectors of original dense symmetric matrix
 117:          /// also.  On entry, Q contains the orthogonal matrix used
 118:          /// to reduce the original matrix to tridiagonal form.
 119:          ///</param>
 120:          /// <param name="K">
 121:          /// (output) INTEGER
 122:          /// The number of non-deflated eigenvalues, and the order of the
 123:          /// related secular equation.
 124:          ///</param>
 125:          /// <param name="N">
 126:          /// (input) INTEGER
 127:          /// The dimension of the symmetric tridiagonal matrix.  N .GE. 0.
 128:          ///</param>
 129:          /// <param name="QSIZ">
 130:          /// (input) INTEGER
 131:          /// The dimension of the orthogonal matrix used to reduce
 132:          /// the full matrix to tridiagonal form.  QSIZ .GE. N if ICOMPQ = 1.
 133:          ///</param>
 134:          /// <param name="D">
 135:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 136:          /// On entry, the eigenvalues of the two submatrices to be
 137:          /// combined.  On exit, the trailing (N-K) updated eigenvalues
 138:          /// (those which were deflated) sorted into increasing order.
 139:          ///</param>
 140:          /// <param name="Q">
 141:          /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
 142:          /// If ICOMPQ = 0, Q is not referenced.  Otherwise,
 143:          /// on entry, Q contains the eigenvectors of the partially solved
 144:          /// system which has been previously updated in matrix
 145:          /// multiplies with other partially solved eigensystems.
 146:          /// On exit, Q contains the trailing (N-K) updated eigenvectors
 147:          /// (those which were deflated) in its last N-K columns.
 148:          ///</param>
 149:          /// <param name="LDQ">
 150:          /// (input) INTEGER
 151:          /// The leading dimension of the array Q.  LDQ .GE. max(1,N).
 152:          ///</param>
 153:          /// <param name="INDXQ">
 154:          /// (input) INTEGER array, dimension (N)
 155:          /// The permutation which separately sorts the two sub-problems
 156:          /// in D into ascending order.  Note that elements in the second
 157:          /// half of this permutation must first have CUTPNT added to
 158:          /// their values in order to be accurate.
 159:          ///</param>
 160:          /// <param name="RHO">
 161:          /// (input/output) DOUBLE PRECISION
 162:          /// On entry, the off-diagonal element associated with the rank-1
 163:          /// cut which originally split the two submatrices which are now
 164:          /// being recombined.
 165:          /// On exit, RHO has been modified to the value required by
 166:          /// DLAED3.
 167:          ///</param>
 168:          /// <param name="CUTPNT">
 169:          /// (input) INTEGER
 170:          /// The location of the last eigenvalue in the leading
 171:          /// sub-matrix.  min(1,N) .LE. CUTPNT .LE. N.
 172:          ///</param>
 173:          /// <param name="Z">
 174:          /// (input) DOUBLE PRECISION array, dimension (N)
 175:          /// On entry, Z contains the updating vector (the last row of
 176:          /// the first sub-eigenvector matrix and the first row of the
 177:          /// second sub-eigenvector matrix).
 178:          /// On exit, the contents of Z are destroyed by the updating
 179:          /// process.
 180:          ///</param>
 181:          /// <param name="DLAMDA">
 182:          /// (output) DOUBLE PRECISION array, dimension (N)
 183:          /// A copy of the first K eigenvalues which will be used by
 184:          /// DLAED3 to form the secular equation.
 185:          ///</param>
 186:          /// <param name="Q2">
 187:          /// (output) DOUBLE PRECISION array, dimension (LDQ2,N)
 188:          /// If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
 189:          /// a copy of the first K eigenvectors which will be used by
 190:          /// DLAED7 in a matrix multiply (DGEMM) to update the new
 191:          /// eigenvectors.
 192:          ///</param>
 193:          /// <param name="LDQ2">
 194:          /// (input) INTEGER
 195:          /// The leading dimension of the array Q2.  LDQ2 .GE. max(1,N).
 196:          ///</param>
 197:          /// <param name="W">
 198:          /// (output) DOUBLE PRECISION array, dimension (N)
 199:          /// The first k values of the final deflation-altered z-vector and
 200:          /// will be passed to DLAED3.
 201:          ///</param>
 202:          /// <param name="PERM">
 203:          /// (output) INTEGER array, dimension (N)
 204:          /// The permutations (from deflation and sorting) to be applied
 205:          /// to each eigenblock.
 206:          ///</param>
 207:          /// <param name="GIVPTR">
 208:          /// (output) INTEGER
 209:          /// The number of Givens rotations which took place in this
 210:          /// subproblem.
 211:          ///</param>
 212:          /// <param name="GIVCOL">
 213:          /// (output) INTEGER array, dimension (2, N)
 214:          /// Each pair of numbers indicates a pair of columns to take place
 215:          /// in a Givens rotation.
 216:          ///</param>
 217:          /// <param name="GIVNUM">
 218:          /// (output) DOUBLE PRECISION array, dimension (2, N)
 219:          /// Each number indicates the S value to be used in the
 220:          /// corresponding Givens rotation.
 221:          ///</param>
 222:          /// <param name="INDXP">
 223:          /// (workspace) INTEGER array, dimension (N)
 224:          /// The permutation used to place deflated values of D at the end
 225:          /// of the array.  INDXP(1:K) points to the nondeflated D-values
 226:          /// and INDXP(K+1:N) points to the deflated eigenvalues.
 227:          ///</param>
 228:          /// <param name="INDX">
 229:          /// (workspace) INTEGER array, dimension (N)
 230:          /// The permutation used to sort the contents of D into ascending
 231:          /// order.
 232:          ///</param>
 233:          /// <param name="INFO">
 234:          /// (output) INTEGER
 235:          /// = 0:  successful exit.
 236:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 237:          ///</param>
 238:          public void Run(int ICOMPQ, ref int K, int N, int QSIZ, ref double[] D, int offset_d, ref double[] Q, int offset_q
 239:                           , int LDQ, ref int[] INDXQ, int offset_indxq, ref double RHO, int CUTPNT, ref double[] Z, int offset_z, ref double[] DLAMDA, int offset_dlamda
 240:                           , ref double[] Q2, int offset_q2, int LDQ2, ref double[] W, int offset_w, ref int[] PERM, int offset_perm, ref int GIVPTR, ref int[] GIVCOL, int offset_givcol
 241:                           , ref double[] GIVNUM, int offset_givnum, ref int[] INDXP, int offset_indxp, ref int[] INDX, int offset_indx, ref int INFO)
 242:          {
 243:   
 244:              #region Array Index Correction
 245:              
 246:               int o_d = -1 + offset_d;  int o_q = -1 - LDQ + offset_q;  int o_indxq = -1 + offset_indxq;  int o_z = -1 + offset_z; 
 247:               int o_dlamda = -1 + offset_dlamda; int o_q2 = -1 - LDQ2 + offset_q2;  int o_w = -1 + offset_w; 
 248:               int o_perm = -1 + offset_perm; int o_givcol = -3 + offset_givcol;  int o_givnum = -3 + offset_givnum; 
 249:               int o_indxp = -1 + offset_indxp; int o_indx = -1 + offset_indx; 
 250:   
 251:              #endregion
 252:   
 253:   
 254:              #region Prolog
 255:              
 256:              // *
 257:              // *  -- LAPACK routine (version 3.1) --
 258:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 259:              // *     November 2006
 260:              // *
 261:              // *     .. Scalar Arguments ..
 262:              // *     ..
 263:              // *     .. Array Arguments ..
 264:              // *     ..
 265:              // *
 266:              // *  Purpose
 267:              // *  =======
 268:              // *
 269:              // *  DLAED8 merges the two sets of eigenvalues together into a single
 270:              // *  sorted set.  Then it tries to deflate the size of the problem.
 271:              // *  There are two ways in which deflation can occur:  when two or more
 272:              // *  eigenvalues are close together or if there is a tiny element in the
 273:              // *  Z vector.  For each such occurrence the order of the related secular
 274:              // *  equation problem is reduced by one.
 275:              // *
 276:              // *  Arguments
 277:              // *  =========
 278:              // *
 279:              // *  ICOMPQ  (input) INTEGER
 280:              // *          = 0:  Compute eigenvalues only.
 281:              // *          = 1:  Compute eigenvectors of original dense symmetric matrix
 282:              // *                also.  On entry, Q contains the orthogonal matrix used
 283:              // *                to reduce the original matrix to tridiagonal form.
 284:              // *
 285:              // *  K      (output) INTEGER
 286:              // *         The number of non-deflated eigenvalues, and the order of the
 287:              // *         related secular equation.
 288:              // *
 289:              // *  N      (input) INTEGER
 290:              // *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
 291:              // *
 292:              // *  QSIZ   (input) INTEGER
 293:              // *         The dimension of the orthogonal matrix used to reduce
 294:              // *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
 295:              // *
 296:              // *  D      (input/output) DOUBLE PRECISION array, dimension (N)
 297:              // *         On entry, the eigenvalues of the two submatrices to be
 298:              // *         combined.  On exit, the trailing (N-K) updated eigenvalues
 299:              // *         (those which were deflated) sorted into increasing order.
 300:              // *
 301:              // *  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
 302:              // *         If ICOMPQ = 0, Q is not referenced.  Otherwise,
 303:              // *         on entry, Q contains the eigenvectors of the partially solved
 304:              // *         system which has been previously updated in matrix
 305:              // *         multiplies with other partially solved eigensystems.
 306:              // *         On exit, Q contains the trailing (N-K) updated eigenvectors
 307:              // *         (those which were deflated) in its last N-K columns.
 308:              // *
 309:              // *  LDQ    (input) INTEGER
 310:              // *         The leading dimension of the array Q.  LDQ >= max(1,N).
 311:              // *
 312:              // *  INDXQ  (input) INTEGER array, dimension (N)
 313:              // *         The permutation which separately sorts the two sub-problems
 314:              // *         in D into ascending order.  Note that elements in the second
 315:              // *         half of this permutation must first have CUTPNT added to
 316:              // *         their values in order to be accurate.
 317:              // *
 318:              // *  RHO    (input/output) DOUBLE PRECISION
 319:              // *         On entry, the off-diagonal element associated with the rank-1
 320:              // *         cut which originally split the two submatrices which are now
 321:              // *         being recombined.
 322:              // *         On exit, RHO has been modified to the value required by
 323:              // *         DLAED3.
 324:              // *
 325:              // *  CUTPNT (input) INTEGER
 326:              // *         The location of the last eigenvalue in the leading
 327:              // *         sub-matrix.  min(1,N) <= CUTPNT <= N.
 328:              // *
 329:              // *  Z      (input) DOUBLE PRECISION array, dimension (N)
 330:              // *         On entry, Z contains the updating vector (the last row of
 331:              // *         the first sub-eigenvector matrix and the first row of the
 332:              // *         second sub-eigenvector matrix).
 333:              // *         On exit, the contents of Z are destroyed by the updating
 334:              // *         process.
 335:              // *
 336:              // *  DLAMDA (output) DOUBLE PRECISION array, dimension (N)
 337:              // *         A copy of the first K eigenvalues which will be used by
 338:              // *         DLAED3 to form the secular equation.
 339:              // *
 340:              // *  Q2     (output) DOUBLE PRECISION array, dimension (LDQ2,N)
 341:              // *         If ICOMPQ = 0, Q2 is not referenced.  Otherwise,
 342:              // *         a copy of the first K eigenvectors which will be used by
 343:              // *         DLAED7 in a matrix multiply (DGEMM) to update the new
 344:              // *         eigenvectors.
 345:              // *
 346:              // *  LDQ2   (input) INTEGER
 347:              // *         The leading dimension of the array Q2.  LDQ2 >= max(1,N).
 348:              // *
 349:              // *  W      (output) DOUBLE PRECISION array, dimension (N)
 350:              // *         The first k values of the final deflation-altered z-vector and
 351:              // *         will be passed to DLAED3.
 352:              // *
 353:              // *  PERM   (output) INTEGER array, dimension (N)
 354:              // *         The permutations (from deflation and sorting) to be applied
 355:              // *         to each eigenblock.
 356:              // *
 357:              // *  GIVPTR (output) INTEGER
 358:              // *         The number of Givens rotations which took place in this
 359:              // *         subproblem.
 360:              // *
 361:              // *  GIVCOL (output) INTEGER array, dimension (2, N)
 362:              // *         Each pair of numbers indicates a pair of columns to take place
 363:              // *         in a Givens rotation.
 364:              // *
 365:              // *  GIVNUM (output) DOUBLE PRECISION array, dimension (2, N)
 366:              // *         Each number indicates the S value to be used in the
 367:              // *         corresponding Givens rotation.
 368:              // *
 369:              // *  INDXP  (workspace) INTEGER array, dimension (N)
 370:              // *         The permutation used to place deflated values of D at the end
 371:              // *         of the array.  INDXP(1:K) points to the nondeflated D-values
 372:              // *         and INDXP(K+1:N) points to the deflated eigenvalues.
 373:              // *
 374:              // *  INDX   (workspace) INTEGER array, dimension (N)
 375:              // *         The permutation used to sort the contents of D into ascending
 376:              // *         order.
 377:              // *
 378:              // *  INFO   (output) INTEGER
 379:              // *          = 0:  successful exit.
 380:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 381:              // *
 382:              // *  Further Details
 383:              // *  ===============
 384:              // *
 385:              // *  Based on contributions by
 386:              // *     Jeff Rutter, Computer Science Division, University of California
 387:              // *     at Berkeley, USA
 388:              // *
 389:              // *  =====================================================================
 390:              // *
 391:              // *     .. Parameters ..
 392:              // *     ..
 393:              // *     .. Local Scalars ..
 394:              // *
 395:              // *     ..
 396:              // *     .. External Functions ..
 397:              // *     ..
 398:              // *     .. External Subroutines ..
 399:              // *     ..
 400:              // *     .. Intrinsic Functions ..
 401:              //      INTRINSIC          ABS, MAX, MIN, SQRT;
 402:              // *     ..
 403:              // *     .. Executable Statements ..
 404:              // *
 405:              // *     Test the input parameters.
 406:              // *
 407:   
 408:              #endregion
 409:   
 410:   
 411:              #region Body
 412:              
 413:              INFO = 0;
 414:              // *
 415:              if (ICOMPQ < 0 || ICOMPQ > 1)
 416:              {
 417:                  INFO =  - 1;
 418:              }
 419:              else
 420:              {
 421:                  if (N < 0)
 422:                  {
 423:                      INFO =  - 3;
 424:                  }
 425:                  else
 426:                  {
 427:                      if (ICOMPQ == 1 && QSIZ < N)
 428:                      {
 429:                          INFO =  - 4;
 430:                      }
 431:                      else
 432:                      {
 433:                          if (LDQ < Math.Max(1, N))
 434:                          {
 435:                              INFO =  - 7;
 436:                          }
 437:                          else
 438:                          {
 439:                              if (CUTPNT < Math.Min(1, N) || CUTPNT > N)
 440:                              {
 441:                                  INFO =  - 10;
 442:                              }
 443:                              else
 444:                              {
 445:                                  if (LDQ2 < Math.Max(1, N))
 446:                                  {
 447:                                      INFO =  - 14;
 448:                                  }
 449:                              }
 450:                          }
 451:                      }
 452:                  }
 453:              }
 454:              if (INFO != 0)
 455:              {
 456:                  this._xerbla.Run("DLAED8",  - INFO);
 457:                  return;
 458:              }
 459:              // *
 460:              // *     Quick return if possible
 461:              // *
 462:              if (N == 0) return;
 463:              // *
 464:              N1 = CUTPNT;
 465:              N2 = N - N1;
 466:              N1P1 = N1 + 1;
 467:              // *
 468:              if (RHO < ZERO)
 469:              {
 470:                  this._dscal.Run(N2, MONE, ref Z, N1P1 + o_z, 1);
 471:              }
 472:              // *
 473:              // *     Normalize z so that norm(z) = 1
 474:              // *
 475:              T = ONE / Math.Sqrt(TWO);
 476:              for (J = 1; J <= N; J++)
 477:              {
 478:                  INDX[J + o_indx] = J;
 479:              }
 480:              this._dscal.Run(N, T, ref Z, offset_z, 1);
 481:              RHO = Math.Abs(TWO * RHO);
 482:              // *
 483:              // *     Sort the eigenvalues into increasing order
 484:              // *
 485:              for (I = CUTPNT + 1; I <= N; I++)
 486:              {
 487:                  INDXQ[I + o_indxq] = INDXQ[I + o_indxq] + CUTPNT;
 488:              }
 489:              for (I = 1; I <= N; I++)
 490:              {
 491:                  DLAMDA[I + o_dlamda] = D[INDXQ[I + o_indxq] + o_d];
 492:                  W[I + o_w] = Z[INDXQ[I + o_indxq] + o_z];
 493:              }
 494:              I = 1;
 495:              J = CUTPNT + 1;
 496:              this._dlamrg.Run(N1, N2, DLAMDA, offset_dlamda, 1, 1, ref INDX, offset_indx);
 497:              for (I = 1; I <= N; I++)
 498:              {
 499:                  D[I + o_d] = DLAMDA[INDX[I + o_indx] + o_dlamda];
 500:                  Z[I + o_z] = W[INDX[I + o_indx] + o_w];
 501:              }
 502:              // *
 503:              // *     Calculate the allowable deflation tolerence
 504:              // *
 505:              IMAX = this._idamax.Run(N, Z, offset_z, 1);
 506:              JMAX = this._idamax.Run(N, D, offset_d, 1);
 507:              EPS = this._dlamch.Run("Epsilon");
 508:              TOL = EIGHT * EPS * Math.Abs(D[JMAX + o_d]);
 509:              // *
 510:              // *     If the rank-1 modifier is small enough, no more needs to be done
 511:              // *     except to reorganize Q so that its columns correspond with the
 512:              // *     elements in D.
 513:              // *
 514:              if (RHO * Math.Abs(Z[IMAX + o_z]) <= TOL)
 515:              {
 516:                  K = 0;
 517:                  if (ICOMPQ == 0)
 518:                  {
 519:                      for (J = 1; J <= N; J++)
 520:                      {
 521:                          PERM[J + o_perm] = INDXQ[INDX[J + o_indx] + o_indxq];
 522:                      }
 523:                  }
 524:                  else
 525:                  {
 526:                      for (J = 1; J <= N; J++)
 527:                      {
 528:                          PERM[J + o_perm] = INDXQ[INDX[J + o_indx] + o_indxq];
 529:                          this._dcopy.Run(QSIZ, Q, 1+(PERM[J + o_perm]) * LDQ + o_q, 1, ref Q2, 1+J * LDQ2 + o_q2, 1);
 530:                      }
 531:                      this._dlacpy.Run("A", QSIZ, N, Q2, 1+1 * LDQ2 + o_q2, LDQ2, ref Q, 1+1 * LDQ + o_q
 532:                                       , LDQ);
 533:                  }
 534:                  return;
 535:              }
 536:              // *
 537:              // *     If there are multiple eigenvalues then the problem deflates.  Here
 538:              // *     the number of equal eigenvalues are found.  As each equal
 539:              // *     eigenvalue is found, an elementary reflector is computed to rotate
 540:              // *     the corresponding eigensubspace so that the corresponding
 541:              // *     components of Z are zero in this new basis.
 542:              // *
 543:              K = 0;
 544:              GIVPTR = 0;
 545:              K2 = N + 1;
 546:              for (J = 1; J <= N; J++)
 547:              {
 548:                  if (RHO * Math.Abs(Z[J + o_z]) <= TOL)
 549:                  {
 550:                      // *
 551:                      // *           Deflate due to small z component.
 552:                      // *
 553:                      K2 = K2 - 1;
 554:                      INDXP[K2 + o_indxp] = J;
 555:                      if (J == N) goto LABEL110;
 556:                  }
 557:                  else
 558:                  {
 559:                      JLAM = J;
 560:                      goto LABEL80;
 561:                  }
 562:              }
 563:          LABEL80:;
 564:              J = J + 1;
 565:              if (J > N) goto LABEL100;
 566:              if (RHO * Math.Abs(Z[J + o_z]) <= TOL)
 567:              {
 568:                  // *
 569:                  // *        Deflate due to small z component.
 570:                  // *
 571:                  K2 = K2 - 1;
 572:                  INDXP[K2 + o_indxp] = J;
 573:              }
 574:              else
 575:              {
 576:                  // *
 577:                  // *        Check if eigenvalues are close enough to allow deflation.
 578:                  // *
 579:                  S = Z[JLAM + o_z];
 580:                  C = Z[J + o_z];
 581:                  // *
 582:                  // *        Find sqrt(a**2+b**2) without overflow or
 583:                  // *        destructive underflow.
 584:                  // *
 585:                  TAU = this._dlapy2.Run(C, S);
 586:                  T = D[J + o_d] - D[JLAM + o_d];
 587:                  C = C / TAU;
 588:                  S =  - S / TAU;
 589:                  if (Math.Abs(T * C * S) <= TOL)
 590:                  {
 591:                      // *
 592:                      // *           Deflation is possible.
 593:                      // *
 594:                      Z[J + o_z] = TAU;
 595:                      Z[JLAM + o_z] = ZERO;
 596:                      // *
 597:                      // *           Record the appropriate Givens rotation
 598:                      // *
 599:                      GIVPTR = GIVPTR + 1;
 600:                      GIVCOL[1+GIVPTR * 2 + o_givcol] = INDXQ[INDX[JLAM + o_indx] + o_indxq];
 601:                      GIVCOL[2+GIVPTR * 2 + o_givcol] = INDXQ[INDX[J + o_indx] + o_indxq];
 602:                      GIVNUM[1+GIVPTR * 2 + o_givnum] = C;
 603:                      GIVNUM[2+GIVPTR * 2 + o_givnum] = S;
 604:                      if (ICOMPQ == 1)
 605:                      {
 606:                          this._drot.Run(QSIZ, ref Q, 1+(INDXQ[INDX[JLAM + o_indx] + o_indxq]) * LDQ + o_q, 1, ref Q, 1+(INDXQ[INDX[J + o_indx] + o_indxq]) * LDQ + o_q, 1, C
 607:                                         , S);
 608:                      }
 609:                      T = D[JLAM + o_d] * C * C + D[J + o_d] * S * S;
 610:                      D[J + o_d] = D[JLAM + o_d] * S * S + D[J + o_d] * C * C;
 611:                      D[JLAM + o_d] = T;
 612:                      K2 = K2 - 1;
 613:                      I = 1;
 614:                  LABEL90:;
 615:                      if (K2 + I <= N)
 616:                      {
 617:                          if (D[JLAM + o_d] < D[INDXP[K2 + I + o_indxp] + o_d])
 618:                          {
 619:                              INDXP[K2 + I - 1 + o_indxp] = INDXP[K2 + I + o_indxp];
 620:                              INDXP[K2 + I + o_indxp] = JLAM;
 621:                              I = I + 1;
 622:                              goto LABEL90;
 623:                          }
 624:                          else
 625:                          {
 626:                              INDXP[K2 + I - 1 + o_indxp] = JLAM;
 627:                          }
 628:                      }
 629:                      else
 630:                      {
 631:                          INDXP[K2 + I - 1 + o_indxp] = JLAM;
 632:                      }
 633:                      JLAM = J;
 634:                  }
 635:                  else
 636:                  {
 637:                      K = K + 1;
 638:                      W[K + o_w] = Z[JLAM + o_z];
 639:                      DLAMDA[K + o_dlamda] = D[JLAM + o_d];
 640:                      INDXP[K + o_indxp] = JLAM;
 641:                      JLAM = J;
 642:                  }
 643:              }
 644:              goto LABEL80;
 645:          LABEL100:;
 646:              // *
 647:              // *     Record the last eigenvalue.
 648:              // *
 649:              K = K + 1;
 650:              W[K + o_w] = Z[JLAM + o_z];
 651:              DLAMDA[K + o_dlamda] = D[JLAM + o_d];
 652:              INDXP[K + o_indxp] = JLAM;
 653:              // *
 654:          LABEL110:;
 655:              // *
 656:              // *     Sort the eigenvalues and corresponding eigenvectors into DLAMDA
 657:              // *     and Q2 respectively.  The eigenvalues/vectors which were not
 658:              // *     deflated go into the first K slots of DLAMDA and Q2 respectively,
 659:              // *     while those which were deflated go into the last N - K slots.
 660:              // *
 661:              if (ICOMPQ == 0)
 662:              {
 663:                  for (J = 1; J <= N; J++)
 664:                  {
 665:                      JP = INDXP[J + o_indxp];
 666:                      DLAMDA[J + o_dlamda] = D[JP + o_d];
 667:                      PERM[J + o_perm] = INDXQ[INDX[JP + o_indx] + o_indxq];
 668:                  }
 669:              }
 670:              else
 671:              {
 672:                  for (J = 1; J <= N; J++)
 673:                  {
 674:                      JP = INDXP[J + o_indxp];
 675:                      DLAMDA[J + o_dlamda] = D[JP + o_d];
 676:                      PERM[J + o_perm] = INDXQ[INDX[JP + o_indx] + o_indxq];
 677:                      this._dcopy.Run(QSIZ, Q, 1+(PERM[J + o_perm]) * LDQ + o_q, 1, ref Q2, 1+J * LDQ2 + o_q2, 1);
 678:                  }
 679:              }
 680:              // *
 681:              // *     The deflated eigenvalues and their corresponding vectors go back
 682:              // *     into the last N - K slots of D and Q respectively.
 683:              // *
 684:              if (K < N)
 685:              {
 686:                  if (ICOMPQ == 0)
 687:                  {
 688:                      this._dcopy.Run(N - K, DLAMDA, K + 1 + o_dlamda, 1, ref D, K + 1 + o_d, 1);
 689:                  }
 690:                  else
 691:                  {
 692:                      this._dcopy.Run(N - K, DLAMDA, K + 1 + o_dlamda, 1, ref D, K + 1 + o_d, 1);
 693:                      this._dlacpy.Run("A", QSIZ, N - K, Q2, 1+(K + 1) * LDQ2 + o_q2, LDQ2, ref Q, 1+(K + 1) * LDQ + o_q
 694:                                       , LDQ);
 695:                  }
 696:              }
 697:              // *
 698:              return;
 699:              // *
 700:              // *     End of DLAED8
 701:              // *
 702:   
 703:              #endregion
 704:   
 705:          }
 706:      }
 707:  }