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:      /// DLAED0 computes all eigenvalues and corresponding eigenvectors of a
  27:      /// symmetric tridiagonal matrix using the divide and conquer method.
  28:      /// 
  29:      ///</summary>
  30:      public class DLAED0
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLAED1 _dlaed1; DLAED7 _dlaed7; DSTEQR _dsteqr; XERBLA _xerbla; 
  37:          ILAENV _ilaenv;
  38:   
  39:          #endregion
  40:   
  41:   
  42:          #region Fields
  43:          
  44:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; int CURLVL = 0; int CURPRB = 0; 
  45:          int CURR = 0;int I = 0; int IGIVCL = 0; int IGIVNM = 0; int IGIVPT = 0; int INDXQ = 0; int IPERM = 0; int IPRMPT = 0; 
  46:          int IQ = 0;int IQPTR = 0; int IWREM = 0; int J = 0; int K = 0; int LGN = 0; int MATSIZ = 0; int MSD2 = 0; int SMLSIZ = 0; 
  47:          int SMM1 = 0;int SPM1 = 0; int SPM2 = 0; int SUBMAT = 0; int SUBPBS = 0; int TLVLS = 0; double TEMP = 0; 
  48:   
  49:          #endregion
  50:   
  51:          public DLAED0(DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLAED1 dlaed1, DLAED7 dlaed7, DSTEQR dsteqr, XERBLA xerbla, ILAENV ilaenv)
  52:          {
  53:      
  54:   
  55:              #region Set Dependencies
  56:              
  57:              this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy; this._dlaed1 = dlaed1; this._dlaed7 = dlaed7; 
  58:              this._dsteqr = dsteqr;this._xerbla = xerbla; this._ilaenv = ilaenv; 
  59:   
  60:              #endregion
  61:   
  62:          }
  63:      
  64:          public DLAED0()
  65:          {
  66:      
  67:   
  68:              #region Dependencies (Initialization)
  69:              
  70:              DCOPY dcopy = new DCOPY();
  71:              LSAME lsame = new LSAME();
  72:              XERBLA xerbla = new XERBLA();
  73:              IDAMAX idamax = new IDAMAX();
  74:              DLAMC3 dlamc3 = new DLAMC3();
  75:              DLAPY2 dlapy2 = new DLAPY2();
  76:              DLAMRG dlamrg = new DLAMRG();
  77:              DROT drot = new DROT();
  78:              DSCAL dscal = new DSCAL();
  79:              DNRM2 dnrm2 = new DNRM2();
  80:              DLAED5 dlaed5 = new DLAED5();
  81:              DLASSQ dlassq = new DLASSQ();
  82:              DLAE2 dlae2 = new DLAE2();
  83:              DLAEV2 dlaev2 = new DLAEV2();
  84:              DSWAP dswap = new DSWAP();
  85:              IEEECK ieeeck = new IEEECK();
  86:              IPARMQ iparmq = new IPARMQ();
  87:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  88:              DLACPY dlacpy = new DLACPY(lsame);
  89:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  90:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  91:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  92:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  93:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  94:              DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
  95:              DLAED6 dlaed6 = new DLAED6(dlamch);
  96:              DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
  97:              DLASET dlaset = new DLASET(lsame);
  98:              DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);
  99:              DLAED1 dlaed1 = new DLAED1(dcopy, dlaed2, dlaed3, dlamrg, xerbla);
 100:              DLAED8 dlaed8 = new DLAED8(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
 101:              DLAED9 dlaed9 = new DLAED9(dlamc3, dnrm2, dcopy, dlaed4, xerbla);
 102:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 103:              DLAEDA dlaeda = new DLAEDA(dcopy, dgemv, drot, xerbla);
 104:              DLAED7 dlaed7 = new DLAED7(dgemm, dlaed8, dlaed9, dlaeda, dlamrg, xerbla);
 105:              DLANST dlanst = new DLANST(lsame, dlassq);
 106:              DLARTG dlartg = new DLARTG(dlamch);
 107:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 108:              DLASR dlasr = new DLASR(lsame, xerbla);
 109:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
 110:              DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
 111:                                         , dlasrt, dswap, xerbla);
 112:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 113:   
 114:              #endregion
 115:   
 116:   
 117:              #region Set Dependencies
 118:              
 119:              this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy; this._dlaed1 = dlaed1; this._dlaed7 = dlaed7; 
 120:              this._dsteqr = dsteqr;this._xerbla = xerbla; this._ilaenv = ilaenv; 
 121:   
 122:              #endregion
 123:   
 124:          }
 125:          /// <summary>
 126:          /// Purpose
 127:          /// =======
 128:          /// 
 129:          /// DLAED0 computes all eigenvalues and corresponding eigenvectors of a
 130:          /// symmetric tridiagonal matrix using the divide and conquer method.
 131:          /// 
 132:          ///</summary>
 133:          /// <param name="ICOMPQ">
 134:          /// (input) INTEGER
 135:          /// = 0:  Compute eigenvalues only.
 136:          /// = 1:  Compute eigenvectors of original dense symmetric matrix
 137:          /// also.  On entry, Q contains the orthogonal matrix used
 138:          /// to reduce the original matrix to tridiagonal form.
 139:          /// = 2:  Compute eigenvalues and eigenvectors of tridiagonal
 140:          /// matrix.
 141:          ///</param>
 142:          /// <param name="QSIZ">
 143:          /// (input) INTEGER
 144:          /// The dimension of the orthogonal matrix used to reduce
 145:          /// the full matrix to tridiagonal form.  QSIZ .GE. N if ICOMPQ = 1.
 146:          ///</param>
 147:          /// <param name="N">
 148:          /// (input) INTEGER
 149:          /// The dimension of the symmetric tridiagonal matrix.  N .GE. 0.
 150:          ///</param>
 151:          /// <param name="D">
 152:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 153:          /// On entry, the main diagonal of the tridiagonal matrix.
 154:          /// On exit, its eigenvalues.
 155:          ///</param>
 156:          /// <param name="E">
 157:          /// (input) DOUBLE PRECISION array, dimension (N-1)
 158:          /// The off-diagonal elements of the tridiagonal matrix.
 159:          /// On exit, E has been destroyed.
 160:          ///</param>
 161:          /// <param name="Q">
 162:          /// (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
 163:          /// On entry, Q must contain an N-by-N orthogonal matrix.
 164:          /// If ICOMPQ = 0    Q is not referenced.
 165:          /// If ICOMPQ = 1    On entry, Q is a subset of the columns of the
 166:          /// orthogonal matrix used to reduce the full
 167:          /// matrix to tridiagonal form corresponding to
 168:          /// the subset of the full matrix which is being
 169:          /// decomposed at this time.
 170:          /// If ICOMPQ = 2    On entry, Q will be the identity matrix.
 171:          /// On exit, Q contains the eigenvectors of the
 172:          /// tridiagonal matrix.
 173:          ///</param>
 174:          /// <param name="LDQ">
 175:          /// (input) INTEGER
 176:          /// The leading dimension of the array Q.  If eigenvectors are
 177:          /// desired, then  LDQ .GE. max(1,N).  In any case,  LDQ .GE. 1.
 178:          ///</param>
 179:          /// <param name="QSTORE">
 180:          /// (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
 181:          /// Referenced only when ICOMPQ = 1.  Used to store parts of
 182:          /// the eigenvector matrix when the updating matrix multiplies
 183:          /// take place.
 184:          ///</param>
 185:          /// <param name="LDQS">
 186:          /// (input) INTEGER
 187:          /// The leading dimension of the array QSTORE.  If ICOMPQ = 1,
 188:          /// then  LDQS .GE. max(1,N).  In any case,  LDQS .GE. 1.
 189:          ///</param>
 190:          /// <param name="WORK">
 191:          /// (workspace) DOUBLE PRECISION array,
 192:          /// If ICOMPQ = 0 or 1, the dimension of WORK must be at least
 193:          /// 1 + 3*N + 2*N*lg N + 2*N**2
 194:          /// ( lg( N ) = smallest integer k
 195:          /// such that 2^k .GE. N )
 196:          /// If ICOMPQ = 2, the dimension of WORK must be at least
 197:          /// 4*N + N**2.
 198:          ///</param>
 199:          /// <param name="IWORK">
 200:          /// (workspace) INTEGER array,
 201:          /// If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
 202:          /// 6 + 6*N + 5*N*lg N.
 203:          /// ( lg( N ) = smallest integer k
 204:          /// such that 2^k .GE. N )
 205:          /// If ICOMPQ = 2, the dimension of IWORK must be at least
 206:          /// 3 + 5*N.
 207:          ///</param>
 208:          /// <param name="INFO">
 209:          /// (output) INTEGER
 210:          /// = 0:  successful exit.
 211:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 212:          /// .GT. 0:  The algorithm failed to compute an eigenvalue while
 213:          /// working on the submatrix lying in rows and columns
 214:          /// INFO/(N+1) through mod(INFO,N+1).
 215:          ///</param>
 216:          public void Run(int ICOMPQ, int QSIZ, int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] Q, int offset_q
 217:                           , int LDQ, ref double[] QSTORE, int offset_qstore, int LDQS, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
 218:          {
 219:   
 220:              #region Array Index Correction
 221:              
 222:               int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_q = -1 - LDQ + offset_q; 
 223:               int o_qstore = -1 - LDQS + offset_qstore; int o_work = -1 + offset_work;  int o_iwork = -1 + offset_iwork; 
 224:   
 225:              #endregion
 226:   
 227:   
 228:              #region Prolog
 229:              
 230:              // *
 231:              // *  -- LAPACK routine (version 3.1) --
 232:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 233:              // *     November 2006
 234:              // *
 235:              // *     .. Scalar Arguments ..
 236:              // *     ..
 237:              // *     .. Array Arguments ..
 238:              // *     ..
 239:              // *
 240:              // *  Purpose
 241:              // *  =======
 242:              // *
 243:              // *  DLAED0 computes all eigenvalues and corresponding eigenvectors of a
 244:              // *  symmetric tridiagonal matrix using the divide and conquer method.
 245:              // *
 246:              // *  Arguments
 247:              // *  =========
 248:              // *
 249:              // *  ICOMPQ  (input) INTEGER
 250:              // *          = 0:  Compute eigenvalues only.
 251:              // *          = 1:  Compute eigenvectors of original dense symmetric matrix
 252:              // *                also.  On entry, Q contains the orthogonal matrix used
 253:              // *                to reduce the original matrix to tridiagonal form.
 254:              // *          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
 255:              // *                matrix.
 256:              // *
 257:              // *  QSIZ   (input) INTEGER
 258:              // *         The dimension of the orthogonal matrix used to reduce
 259:              // *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
 260:              // *
 261:              // *  N      (input) INTEGER
 262:              // *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
 263:              // *
 264:              // *  D      (input/output) DOUBLE PRECISION array, dimension (N)
 265:              // *         On entry, the main diagonal of the tridiagonal matrix.
 266:              // *         On exit, its eigenvalues.
 267:              // *
 268:              // *  E      (input) DOUBLE PRECISION array, dimension (N-1)
 269:              // *         The off-diagonal elements of the tridiagonal matrix.
 270:              // *         On exit, E has been destroyed.
 271:              // *
 272:              // *  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
 273:              // *         On entry, Q must contain an N-by-N orthogonal matrix.
 274:              // *         If ICOMPQ = 0    Q is not referenced.
 275:              // *         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
 276:              // *                          orthogonal matrix used to reduce the full
 277:              // *                          matrix to tridiagonal form corresponding to
 278:              // *                          the subset of the full matrix which is being
 279:              // *                          decomposed at this time.
 280:              // *         If ICOMPQ = 2    On entry, Q will be the identity matrix.
 281:              // *                          On exit, Q contains the eigenvectors of the
 282:              // *                          tridiagonal matrix.
 283:              // *
 284:              // *  LDQ    (input) INTEGER
 285:              // *         The leading dimension of the array Q.  If eigenvectors are
 286:              // *         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
 287:              // *
 288:              // *  QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N)
 289:              // *         Referenced only when ICOMPQ = 1.  Used to store parts of
 290:              // *         the eigenvector matrix when the updating matrix multiplies
 291:              // *         take place.
 292:              // *
 293:              // *  LDQS   (input) INTEGER
 294:              // *         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
 295:              // *         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
 296:              // *
 297:              // *  WORK   (workspace) DOUBLE PRECISION array,
 298:              // *         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
 299:              // *                     1 + 3*N + 2*N*lg N + 2*N**2
 300:              // *                     ( lg( N ) = smallest integer k
 301:              // *                                 such that 2^k >= N )
 302:              // *         If ICOMPQ = 2, the dimension of WORK must be at least
 303:              // *                     4*N + N**2.
 304:              // *
 305:              // *  IWORK  (workspace) INTEGER array,
 306:              // *         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
 307:              // *                        6 + 6*N + 5*N*lg N.
 308:              // *                        ( lg( N ) = smallest integer k
 309:              // *                                    such that 2^k >= N )
 310:              // *         If ICOMPQ = 2, the dimension of IWORK must be at least
 311:              // *                        3 + 5*N.
 312:              // *
 313:              // *  INFO   (output) INTEGER
 314:              // *          = 0:  successful exit.
 315:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 316:              // *          > 0:  The algorithm failed to compute an eigenvalue while
 317:              // *                working on the submatrix lying in rows and columns
 318:              // *                INFO/(N+1) through mod(INFO,N+1).
 319:              // *
 320:              // *  Further Details
 321:              // *  ===============
 322:              // *
 323:              // *  Based on contributions by
 324:              // *     Jeff Rutter, Computer Science Division, University of California
 325:              // *     at Berkeley, USA
 326:              // *
 327:              // *  =====================================================================
 328:              // *
 329:              // *     .. Parameters ..
 330:              // *     ..
 331:              // *     .. Local Scalars ..
 332:              // *     ..
 333:              // *     .. External Subroutines ..
 334:              // *     ..
 335:              // *     .. External Functions ..
 336:              // *     ..
 337:              // *     .. Intrinsic Functions ..
 338:              //      INTRINSIC          ABS, DBLE, INT, LOG, MAX;
 339:              // *     ..
 340:              // *     .. Executable Statements ..
 341:              // *
 342:              // *     Test the input parameters.
 343:              // *
 344:   
 345:              #endregion
 346:   
 347:   
 348:              #region Body
 349:              
 350:              INFO = 0;
 351:              // *
 352:              if (ICOMPQ < 0 || ICOMPQ > 2)
 353:              {
 354:                  INFO =  - 1;
 355:              }
 356:              else
 357:              {
 358:                  if ((ICOMPQ == 1) && (QSIZ < Math.Max(0, N)))
 359:                  {
 360:                      INFO =  - 2;
 361:                  }
 362:                  else
 363:                  {
 364:                      if (N < 0)
 365:                      {
 366:                          INFO =  - 3;
 367:                      }
 368:                      else
 369:                      {
 370:                          if (LDQ < Math.Max(1, N))
 371:                          {
 372:                              INFO =  - 7;
 373:                          }
 374:                          else
 375:                          {
 376:                              if (LDQS < Math.Max(1, N))
 377:                              {
 378:                                  INFO =  - 9;
 379:                              }
 380:                          }
 381:                      }
 382:                  }
 383:              }
 384:              if (INFO != 0)
 385:              {
 386:                  this._xerbla.Run("DLAED0",  - INFO);
 387:                  return;
 388:              }
 389:              // *
 390:              // *     Quick return if possible
 391:              // *
 392:              if (N == 0) return;
 393:              // *
 394:              SMLSIZ = this._ilaenv.Run(9, "DLAED0", " ", 0, 0, 0, 0);
 395:              // *
 396:              // *     Determine the size and placement of the submatrices, and save in
 397:              // *     the leading elements of IWORK.
 398:              // *
 399:              IWORK[1 + o_iwork] = N;
 400:              SUBPBS = 1;
 401:              TLVLS = 0;
 402:          LABEL10:;
 403:              if (IWORK[SUBPBS + o_iwork] > SMLSIZ)
 404:              {
 405:                  for (J = SUBPBS; J >= 1; J +=  - 1)
 406:                  {
 407:                      IWORK[2 * J + o_iwork] = (IWORK[J + o_iwork] + 1) / 2;
 408:                      IWORK[2 * J - 1 + o_iwork] = IWORK[J + o_iwork] / 2;
 409:                  }
 410:                  TLVLS = TLVLS + 1;
 411:                  SUBPBS = 2 * SUBPBS;
 412:                  goto LABEL10;
 413:              }
 414:              for (J = 2; J <= SUBPBS; J++)
 415:              {
 416:                  IWORK[J + o_iwork] = IWORK[J + o_iwork] + IWORK[J - 1 + o_iwork];
 417:              }
 418:              // *
 419:              // *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
 420:              // *     using rank-1 modifications (cuts).
 421:              // *
 422:              SPM1 = SUBPBS - 1;
 423:              for (I = 1; I <= SPM1; I++)
 424:              {
 425:                  SUBMAT = IWORK[I + o_iwork] + 1;
 426:                  SMM1 = SUBMAT - 1;
 427:                  D[SMM1 + o_d] = D[SMM1 + o_d] - Math.Abs(E[SMM1 + o_e]);
 428:                  D[SUBMAT + o_d] = D[SUBMAT + o_d] - Math.Abs(E[SMM1 + o_e]);
 429:              }
 430:              // *
 431:              INDXQ = 4 * N + 3;
 432:              if (ICOMPQ != 2)
 433:              {
 434:                  // *
 435:                  // *        Set up workspaces for eigenvalues only/accumulate new vectors
 436:                  // *        routine
 437:                  // *
 438:                  TEMP = Math.Log(Convert.ToDouble(N)) / Math.Log(TWO);
 439:                  LGN = Convert.ToInt32(Math.Truncate(TEMP));
 440:                  if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
 441:                  if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
 442:                  IPRMPT = INDXQ + N + 1;
 443:                  IPERM = IPRMPT + N * LGN;
 444:                  IQPTR = IPERM + N * LGN;
 445:                  IGIVPT = IQPTR + N + 2;
 446:                  IGIVCL = IGIVPT + N * LGN;
 447:                  // *
 448:                  IGIVNM = 1;
 449:                  IQ = IGIVNM + 2 * N * LGN;
 450:                  IWREM = IQ + (int)Math.Pow(N, 2) + 1;
 451:                  // *
 452:                  // *        Initialize pointers
 453:                  // *
 454:                  for (I = 0; I <= SUBPBS; I++)
 455:                  {
 456:                      IWORK[IPRMPT + I + o_iwork] = 1;
 457:                      IWORK[IGIVPT + I + o_iwork] = 1;
 458:                  }
 459:                  IWORK[IQPTR + o_iwork] = 1;
 460:              }
 461:              // *
 462:              // *     Solve each submatrix eigenproblem at the bottom of the divide and
 463:              // *     conquer tree.
 464:              // *
 465:              CURR = 0;
 466:              for (I = 0; I <= SPM1; I++)
 467:              {
 468:                  if (I == 0)
 469:                  {
 470:                      SUBMAT = 1;
 471:                      MATSIZ = IWORK[1 + o_iwork];
 472:                  }
 473:                  else
 474:                  {
 475:                      SUBMAT = IWORK[I + o_iwork] + 1;
 476:                      MATSIZ = IWORK[I + 1 + o_iwork] - IWORK[I + o_iwork];
 477:                  }
 478:                  if (ICOMPQ == 2)
 479:                  {
 480:                      this._dsteqr.Run("I", MATSIZ, ref D, SUBMAT + o_d, ref E, SUBMAT + o_e, ref Q, SUBMAT+SUBMAT * LDQ + o_q, LDQ
 481:                                       , ref WORK, offset_work, ref INFO);
 482:                      if (INFO != 0) goto LABEL130;
 483:                  }
 484:                  else
 485:                  {
 486:                      this._dsteqr.Run("I", MATSIZ, ref D, SUBMAT + o_d, ref E, SUBMAT + o_e, ref WORK, IQ - 1 + IWORK[IQPTR + CURR + o_iwork] + o_work, MATSIZ
 487:                                       , ref WORK, offset_work, ref INFO);
 488:                      if (INFO != 0) goto LABEL130;
 489:                      if (ICOMPQ == 1)
 490:                      {
 491:                          this._dgemm.Run("N", "N", QSIZ, MATSIZ, MATSIZ, ONE
 492:                                          , Q, 1+SUBMAT * LDQ + o_q, LDQ, WORK, IQ - 1 + IWORK[IQPTR + CURR + o_iwork] + o_work, MATSIZ, ZERO, ref QSTORE, 1+SUBMAT * LDQS + o_qstore
 493:                                          , LDQS);
 494:                      }
 495:                      IWORK[IQPTR + CURR + 1 + o_iwork] = IWORK[IQPTR + CURR + o_iwork] + (int)Math.Pow(MATSIZ, 2);
 496:                      CURR = CURR + 1;
 497:                  }
 498:                  K = 1;
 499:                  for (J = SUBMAT; J <= IWORK[I + 1 + o_iwork]; J++)
 500:                  {
 501:                      IWORK[INDXQ + J + o_iwork] = K;
 502:                      K = K + 1;
 503:                  }
 504:              }
 505:              // *
 506:              // *     Successively merge eigensystems of adjacent submatrices
 507:              // *     into eigensystem for the corresponding larger matrix.
 508:              // *
 509:              // *     while ( SUBPBS > 1 )
 510:              // *
 511:              CURLVL = 1;
 512:          LABEL80:;
 513:              if (SUBPBS > 1)
 514:              {
 515:                  SPM2 = SUBPBS - 2;
 516:                  for (I = 0; I <= SPM2; I += 2)
 517:                  {
 518:                      if (I == 0)
 519:                      {
 520:                          SUBMAT = 1;
 521:                          MATSIZ = IWORK[2 + o_iwork];
 522:                          MSD2 = IWORK[1 + o_iwork];
 523:                          CURPRB = 0;
 524:                      }
 525:                      else
 526:                      {
 527:                          SUBMAT = IWORK[I + o_iwork] + 1;
 528:                          MATSIZ = IWORK[I + 2 + o_iwork] - IWORK[I + o_iwork];
 529:                          MSD2 = MATSIZ / 2;
 530:                          CURPRB = CURPRB + 1;
 531:                      }
 532:                      // *
 533:                      // *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
 534:                      // *     into an eigensystem of size MATSIZ.
 535:                      // *     DLAED1 is used only for the full eigensystem of a tridiagonal
 536:                      // *     matrix.
 537:                      // *     DLAED7 handles the cases in which eigenvalues only or eigenvalues
 538:                      // *     and eigenvectors of a full symmetric matrix (which was reduced to
 539:                      // *     tridiagonal form) are desired.
 540:                      // *
 541:                      if (ICOMPQ == 2)
 542:                      {
 543:                          this._dlaed1.Run(MATSIZ, ref D, SUBMAT + o_d, ref Q, SUBMAT+SUBMAT * LDQ + o_q, LDQ, ref IWORK, INDXQ + SUBMAT + o_iwork, ref E[SUBMAT + MSD2 - 1 + o_e]
 544:                                           , MSD2, ref WORK, offset_work, ref IWORK, SUBPBS + 1 + o_iwork, ref INFO);
 545:                      }
 546:                      else
 547:                      {
 548:                          this._dlaed7.Run(ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB
 549:                                           , ref D, SUBMAT + o_d, ref QSTORE, 1+SUBMAT * LDQS + o_qstore, LDQS, ref IWORK, INDXQ + SUBMAT + o_iwork, ref E[SUBMAT + MSD2 - 1 + o_e], MSD2
 550:                                           , ref WORK, IQ + o_work, ref IWORK, IQPTR + o_iwork, ref IWORK, IPRMPT + o_iwork, ref IWORK, IPERM + o_iwork, ref IWORK, IGIVPT + o_iwork, ref IWORK, IGIVCL + o_iwork
 551:                                           , ref WORK, IGIVNM + o_work, ref WORK, IWREM + o_work, ref IWORK, SUBPBS + 1 + o_iwork, ref INFO);
 552:                      }
 553:                      if (INFO != 0) goto LABEL130;
 554:                      IWORK[I / 2 + 1 + o_iwork] = IWORK[I + 2 + o_iwork];
 555:                  }
 556:                  SUBPBS = SUBPBS / 2;
 557:                  CURLVL = CURLVL + 1;
 558:                  goto LABEL80;
 559:              }
 560:              // *
 561:              // *     end while
 562:              // *
 563:              // *     Re-merge the eigenvalues/vectors which were deflated at the final
 564:              // *     merge step.
 565:              // *
 566:              if (ICOMPQ == 1)
 567:              {
 568:                  for (I = 1; I <= N; I++)
 569:                  {
 570:                      J = IWORK[INDXQ + I + o_iwork];
 571:                      WORK[I + o_work] = D[J + o_d];
 572:                      this._dcopy.Run(QSIZ, QSTORE, 1+J * LDQS + o_qstore, 1, ref Q, 1+I * LDQ + o_q, 1);
 573:                  }
 574:                  this._dcopy.Run(N, WORK, offset_work, 1, ref D, offset_d, 1);
 575:              }
 576:              else
 577:              {
 578:                  if (ICOMPQ == 2)
 579:                  {
 580:                      for (I = 1; I <= N; I++)
 581:                      {
 582:                          J = IWORK[INDXQ + I + o_iwork];
 583:                          WORK[I + o_work] = D[J + o_d];
 584:                          this._dcopy.Run(N, Q, 1+J * LDQ + o_q, 1, ref WORK, N * I + 1 + o_work, 1);
 585:                      }
 586:                      this._dcopy.Run(N, WORK, offset_work, 1, ref D, offset_d, 1);
 587:                      this._dlacpy.Run("A", N, N, WORK, N + 1 + o_work, N, ref Q, offset_q
 588:                                       , LDQ);
 589:                  }
 590:                  else
 591:                  {
 592:                      for (I = 1; I <= N; I++)
 593:                      {
 594:                          J = IWORK[INDXQ + I + o_iwork];
 595:                          WORK[I + o_work] = D[J + o_d];
 596:                      }
 597:                      this._dcopy.Run(N, WORK, offset_work, 1, ref D, offset_d, 1);
 598:                  }
 599:              }
 600:              goto LABEL140;
 601:              // *
 602:          LABEL130:;
 603:              INFO = SUBMAT * (N + 1) + SUBMAT + MATSIZ - 1;
 604:              // *
 605:          LABEL140:;
 606:              return;
 607:              // *
 608:              // *     End of DLAED0
 609:              // *
 610:   
 611:              #endregion
 612:   
 613:          }
 614:      }
 615:  }