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:      /// DLALSD uses the singular value decomposition of A to solve the least
  27:      /// squares problem of finding X to minimize the Euclidean norm of each
  28:      /// column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
  29:      /// are N-by-NRHS. The solution X overwrites B.
  30:      /// 
  31:      /// The singular values of A smaller than RCOND times the largest
  32:      /// singular value are treated as zero in solving the least squares
  33:      /// problem; in this case a minimum norm solution is returned.
  34:      /// The actual singular values are returned in D in ascending order.
  35:      /// 
  36:      /// This code makes very mild assumptions about floating point
  37:      /// arithmetic. It will work on machines with a guard digit in
  38:      /// add/subtract, or on those binary machines without guard digits
  39:      /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
  40:      /// It could conceivably fail on hexadecimal or decimal machines
  41:      /// without guard digits, but we know of none.
  42:      /// 
  43:      ///</summary>
  44:      public class DLALSD
  45:      {
  46:      
  47:   
  48:          #region Dependencies
  49:          
  50:          IDAMAX _idamax; DLAMCH _dlamch; DLANST _dlanst; DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLALSA _dlalsa; 
  51:          DLARTG _dlartg;DLASCL _dlascl; DLASDA _dlasda; DLASDQ _dlasdq; DLASET _dlaset; DLASRT _dlasrt; DROT _drot; XERBLA _xerbla; 
  52:   
  53:          #endregion
  54:   
  55:   
  56:          #region Fields
  57:          
  58:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; int BX = 0; int BXST = 0; int C = 0; 
  59:          int DIFL = 0;int DIFR = 0; int GIVCOL = 0; int GIVNUM = 0; int GIVPTR = 0; int I = 0; int ICMPQ1 = 0; int ICMPQ2 = 0; 
  60:          int IWK = 0;int J = 0; int K = 0; int NLVL = 0; int NM1 = 0; int NSIZE = 0; int NSUB = 0; int NWORK = 0; int PERM = 0; 
  61:          int POLES = 0;int S = 0; int SIZEI = 0; int SMLSZP = 0; int SQRE = 0; int ST = 0; int ST1 = 0; int U = 0; int VT = 0; 
  62:          int Z = 0;double CS = 0; double EPS = 0; double ORGNRM = 0; double R = 0; double RCND = 0; double SN = 0; double TOL = 0; 
  63:   
  64:          #endregion
  65:   
  66:          public DLALSD(IDAMAX idamax, DLAMCH dlamch, DLANST dlanst, DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLALSA dlalsa, DLARTG dlartg, DLASCL dlascl, DLASDA dlasda
  67:                        , DLASDQ dlasdq, DLASET dlaset, DLASRT dlasrt, DROT drot, XERBLA xerbla)
  68:          {
  69:      
  70:   
  71:              #region Set Dependencies
  72:              
  73:              this._idamax = idamax; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy; this._dgemm = dgemm; 
  74:              this._dlacpy = dlacpy;this._dlalsa = dlalsa; this._dlartg = dlartg; this._dlascl = dlascl; this._dlasda = dlasda; 
  75:              this._dlasdq = dlasdq;this._dlaset = dlaset; this._dlasrt = dlasrt; this._drot = drot; this._xerbla = xerbla; 
  76:   
  77:              #endregion
  78:   
  79:          }
  80:      
  81:          public DLALSD()
  82:          {
  83:      
  84:   
  85:              #region Dependencies (Initialization)
  86:              
  87:              IDAMAX idamax = new IDAMAX();
  88:              LSAME lsame = new LSAME();
  89:              DLAMC3 dlamc3 = new DLAMC3();
  90:              DLASSQ dlassq = new DLASSQ();
  91:              DCOPY dcopy = new DCOPY();
  92:              XERBLA xerbla = new XERBLA();
  93:              DROT drot = new DROT();
  94:              DSCAL dscal = new DSCAL();
  95:              DNRM2 dnrm2 = new DNRM2();
  96:              DLASDT dlasdt = new DLASDT();
  97:              DLAMRG dlamrg = new DLAMRG();
  98:              DLAPY2 dlapy2 = new DLAPY2();
  99:              DLASD5 dlasd5 = new DLASD5();
 100:              DDOT ddot = new DDOT();
 101:              DLAS2 dlas2 = new DLAS2();
 102:              DLASQ5 dlasq5 = new DLASQ5();
 103:              DLAZQ4 dlazq4 = new DLAZQ4();
 104:              IEEECK ieeeck = new IEEECK();
 105:              IPARMQ iparmq = new IPARMQ();
 106:              DSWAP dswap = new DSWAP();
 107:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 108:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 109:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 110:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 111:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 112:              DLANST dlanst = new DLANST(lsame, dlassq);
 113:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 114:              DLACPY dlacpy = new DLACPY(lsame);
 115:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 116:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 117:              DLALS0 dlals0 = new DLALS0(dcopy, dgemv, dlacpy, dlascl, drot, dscal, xerbla, dlamc3, dnrm2);
 118:              DLALSA dlalsa = new DLALSA(dcopy, dgemm, dlals0, dlasdt, xerbla);
 119:              DLARTG dlartg = new DLARTG(dlamch);
 120:              DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);
 121:              DLAED6 dlaed6 = new DLAED6(dlamch);
 122:              DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
 123:              DLASET dlaset = new DLASET(lsame);
 124:              DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);
 125:              DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);
 126:              DLASQ6 dlasq6 = new DLASQ6(dlamch);
 127:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
 128:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
 129:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 130:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
 131:              DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
 132:              DLASR dlasr = new DLASR(lsame, xerbla);
 133:              DLASV2 dlasv2 = new DLASV2(dlamch);
 134:              DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
 135:                                         , xerbla);
 136:              DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);
 137:              DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);
 138:   
 139:              #endregion
 140:   
 141:   
 142:              #region Set Dependencies
 143:              
 144:              this._idamax = idamax; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy; this._dgemm = dgemm; 
 145:              this._dlacpy = dlacpy;this._dlalsa = dlalsa; this._dlartg = dlartg; this._dlascl = dlascl; this._dlasda = dlasda; 
 146:              this._dlasdq = dlasdq;this._dlaset = dlaset; this._dlasrt = dlasrt; this._drot = drot; this._xerbla = xerbla; 
 147:   
 148:              #endregion
 149:   
 150:          }
 151:          /// <summary>
 152:          /// Purpose
 153:          /// =======
 154:          /// 
 155:          /// DLALSD uses the singular value decomposition of A to solve the least
 156:          /// squares problem of finding X to minimize the Euclidean norm of each
 157:          /// column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
 158:          /// are N-by-NRHS. The solution X overwrites B.
 159:          /// 
 160:          /// The singular values of A smaller than RCOND times the largest
 161:          /// singular value are treated as zero in solving the least squares
 162:          /// problem; in this case a minimum norm solution is returned.
 163:          /// The actual singular values are returned in D in ascending order.
 164:          /// 
 165:          /// This code makes very mild assumptions about floating point
 166:          /// arithmetic. It will work on machines with a guard digit in
 167:          /// add/subtract, or on those binary machines without guard digits
 168:          /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 169:          /// It could conceivably fail on hexadecimal or decimal machines
 170:          /// without guard digits, but we know of none.
 171:          /// 
 172:          ///</summary>
 173:          /// <param name="UPLO">
 174:          /// (input) CHARACTER*1
 175:          /// = 'U': D and E define an upper bidiagonal matrix.
 176:          /// = 'L': D and E define a  lower bidiagonal matrix.
 177:          ///</param>
 178:          /// <param name="SMLSIZ">
 179:          /// (input) INTEGER
 180:          /// The maximum size of the subproblems at the bottom of the
 181:          /// computation tree.
 182:          ///</param>
 183:          /// <param name="N">
 184:          /// (input) INTEGER
 185:          /// The dimension of the  bidiagonal matrix.  N .GE. 0.
 186:          ///</param>
 187:          /// <param name="NRHS">
 188:          /// (input) INTEGER
 189:          /// The number of columns of B. NRHS must be at least 1.
 190:          ///</param>
 191:          /// <param name="D">
 192:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 193:          /// On entry D contains the main diagonal of the bidiagonal
 194:          /// matrix. On exit, if INFO = 0, D contains its singular values.
 195:          ///</param>
 196:          /// <param name="E">
 197:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)
 198:          /// Contains the super-diagonal entries of the bidiagonal matrix.
 199:          /// On exit, E has been destroyed.
 200:          ///</param>
 201:          /// <param name="B">
 202:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 203:          /// On input, B contains the right hand sides of the least
 204:          /// squares problem. On output, B contains the solution X.
 205:          ///</param>
 206:          /// <param name="LDB">
 207:          /// (input) INTEGER
 208:          /// The leading dimension of B in the calling subprogram.
 209:          /// LDB must be at least max(1,N).
 210:          ///</param>
 211:          /// <param name="RCOND">
 212:          /// (input) DOUBLE PRECISION
 213:          /// The singular values of A less than or equal to RCOND times
 214:          /// the largest singular value are treated as zero in solving
 215:          /// the least squares problem. If RCOND is negative,
 216:          /// machine precision is used instead.
 217:          /// For example, if diag(S)*X=B were the least squares problem,
 218:          /// where diag(S) is a diagonal matrix of singular values, the
 219:          /// solution would be X(i) = B(i) / S(i) if S(i) is greater than
 220:          /// RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
 221:          /// RCOND*max(S).
 222:          ///</param>
 223:          /// <param name="RANK">
 224:          /// (output) INTEGER
 225:          /// The number of singular values of A greater than RCOND times
 226:          /// the largest singular value.
 227:          ///</param>
 228:          /// <param name="WORK">
 229:          /// (workspace) DOUBLE PRECISION array, dimension at least
 230:          /// (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
 231:          /// where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
 232:          ///</param>
 233:          /// <param name="IWORK">
 234:          /// (workspace) INTEGER array, dimension at least
 235:          /// (3*N*NLVL + 11*N)
 236:          ///</param>
 237:          /// <param name="INFO">
 238:          /// (output) INTEGER
 239:          /// = 0:  successful exit.
 240:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 241:          /// .GT. 0:  The algorithm failed to compute an singular value while
 242:          /// working on the submatrix lying in rows and columns
 243:          /// INFO/(N+1) through MOD(INFO,N+1).
 244:          ///</param>
 245:          public void Run(string UPLO, int SMLSIZ, int N, int NRHS, ref double[] D, int offset_d, ref double[] E, int offset_e
 246:                           , ref double[] B, int offset_b, int LDB, double RCOND, ref int RANK, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork
 247:                           , ref int INFO)
 248:          {
 249:   
 250:              #region Array Index Correction
 251:              
 252:               int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_b = -1 - LDB + offset_b;  int o_work = -1 + offset_work; 
 253:               int o_iwork = -1 + offset_iwork;
 254:   
 255:              #endregion
 256:   
 257:   
 258:              #region Strings
 259:              
 260:              UPLO = UPLO.Substring(0, 1);  
 261:   
 262:              #endregion
 263:   
 264:   
 265:              #region Prolog
 266:              
 267:              // *
 268:              // *  -- LAPACK routine (version 3.1) --
 269:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 270:              // *     November 2006
 271:              // *
 272:              // *     .. Scalar Arguments ..
 273:              // *     ..
 274:              // *     .. Array Arguments ..
 275:              // *     ..
 276:              // *
 277:              // *  Purpose
 278:              // *  =======
 279:              // *
 280:              // *  DLALSD uses the singular value decomposition of A to solve the least
 281:              // *  squares problem of finding X to minimize the Euclidean norm of each
 282:              // *  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
 283:              // *  are N-by-NRHS. The solution X overwrites B.
 284:              // *
 285:              // *  The singular values of A smaller than RCOND times the largest
 286:              // *  singular value are treated as zero in solving the least squares
 287:              // *  problem; in this case a minimum norm solution is returned.
 288:              // *  The actual singular values are returned in D in ascending order.
 289:              // *
 290:              // *  This code makes very mild assumptions about floating point
 291:              // *  arithmetic. It will work on machines with a guard digit in
 292:              // *  add/subtract, or on those binary machines without guard digits
 293:              // *  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 294:              // *  It could conceivably fail on hexadecimal or decimal machines
 295:              // *  without guard digits, but we know of none.
 296:              // *
 297:              // *  Arguments
 298:              // *  =========
 299:              // *
 300:              // *  UPLO   (input) CHARACTER*1
 301:              // *         = 'U': D and E define an upper bidiagonal matrix.
 302:              // *         = 'L': D and E define a  lower bidiagonal matrix.
 303:              // *
 304:              // *  SMLSIZ (input) INTEGER
 305:              // *         The maximum size of the subproblems at the bottom of the
 306:              // *         computation tree.
 307:              // *
 308:              // *  N      (input) INTEGER
 309:              // *         The dimension of the  bidiagonal matrix.  N >= 0.
 310:              // *
 311:              // *  NRHS   (input) INTEGER
 312:              // *         The number of columns of B. NRHS must be at least 1.
 313:              // *
 314:              // *  D      (input/output) DOUBLE PRECISION array, dimension (N)
 315:              // *         On entry D contains the main diagonal of the bidiagonal
 316:              // *         matrix. On exit, if INFO = 0, D contains its singular values.
 317:              // *
 318:              // *  E      (input/output) DOUBLE PRECISION array, dimension (N-1)
 319:              // *         Contains the super-diagonal entries of the bidiagonal matrix.
 320:              // *         On exit, E has been destroyed.
 321:              // *
 322:              // *  B      (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 323:              // *         On input, B contains the right hand sides of the least
 324:              // *         squares problem. On output, B contains the solution X.
 325:              // *
 326:              // *  LDB    (input) INTEGER
 327:              // *         The leading dimension of B in the calling subprogram.
 328:              // *         LDB must be at least max(1,N).
 329:              // *
 330:              // *  RCOND  (input) DOUBLE PRECISION
 331:              // *         The singular values of A less than or equal to RCOND times
 332:              // *         the largest singular value are treated as zero in solving
 333:              // *         the least squares problem. If RCOND is negative,
 334:              // *         machine precision is used instead.
 335:              // *         For example, if diag(S)*X=B were the least squares problem,
 336:              // *         where diag(S) is a diagonal matrix of singular values, the
 337:              // *         solution would be X(i) = B(i) / S(i) if S(i) is greater than
 338:              // *         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
 339:              // *         RCOND*max(S).
 340:              // *
 341:              // *  RANK   (output) INTEGER
 342:              // *         The number of singular values of A greater than RCOND times
 343:              // *         the largest singular value.
 344:              // *
 345:              // *  WORK   (workspace) DOUBLE PRECISION array, dimension at least
 346:              // *         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
 347:              // *         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
 348:              // *
 349:              // *  IWORK  (workspace) INTEGER array, dimension at least
 350:              // *         (3*N*NLVL + 11*N)
 351:              // *
 352:              // *  INFO   (output) INTEGER
 353:              // *         = 0:  successful exit.
 354:              // *         < 0:  if INFO = -i, the i-th argument had an illegal value.
 355:              // *         > 0:  The algorithm failed to compute an singular value while
 356:              // *               working on the submatrix lying in rows and columns
 357:              // *               INFO/(N+1) through MOD(INFO,N+1).
 358:              // *
 359:              // *  Further Details
 360:              // *  ===============
 361:              // *
 362:              // *  Based on contributions by
 363:              // *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
 364:              // *       California at Berkeley, USA
 365:              // *     Osni Marques, LBNL/NERSC, USA
 366:              // *
 367:              // *  =====================================================================
 368:              // *
 369:              // *     .. Parameters ..
 370:              // *     ..
 371:              // *     .. Local Scalars ..
 372:              // *     ..
 373:              // *     .. External Functions ..
 374:              // *     ..
 375:              // *     .. External Subroutines ..
 376:              // *     ..
 377:              // *     .. Intrinsic Functions ..
 378:              //      INTRINSIC          ABS, DBLE, INT, LOG, SIGN;
 379:              // *     ..
 380:              // *     .. Executable Statements ..
 381:              // *
 382:              // *     Test the input parameters.
 383:              // *
 384:   
 385:              #endregion
 386:   
 387:   
 388:              #region Body
 389:              
 390:              INFO = 0;
 391:              // *
 392:              if (N < 0)
 393:              {
 394:                  INFO =  - 3;
 395:              }
 396:              else
 397:              {
 398:                  if (NRHS < 1)
 399:                  {
 400:                      INFO =  - 4;
 401:                  }
 402:                  else
 403:                  {
 404:                      if ((LDB < 1) || (LDB < N))
 405:                      {
 406:                          INFO =  - 8;
 407:                      }
 408:                  }
 409:              }
 410:              if (INFO != 0)
 411:              {
 412:                  this._xerbla.Run("DLALSD",  - INFO);
 413:                  return;
 414:              }
 415:              // *
 416:              EPS = this._dlamch.Run("Epsilon");
 417:              // *
 418:              // *     Set up the tolerance.
 419:              // *
 420:              if ((RCOND <= ZERO) || (RCOND >= ONE))
 421:              {
 422:                  RCND = EPS;
 423:              }
 424:              else
 425:              {
 426:                  RCND = RCOND;
 427:              }
 428:              // *
 429:              RANK = 0;
 430:              // *
 431:              // *     Quick return if possible.
 432:              // *
 433:              if (N == 0)
 434:              {
 435:                  return;
 436:              }
 437:              else
 438:              {
 439:                  if (N == 1)
 440:                  {
 441:                      if (D[1 + o_d] == ZERO)
 442:                      {
 443:                          this._dlaset.Run("A", 1, NRHS, ZERO, ZERO, ref B, offset_b
 444:                                           , LDB);
 445:                      }
 446:                      else
 447:                      {
 448:                          RANK = 1;
 449:                          this._dlascl.Run("G", 0, 0, D[1 + o_d], ONE, 1
 450:                                           , NRHS, ref B, offset_b, LDB, ref INFO);
 451:                          D[1 + o_d] = Math.Abs(D[1 + o_d]);
 452:                      }
 453:                      return;
 454:                  }
 455:              }
 456:              // *
 457:              // *     Rotate the matrix if it is lower bidiagonal.
 458:              // *
 459:              if (UPLO == "L")
 460:              {
 461:                  for (I = 1; I <= N - 1; I++)
 462:                  {
 463:                      this._dlartg.Run(D[I + o_d], E[I + o_e], ref CS, ref SN, ref R);
 464:                      D[I + o_d] = R;
 465:                      E[I + o_e] = SN * D[I + 1 + o_d];
 466:                      D[I + 1 + o_d] = CS * D[I + 1 + o_d];
 467:                      if (NRHS == 1)
 468:                      {
 469:                          this._drot.Run(1, ref B, I+1 * LDB + o_b, 1, ref B, I + 1+1 * LDB + o_b, 1, CS
 470:                                         , SN);
 471:                      }
 472:                      else
 473:                      {
 474:                          WORK[I * 2 - 1 + o_work] = CS;
 475:                          WORK[I * 2 + o_work] = SN;
 476:                      }
 477:                  }
 478:                  if (NRHS > 1)
 479:                  {
 480:                      for (I = 1; I <= NRHS; I++)
 481:                      {
 482:                          for (J = 1; J <= N - 1; J++)
 483:                          {
 484:                              CS = WORK[J * 2 - 1 + o_work];
 485:                              SN = WORK[J * 2 + o_work];
 486:                              this._drot.Run(1, ref B, J+I * LDB + o_b, 1, ref B, J + 1+I * LDB + o_b, 1, CS
 487:                                             , SN);
 488:                          }
 489:                      }
 490:                  }
 491:              }
 492:              // *
 493:              // *     Scale.
 494:              // *
 495:              NM1 = N - 1;
 496:              ORGNRM = this._dlanst.Run("M", N, D, offset_d, E, offset_e);
 497:              if (ORGNRM == ZERO)
 498:              {
 499:                  this._dlaset.Run("A", N, NRHS, ZERO, ZERO, ref B, offset_b
 500:                                   , LDB);
 501:                  return;
 502:              }
 503:              // *
 504:              this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
 505:                               , 1, ref D, offset_d, N, ref INFO);
 506:              this._dlascl.Run("G", 0, 0, ORGNRM, ONE, NM1
 507:                               , 1, ref E, offset_e, NM1, ref INFO);
 508:              // *
 509:              // *     If N is smaller than the minimum divide size SMLSIZ, then solve
 510:              // *     the problem with another solver.
 511:              // *
 512:              if (N <= SMLSIZ)
 513:              {
 514:                  NWORK = 1 + N * N;
 515:                  this._dlaset.Run("A", N, N, ZERO, ONE, ref WORK, offset_work
 516:                                   , N);
 517:                  this._dlasdq.Run("U", 0, N, N, 0, NRHS
 518:                                   , ref D, offset_d, ref E, offset_e, ref WORK, offset_work, N, ref WORK, offset_work, N
 519:                                   , ref B, offset_b, LDB, ref WORK, NWORK + o_work, ref INFO);
 520:                  if (INFO != 0)
 521:                  {
 522:                      return;
 523:                  }
 524:                  TOL = RCND * Math.Abs(D[this._idamax.Run(N, D, offset_d, 1) + o_d]);
 525:                  for (I = 1; I <= N; I++)
 526:                  {
 527:                      if (D[I + o_d] <= TOL)
 528:                      {
 529:                          this._dlaset.Run("A", 1, NRHS, ZERO, ZERO, ref B, I+1 * LDB + o_b
 530:                                           , LDB);
 531:                      }
 532:                      else
 533:                      {
 534:                          this._dlascl.Run("G", 0, 0, D[I + o_d], ONE, 1
 535:                                           , NRHS, ref B, I+1 * LDB + o_b, LDB, ref INFO);
 536:                          RANK = RANK + 1;
 537:                      }
 538:                  }
 539:                  this._dgemm.Run("T", "N", N, NRHS, N, ONE
 540:                                  , WORK, offset_work, N, B, offset_b, LDB, ZERO, ref WORK, NWORK + o_work
 541:                                  , N);
 542:                  this._dlacpy.Run("A", N, NRHS, WORK, NWORK + o_work, N, ref B, offset_b
 543:                                   , LDB);
 544:                  // *
 545:                  // *        Unscale.
 546:                  // *
 547:                  this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
 548:                                   , 1, ref D, offset_d, N, ref INFO);
 549:                  this._dlasrt.Run("D", N, ref D, offset_d, ref INFO);
 550:                  this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
 551:                                   , NRHS, ref B, offset_b, LDB, ref INFO);
 552:                  // *
 553:                  return;
 554:              }
 555:              // *
 556:              // *     Book-keeping and setting up some constants.
 557:              // *
 558:              NLVL = Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(N) / Convert.ToDouble(SMLSIZ + 1)) / Math.Log(TWO))) + 1;
 559:              // *
 560:              SMLSZP = SMLSIZ + 1;
 561:              // *
 562:              U = 1;
 563:              VT = 1 + SMLSIZ * N;
 564:              DIFL = VT + SMLSZP * N;
 565:              DIFR = DIFL + NLVL * N;
 566:              Z = DIFR + NLVL * N * 2;
 567:              C = Z + NLVL * N;
 568:              S = C + N;
 569:              POLES = S + N;
 570:              GIVNUM = POLES + 2 * NLVL * N;
 571:              BX = GIVNUM + 2 * NLVL * N;
 572:              NWORK = BX + N * NRHS;
 573:              // *
 574:              SIZEI = 1 + N;
 575:              K = SIZEI + N;
 576:              GIVPTR = K + N;
 577:              PERM = GIVPTR + N;
 578:              GIVCOL = PERM + NLVL * N;
 579:              IWK = GIVCOL + NLVL * N * 2;
 580:              // *
 581:              ST = 1;
 582:              SQRE = 0;
 583:              ICMPQ1 = 1;
 584:              ICMPQ2 = 0;
 585:              NSUB = 0;
 586:              // *
 587:              for (I = 1; I <= N; I++)
 588:              {
 589:                  if (Math.Abs(D[I + o_d]) < EPS)
 590:                  {
 591:                      D[I + o_d] = FortranLib.Sign(EPS,D[I + o_d]);
 592:                  }
 593:              }
 594:              // *
 595:              for (I = 1; I <= NM1; I++)
 596:              {
 597:                  if ((Math.Abs(E[I + o_e]) < EPS) || (I == NM1))
 598:                  {
 599:                      NSUB = NSUB + 1;
 600:                      IWORK[NSUB + o_iwork] = ST;
 601:                      // *
 602:                      // *           Subproblem found. First determine its size and then
 603:                      // *           apply divide and conquer on it.
 604:                      // *
 605:                      if (I < NM1)
 606:                      {
 607:                          // *
 608:                          // *              A subproblem with E(I) small for I < NM1.
 609:                          // *
 610:                          NSIZE = I - ST + 1;
 611:                          IWORK[SIZEI + NSUB - 1 + o_iwork] = NSIZE;
 612:                      }
 613:                      else
 614:                      {
 615:                          if (Math.Abs(E[I + o_e]) >= EPS)
 616:                          {
 617:                              // *
 618:                              // *              A subproblem with E(NM1) not too small but I = NM1.
 619:                              // *
 620:                              NSIZE = N - ST + 1;
 621:                              IWORK[SIZEI + NSUB - 1 + o_iwork] = NSIZE;
 622:                          }
 623:                          else
 624:                          {
 625:                              // *
 626:                              // *              A subproblem with E(NM1) small. This implies an
 627:                              // *              1-by-1 subproblem at D(N), which is not solved
 628:                              // *              explicitly.
 629:                              // *
 630:                              NSIZE = I - ST + 1;
 631:                              IWORK[SIZEI + NSUB - 1 + o_iwork] = NSIZE;
 632:                              NSUB = NSUB + 1;
 633:                              IWORK[NSUB + o_iwork] = N;
 634:                              IWORK[SIZEI + NSUB - 1 + o_iwork] = 1;
 635:                              this._dcopy.Run(NRHS, B, N+1 * LDB + o_b, LDB, ref WORK, BX + NM1 + o_work, N);
 636:                          }
 637:                      }
 638:                      ST1 = ST - 1;
 639:                      if (NSIZE == 1)
 640:                      {
 641:                          // *
 642:                          // *              This is a 1-by-1 subproblem and is not solved
 643:                          // *              explicitly.
 644:                          // *
 645:                          this._dcopy.Run(NRHS, B, ST+1 * LDB + o_b, LDB, ref WORK, BX + ST1 + o_work, N);
 646:                      }
 647:                      else
 648:                      {
 649:                          if (NSIZE <= SMLSIZ)
 650:                          {
 651:                              // *
 652:                              // *              This is a small subproblem and is solved by DLASDQ.
 653:                              // *
 654:                              this._dlaset.Run("A", NSIZE, NSIZE, ZERO, ONE, ref WORK, VT + ST1 + o_work
 655:                                               , N);
 656:                              this._dlasdq.Run("U", 0, NSIZE, NSIZE, 0, NRHS
 657:                                               , ref D, ST + o_d, ref E, ST + o_e, ref WORK, VT + ST1 + o_work, N, ref WORK, NWORK + o_work, N
 658:                                               , ref B, ST+1 * LDB + o_b, LDB, ref WORK, NWORK + o_work, ref INFO);
 659:                              if (INFO != 0)
 660:                              {
 661:                                  return;
 662:                              }
 663:                              this._dlacpy.Run("A", NSIZE, NRHS, B, ST+1 * LDB + o_b, LDB, ref WORK, BX + ST1 + o_work
 664:                                               , N);
 665:                          }
 666:                          else
 667:                          {
 668:                              // *
 669:                              // *              A large problem. Solve it using divide and conquer.
 670:                              // *
 671:                              this._dlasda.Run(ICMPQ1, SMLSIZ, NSIZE, SQRE, ref D, ST + o_d, ref E, ST + o_e
 672:                                               , ref WORK, U + ST1 + o_work, N, ref WORK, VT + ST1 + o_work, ref IWORK, K + ST1 + o_iwork, ref WORK, DIFL + ST1 + o_work, ref WORK, DIFR + ST1 + o_work
 673:                                               , ref WORK, Z + ST1 + o_work, ref WORK, POLES + ST1 + o_work, ref IWORK, GIVPTR + ST1 + o_iwork, ref IWORK, GIVCOL + ST1 + o_iwork, N, ref IWORK, PERM + ST1 + o_iwork
 674:                                               , ref WORK, GIVNUM + ST1 + o_work, ref WORK, C + ST1 + o_work, ref WORK, S + ST1 + o_work, ref WORK, NWORK + o_work, ref IWORK, IWK + o_iwork, ref INFO);
 675:                              if (INFO != 0)
 676:                              {
 677:                                  return;
 678:                              }
 679:                              BXST = BX + ST1;
 680:                              this._dlalsa.Run(ICMPQ2, SMLSIZ, NSIZE, NRHS, ref B, ST+1 * LDB + o_b, LDB
 681:                                               , ref WORK, BXST + o_work, N, WORK, U + ST1 + o_work, N, WORK, VT + ST1 + o_work, IWORK, K + ST1 + o_iwork
 682:                                               , WORK, DIFL + ST1 + o_work, WORK, DIFR + ST1 + o_work, WORK, Z + ST1 + o_work, WORK, POLES + ST1 + o_work, IWORK, GIVPTR + ST1 + o_iwork, IWORK, GIVCOL + ST1 + o_iwork
 683:                                               , N, IWORK, PERM + ST1 + o_iwork, WORK, GIVNUM + ST1 + o_work, WORK, C + ST1 + o_work, WORK, S + ST1 + o_work, ref WORK, NWORK + o_work
 684:                                               , ref IWORK, IWK + o_iwork, ref INFO);
 685:                              if (INFO != 0)
 686:                              {
 687:                                  return;
 688:                              }
 689:                          }
 690:                      }
 691:                      ST = I + 1;
 692:                  }
 693:              }
 694:              // *
 695:              // *     Apply the singular values and treat the tiny ones as zero.
 696:              // *
 697:              TOL = RCND * Math.Abs(D[this._idamax.Run(N, D, offset_d, 1) + o_d]);
 698:              // *
 699:              for (I = 1; I <= N; I++)
 700:              {
 701:                  // *
 702:                  // *        Some of the elements in D can be negative because 1-by-1
 703:                  // *        subproblems were not solved explicitly.
 704:                  // *
 705:                  if (Math.Abs(D[I + o_d]) <= TOL)
 706:                  {
 707:                      this._dlaset.Run("A", 1, NRHS, ZERO, ZERO, ref WORK, BX + I - 1 + o_work
 708:                                       , N);
 709:                  }
 710:                  else
 711:                  {
 712:                      RANK = RANK + 1;
 713:                      this._dlascl.Run("G", 0, 0, D[I + o_d], ONE, 1
 714:                                       , NRHS, ref WORK, BX + I - 1 + o_work, N, ref INFO);
 715:                  }
 716:                  D[I + o_d] = Math.Abs(D[I + o_d]);
 717:              }
 718:              // *
 719:              // *     Now apply back the right singular vectors.
 720:              // *
 721:              ICMPQ2 = 1;
 722:              for (I = 1; I <= NSUB; I++)
 723:              {
 724:                  ST = IWORK[I + o_iwork];
 725:                  ST1 = ST - 1;
 726:                  NSIZE = IWORK[SIZEI + I - 1 + o_iwork];
 727:                  BXST = BX + ST1;
 728:                  if (NSIZE == 1)
 729:                  {
 730:                      this._dcopy.Run(NRHS, WORK, BXST + o_work, N, ref B, ST+1 * LDB + o_b, LDB);
 731:                  }
 732:                  else
 733:                  {
 734:                      if (NSIZE <= SMLSIZ)
 735:                      {
 736:                          this._dgemm.Run("T", "N", NSIZE, NRHS, NSIZE, ONE
 737:                                          , WORK, VT + ST1 + o_work, N, WORK, BXST + o_work, N, ZERO, ref B, ST+1 * LDB + o_b
 738:                                          , LDB);
 739:                      }
 740:                      else
 741:                      {
 742:                          this._dlalsa.Run(ICMPQ2, SMLSIZ, NSIZE, NRHS, ref WORK, BXST + o_work, N
 743:                                           , ref B, ST+1 * LDB + o_b, LDB, WORK, U + ST1 + o_work, N, WORK, VT + ST1 + o_work, IWORK, K + ST1 + o_iwork
 744:                                           , WORK, DIFL + ST1 + o_work, WORK, DIFR + ST1 + o_work, WORK, Z + ST1 + o_work, WORK, POLES + ST1 + o_work, IWORK, GIVPTR + ST1 + o_iwork, IWORK, GIVCOL + ST1 + o_iwork
 745:                                           , N, IWORK, PERM + ST1 + o_iwork, WORK, GIVNUM + ST1 + o_work, WORK, C + ST1 + o_work, WORK, S + ST1 + o_work, ref WORK, NWORK + o_work
 746:                                           , ref IWORK, IWK + o_iwork, ref INFO);
 747:                          if (INFO != 0)
 748:                          {
 749:                              return;
 750:                          }
 751:                      }
 752:                  }
 753:              }
 754:              // *
 755:              // *     Unscale and sort the singular values.
 756:              // *
 757:              this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
 758:                               , 1, ref D, offset_d, N, ref INFO);
 759:              this._dlasrt.Run("D", N, ref D, offset_d, ref INFO);
 760:              this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
 761:                               , NRHS, ref B, offset_b, LDB, ref INFO);
 762:              // *
 763:              return;
 764:              // *
 765:              // *     End of DLALSD
 766:              // *
 767:   
 768:              #endregion
 769:   
 770:          }
 771:      }
 772:  }