`   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:  }`