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