`   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 driver routine (version 3.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// November 2006`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DGELSD computes the minimum-norm solution to a real linear least`
`  27:      /// squares problem:`
`  28:      /// minimize 2-norm(| b - A*x |)`
`  29:      /// using the singular value decomposition (SVD) of A. A is an M-by-N`
`  30:      /// matrix which may be rank-deficient.`
`  31:      /// `
`  32:      /// Several right hand side vectors b and solution vectors x can be`
`  33:      /// handled in a single call; they are stored as the columns of the`
`  34:      /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution`
`  35:      /// matrix X.`
`  36:      /// `
`  37:      /// The problem is solved in three steps:`
`  38:      /// (1) Reduce the coefficient matrix A to bidiagonal form with`
`  39:      /// Householder transformations, reducing the original problem`
`  40:      /// into a "bidiagonal least squares problem" (BLS)`
`  41:      /// (2) Solve the BLS using a divide and conquer approach.`
`  42:      /// (3) Apply back all the Householder tranformations to solve`
`  43:      /// the original least squares problem.`
`  44:      /// `
`  45:      /// The effective rank of A is determined by treating as zero those`
`  46:      /// singular values which are less than RCOND times the largest singular`
`  47:      /// value.`
`  48:      /// `
`  49:      /// The divide and conquer algorithm makes very mild assumptions about`
`  50:      /// floating point arithmetic. It will work on machines with a guard`
`  51:      /// digit in add/subtract, or on those binary machines without guard`
`  52:      /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or`
`  53:      /// Cray-2. It could conceivably fail on hexadecimal or decimal machines`
`  54:      /// without guard digits, but we know of none.`
`  55:      /// `
`  56:      ///</summary>`
`  57:      public class DGELSD`
`  58:      {`
`  59:      `
`  60:   `
`  61:          #region Dependencies`
`  62:          `
`  63:          DGEBRD _dgebrd; DGELQF _dgelqf; DGEQRF _dgeqrf; DLABAD _dlabad; DLACPY _dlacpy; DLALSD _dlalsd; DLASCL _dlascl; `
`  64:          DLASET _dlaset;DORMBR _dormbr; DORMLQ _dormlq; DORMQR _dormqr; XERBLA _xerbla; ILAENV _ilaenv; DLAMCH _dlamch; `
`  65:          DLANGE _dlange;`
`  66:   `
`  67:          #endregion`
`  68:   `
`  69:   `
`  70:          #region Fields`
`  71:          `
`  72:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; bool LQUERY = false; int IASCL = 0; `
`  73:          int IBSCL = 0;int IE = 0; int IL = 0; int ITAU = 0; int ITAUP = 0; int ITAUQ = 0; int LDWORK = 0; int MAXMN = 0; `
`  74:          int MAXWRK = 0;int MINMN = 0; int MINWRK = 0; int MM = 0; int MNTHR = 0; int NLVL = 0; int NWORK = 0; int SMLSIZ = 0; `
`  75:          int WLALSD = 0;double ANRM = 0; double BIGNUM = 0; double BNRM = 0; double EPS = 0; double SFMIN = 0; double SMLNUM = 0; `
`  76:   `
`  77:          #endregion`
`  78:   `
`  79:          public DGELSD(DGEBRD dgebrd, DGELQF dgelqf, DGEQRF dgeqrf, DLABAD dlabad, DLACPY dlacpy, DLALSD dlalsd, DLASCL dlascl, DLASET dlaset, DORMBR dormbr, DORMLQ dormlq`
`  80:                        , DORMQR dormqr, XERBLA xerbla, ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange)`
`  81:          {`
`  82:      `
`  83:   `
`  84:              #region Set Dependencies`
`  85:              `
`  86:              this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgeqrf = dgeqrf; this._dlabad = dlabad; this._dlacpy = dlacpy; `
`  87:              this._dlalsd = dlalsd;this._dlascl = dlascl; this._dlaset = dlaset; this._dormbr = dormbr; this._dormlq = dormlq; `
`  88:              this._dormqr = dormqr;this._xerbla = xerbla; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange; `
`  89:   `
`  90:              #endregion`
`  91:   `
`  92:          }`
`  93:      `
`  94:          public DGELSD()`
`  95:          {`
`  96:      `
`  97:   `
`  98:              #region Dependencies (Initialization)`
`  99:              `
` 100:              LSAME lsame = new LSAME();`
` 101:              XERBLA xerbla = new XERBLA();`
` 102:              DLAMC3 dlamc3 = new DLAMC3();`
` 103:              DLAPY2 dlapy2 = new DLAPY2();`
` 104:              DNRM2 dnrm2 = new DNRM2();`
` 105:              DSCAL dscal = new DSCAL();`
` 106:              IEEECK ieeeck = new IEEECK();`
` 107:              IPARMQ iparmq = new IPARMQ();`
` 108:              DCOPY dcopy = new DCOPY();`
` 109:              DLABAD dlabad = new DLABAD();`
` 110:              IDAMAX idamax = new IDAMAX();`
` 111:              DLASSQ dlassq = new DLASSQ();`
` 112:              DROT drot = new DROT();`
` 113:              DLASDT dlasdt = new DLASDT();`
` 114:              DLAMRG dlamrg = new DLAMRG();`
` 115:              DLASD5 dlasd5 = new DLASD5();`
` 116:              DDOT ddot = new DDOT();`
` 117:              DLAS2 dlas2 = new DLAS2();`
` 118:              DLASQ5 dlasq5 = new DLASQ5();`
` 119:              DLAZQ4 dlazq4 = new DLAZQ4();`
` 120:              DSWAP dswap = new DSWAP();`
` 121:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
` 122:              DGER dger = new DGER(xerbla);`
` 123:              DLARF dlarf = new DLARF(dgemv, dger, lsame);`
` 124:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 125:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 126:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 127:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 128:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 129:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
` 130:              DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);`
` 131:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
` 132:              DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);`
` 133:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 134:              DGEBRD dgebrd = new DGEBRD(dgebd2, dgemm, dlabrd, xerbla, ilaenv);`
` 135:              DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);`
` 136:              DTRMM dtrmm = new DTRMM(lsame, xerbla);`
` 137:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);`
` 138:              DTRMV dtrmv = new DTRMV(lsame, xerbla);`
` 139:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);`
` 140:              DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);`
` 141:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);`
` 142:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);`
` 143:              DLACPY dlacpy = new DLACPY(lsame);`
` 144:              DLANST dlanst = new DLANST(lsame, dlassq);`
` 145:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);`
` 146:              DLALS0 dlals0 = new DLALS0(dcopy, dgemv, dlacpy, dlascl, drot, dscal, xerbla, dlamc3, dnrm2);`
` 147:              DLALSA dlalsa = new DLALSA(dcopy, dgemm, dlals0, dlasdt, xerbla);`
` 148:              DLARTG dlartg = new DLARTG(dlamch);`
` 149:              DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);`
` 150:              DLAED6 dlaed6 = new DLAED6(dlamch);`
` 151:              DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);`
` 152:              DLASET dlaset = new DLASET(lsame);`
` 153:              DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);`
` 154:              DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);`
` 155:              DLASQ6 dlasq6 = new DLASQ6(dlamch);`
` 156:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);`
` 157:              DLASRT dlasrt = new DLASRT(lsame, xerbla);`
` 158:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);`
` 159:              DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);`
` 160:              DLASR dlasr = new DLASR(lsame, xerbla);`
` 161:              DLASV2 dlasv2 = new DLASV2(dlamch);`
` 162:              DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap`
` 163:                                         , xerbla);`
` 164:              DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);`
` 165:              DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);`
` 166:              DLALSD dlalsd = new DLALSD(idamax, dlamch, dlanst, dcopy, dgemm, dlacpy, dlalsa, dlartg, dlascl, dlasda`
` 167:                                         , dlasdq, dlaset, dlasrt, drot, xerbla);`
` 168:              DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);`
` 169:              DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);`
` 170:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);`
` 171:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);`
` 172:              DORMBR dormbr = new DORMBR(lsame, ilaenv, dormlq, dormqr, xerbla);`
` 173:              DLANGE dlange = new DLANGE(dlassq, lsame);`
` 174:   `
` 175:              #endregion`
` 176:   `
` 177:   `
` 178:              #region Set Dependencies`
` 179:              `
` 180:              this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgeqrf = dgeqrf; this._dlabad = dlabad; this._dlacpy = dlacpy; `
` 181:              this._dlalsd = dlalsd;this._dlascl = dlascl; this._dlaset = dlaset; this._dormbr = dormbr; this._dormlq = dormlq; `
` 182:              this._dormqr = dormqr;this._xerbla = xerbla; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlange = dlange; `
` 183:   `
` 184:              #endregion`
` 185:   `
` 186:          }`
` 187:          /// <summary>`
` 188:          /// Purpose`
` 189:          /// =======`
` 190:          /// `
` 191:          /// DGELSD computes the minimum-norm solution to a real linear least`
` 192:          /// squares problem:`
` 193:          /// minimize 2-norm(| b - A*x |)`
` 194:          /// using the singular value decomposition (SVD) of A. A is an M-by-N`
` 195:          /// matrix which may be rank-deficient.`
` 196:          /// `
` 197:          /// Several right hand side vectors b and solution vectors x can be`
` 198:          /// handled in a single call; they are stored as the columns of the`
` 199:          /// M-by-NRHS right hand side matrix B and the N-by-NRHS solution`
` 200:          /// matrix X.`
` 201:          /// `
` 202:          /// The problem is solved in three steps:`
` 203:          /// (1) Reduce the coefficient matrix A to bidiagonal form with`
` 204:          /// Householder transformations, reducing the original problem`
` 205:          /// into a "bidiagonal least squares problem" (BLS)`
` 206:          /// (2) Solve the BLS using a divide and conquer approach.`
` 207:          /// (3) Apply back all the Householder tranformations to solve`
` 208:          /// the original least squares problem.`
` 209:          /// `
` 210:          /// The effective rank of A is determined by treating as zero those`
` 211:          /// singular values which are less than RCOND times the largest singular`
` 212:          /// value.`
` 213:          /// `
` 214:          /// The divide and conquer algorithm makes very mild assumptions about`
` 215:          /// floating point arithmetic. It will work on machines with a guard`
` 216:          /// digit in add/subtract, or on those binary machines without guard`
` 217:          /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or`
` 218:          /// Cray-2. It could conceivably fail on hexadecimal or decimal machines`
` 219:          /// without guard digits, but we know of none.`
` 220:          /// `
` 221:          ///</summary>`
` 222:          /// <param name="M">`
` 223:          /// (input) INTEGER`
` 224:          /// The number of rows of A. M .GE. 0.`
` 225:          ///</param>`
` 226:          /// <param name="N">`
` 227:          /// (input) INTEGER`
` 228:          /// The number of columns of A. N .GE. 0.`
` 229:          ///</param>`
` 230:          /// <param name="NRHS">`
` 231:          /// (input) INTEGER`
` 232:          /// The number of right hand sides, i.e., the number of columns`
` 233:          /// of the matrices B and X. NRHS .GE. 0.`
` 234:          ///</param>`
` 235:          /// <param name="A">`
` 236:          /// (input) DOUBLE PRECISION array, dimension (LDA,N)`
` 237:          /// On entry, the M-by-N matrix A.`
` 238:          /// On exit, A has been destroyed.`
` 239:          ///</param>`
` 240:          /// <param name="LDA">`
` 241:          /// (input) INTEGER`
` 242:          /// The leading dimension of the array A.  LDA .GE. max(1,M).`
` 243:          ///</param>`
` 244:          /// <param name="B">`
` 245:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)`
` 246:          /// On entry, the M-by-NRHS right hand side matrix B.`
` 247:          /// On exit, B is overwritten by the N-by-NRHS solution`
` 248:          /// matrix X.  If m .GE. n and RANK = n, the residual`
` 249:          /// sum-of-squares for the solution in the i-th column is given`
` 250:          /// by the sum of squares of elements n+1:m in that column.`
` 251:          ///</param>`
` 252:          /// <param name="LDB">`
` 253:          /// (input) INTEGER`
` 254:          /// The leading dimension of the array B. LDB .GE. max(1,max(M,N)).`
` 255:          ///</param>`
` 256:          /// <param name="S">`
` 257:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))`
` 258:          /// The singular values of A in decreasing order.`
` 259:          /// The condition number of A in the 2-norm = S(1)/S(min(m,n)).`
` 260:          ///</param>`
` 261:          /// <param name="RCOND">`
` 262:          /// (input) DOUBLE PRECISION`
` 263:          /// RCOND is used to determine the effective rank of A.`
` 264:          /// Singular values S(i) .LE. RCOND*S(1) are treated as zero.`
` 265:          /// If RCOND .LT. 0, machine precision is used instead.`
` 266:          ///</param>`
` 267:          /// <param name="RANK">`
` 268:          /// (output) INTEGER`
` 269:          /// The effective rank of A, i.e., the number of singular values`
` 270:          /// which are greater than RCOND*S(1).`
` 271:          ///</param>`
` 272:          /// <param name="WORK">`
` 273:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 274:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.`
` 275:          ///</param>`
` 276:          /// <param name="LWORK">`
` 277:          /// (input) INTEGER`
` 278:          /// The dimension of the array WORK. LWORK must be at least 1.`
` 279:          /// The exact minimum amount of workspace needed depends on M,`
` 280:          /// N and NRHS. As long as LWORK is at least`
` 281:          /// 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,`
` 282:          /// if M is greater than or equal to N or`
` 283:          /// 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,`
` 284:          /// if M is less than N, the code will execute correctly.`
` 285:          /// SMLSIZ is returned by ILAENV and is equal to the maximum`
` 286:          /// size of the subproblems at the bottom of the computation`
` 287:          /// tree (usually about 25), and`
` 288:          /// NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )`
` 289:          /// For good performance, LWORK should generally be larger.`
` 290:          /// `
` 291:          /// If LWORK = -1, then a workspace query is assumed; the routine`
` 292:          /// only calculates the optimal size of the WORK array, returns`
` 293:          /// this value as the first entry of the WORK array, and no error`
` 294:          /// message related to LWORK is issued by XERBLA.`
` 295:          ///</param>`
` 296:          /// <param name="IWORK">`
` 297:          /// (workspace) INTEGER array, dimension (MAX(1,LIWORK))`
` 298:          /// LIWORK .GE. 3 * MINMN * NLVL + 11 * MINMN,`
` 299:          /// where MINMN = MIN( M,N ).`
` 300:          ///</param>`
` 301:          /// <param name="INFO">`
` 302:          /// (output) INTEGER`
` 303:          /// = 0:  successful exit`
` 304:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 305:          /// .GT. 0:  the algorithm for computing the SVD failed to converge;`
` 306:          /// if INFO = i, i off-diagonal elements of an intermediate`
` 307:          /// bidiagonal form did not converge to zero.`
` 308:          ///</param>`
` 309:          public void Run(int M, int N, int NRHS, ref double[] A, int offset_a, int LDA, ref double[] B, int offset_b`
` 310:                           , int LDB, ref double[] S, int offset_s, double RCOND, ref int RANK, ref double[] WORK, int offset_work, int LWORK`
` 311:                           , ref int[] IWORK, int offset_iwork, ref int INFO)`
` 312:          {`
` 313:   `
` 314:              #region Array Index Correction`
` 315:              `
` 316:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b;  int o_s = -1 + offset_s; `
` 317:               int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork; `
` 318:   `
` 319:              #endregion`
` 320:   `
` 321:   `
` 322:              #region Prolog`
` 323:              `
` 324:              // *`
` 325:              // *  -- LAPACK driver routine (version 3.1) --`
` 326:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 327:              // *     November 2006`
` 328:              // *`
` 329:              // *     .. Scalar Arguments ..`
` 330:              // *     ..`
` 331:              // *     .. Array Arguments ..`
` 332:              // *     ..`
` 333:              // *`
` 334:              // *  Purpose`
` 335:              // *  =======`
` 336:              // *`
` 337:              // *  DGELSD computes the minimum-norm solution to a real linear least`
` 338:              // *  squares problem:`
` 339:              // *      minimize 2-norm(| b - A*x |)`
` 340:              // *  using the singular value decomposition (SVD) of A. A is an M-by-N`
` 341:              // *  matrix which may be rank-deficient.`
` 342:              // *`
` 343:              // *  Several right hand side vectors b and solution vectors x can be`
` 344:              // *  handled in a single call; they are stored as the columns of the`
` 345:              // *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution`
` 346:              // *  matrix X.`
` 347:              // *`
` 348:              // *  The problem is solved in three steps:`
` 349:              // *  (1) Reduce the coefficient matrix A to bidiagonal form with`
` 350:              // *      Householder transformations, reducing the original problem`
` 351:              // *      into a "bidiagonal least squares problem" (BLS)`
` 352:              // *  (2) Solve the BLS using a divide and conquer approach.`
` 353:              // *  (3) Apply back all the Householder tranformations to solve`
` 354:              // *      the original least squares problem.`
` 355:              // *`
` 356:              // *  The effective rank of A is determined by treating as zero those`
` 357:              // *  singular values which are less than RCOND times the largest singular`
` 358:              // *  value.`
` 359:              // *`
` 360:              // *  The divide and conquer algorithm makes very mild assumptions about`
` 361:              // *  floating point arithmetic. It will work on machines with a guard`
` 362:              // *  digit in add/subtract, or on those binary machines without guard`
` 363:              // *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or`
` 364:              // *  Cray-2. It could conceivably fail on hexadecimal or decimal machines`
` 365:              // *  without guard digits, but we know of none.`
` 366:              // *`
` 367:              // *  Arguments`
` 368:              // *  =========`
` 369:              // *`
` 370:              // *  M       (input) INTEGER`
` 371:              // *          The number of rows of A. M >= 0.`
` 372:              // *`
` 373:              // *  N       (input) INTEGER`
` 374:              // *          The number of columns of A. N >= 0.`
` 375:              // *`
` 376:              // *  NRHS    (input) INTEGER`
` 377:              // *          The number of right hand sides, i.e., the number of columns`
` 378:              // *          of the matrices B and X. NRHS >= 0.`
` 379:              // *`
` 380:              // *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)`
` 381:              // *          On entry, the M-by-N matrix A.`
` 382:              // *          On exit, A has been destroyed.`
` 383:              // *`
` 384:              // *  LDA     (input) INTEGER`
` 385:              // *          The leading dimension of the array A.  LDA >= max(1,M).`
` 386:              // *`
` 387:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)`
` 388:              // *          On entry, the M-by-NRHS right hand side matrix B.`
` 389:              // *          On exit, B is overwritten by the N-by-NRHS solution`
` 390:              // *          matrix X.  If m >= n and RANK = n, the residual`
` 391:              // *          sum-of-squares for the solution in the i-th column is given`
` 392:              // *          by the sum of squares of elements n+1:m in that column.`
` 393:              // *`
` 394:              // *  LDB     (input) INTEGER`
` 395:              // *          The leading dimension of the array B. LDB >= max(1,max(M,N)).`
` 396:              // *`
` 397:              // *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))`
` 398:              // *          The singular values of A in decreasing order.`
` 399:              // *          The condition number of A in the 2-norm = S(1)/S(min(m,n)).`
` 400:              // *`
` 401:              // *  RCOND   (input) DOUBLE PRECISION`
` 402:              // *          RCOND is used to determine the effective rank of A.`
` 403:              // *          Singular values S(i) <= RCOND*S(1) are treated as zero.`
` 404:              // *          If RCOND < 0, machine precision is used instead.`
` 405:              // *`
` 406:              // *  RANK    (output) INTEGER`
` 407:              // *          The effective rank of A, i.e., the number of singular values`
` 408:              // *          which are greater than RCOND*S(1).`
` 409:              // *`
` 410:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 411:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.`
` 412:              // *`
` 413:              // *  LWORK   (input) INTEGER`
` 414:              // *          The dimension of the array WORK. LWORK must be at least 1.`
` 415:              // *          The exact minimum amount of workspace needed depends on M,`
` 416:              // *          N and NRHS. As long as LWORK is at least`
` 417:              // *              12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,`
` 418:              // *          if M is greater than or equal to N or`
` 419:              // *              12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,`
` 420:              // *          if M is less than N, the code will execute correctly.`
` 421:              // *          SMLSIZ is returned by ILAENV and is equal to the maximum`
` 422:              // *          size of the subproblems at the bottom of the computation`
` 423:              // *          tree (usually about 25), and`
` 424:              // *             NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )`
` 425:              // *          For good performance, LWORK should generally be larger.`
` 426:              // *`
` 427:              // *          If LWORK = -1, then a workspace query is assumed; the routine`
` 428:              // *          only calculates the optimal size of the WORK array, returns`
` 429:              // *          this value as the first entry of the WORK array, and no error`
` 430:              // *          message related to LWORK is issued by XERBLA.`
` 431:              // *`
` 432:              // *  IWORK   (workspace) INTEGER array, dimension (MAX(1,LIWORK))`
` 433:              // *          LIWORK >= 3 * MINMN * NLVL + 11 * MINMN,`
` 434:              // *          where MINMN = MIN( M,N ).`
` 435:              // *`
` 436:              // *  INFO    (output) INTEGER`
` 437:              // *          = 0:  successful exit`
` 438:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 439:              // *          > 0:  the algorithm for computing the SVD failed to converge;`
` 440:              // *                if INFO = i, i off-diagonal elements of an intermediate`
` 441:              // *                bidiagonal form did not converge to zero.`
` 442:              // *`
` 443:              // *  Further Details`
` 444:              // *  ===============`
` 445:              // *`
` 446:              // *  Based on contributions by`
` 447:              // *     Ming Gu and Ren-Cang Li, Computer Science Division, University of`
` 448:              // *       California at Berkeley, USA`
` 449:              // *     Osni Marques, LBNL/NERSC, USA`
` 450:              // *`
` 451:              // *  =====================================================================`
` 452:              // *`
` 453:              // *     .. Parameters ..`
` 454:              // *     ..`
` 455:              // *     .. Local Scalars ..`
` 456:              // *     ..`
` 457:              // *     .. External Subroutines ..`
` 458:              // *     ..`
` 459:              // *     .. External Functions ..`
` 460:              // *     ..`
` 461:              // *     .. Intrinsic Functions ..`
` 462:              //      INTRINSIC          DBLE, INT, LOG, MAX, MIN;`
` 463:              // *     ..`
` 464:              // *     .. Executable Statements ..`
` 465:              // *`
` 466:              // *     Test the input arguments.`
` 467:              // *`
` 468:   `
` 469:              #endregion`
` 470:   `
` 471:   `
` 472:              #region Body`
` 473:              `
` 474:              INFO = 0;`
` 475:              MINMN = Math.Min(M, N);`
` 476:              MAXMN = Math.Max(M, N);`
` 477:              MNTHR = this._ilaenv.Run(6, "DGELSD", " ", M, N, NRHS,  - 1);`
` 478:              LQUERY = (LWORK ==  - 1);`
` 479:              if (M < 0)`
` 480:              {`
` 481:                  INFO =  - 1;`
` 482:              }`
` 483:              else`
` 484:              {`
` 485:                  if (N < 0)`
` 486:                  {`
` 487:                      INFO =  - 2;`
` 488:                  }`
` 489:                  else`
` 490:                  {`
` 491:                      if (NRHS < 0)`
` 492:                      {`
` 493:                          INFO =  - 3;`
` 494:                      }`
` 495:                      else`
` 496:                      {`
` 497:                          if (LDA < Math.Max(1, M))`
` 498:                          {`
` 499:                              INFO =  - 5;`
` 500:                          }`
` 501:                          else`
` 502:                          {`
` 503:                              if (LDB < Math.Max(1, MAXMN))`
` 504:                              {`
` 505:                                  INFO =  - 7;`
` 506:                              }`
` 507:                          }`
` 508:                      }`
` 509:                  }`
` 510:              }`
` 511:              // *`
` 512:              SMLSIZ = this._ilaenv.Run(9, "DGELSD", " ", 0, 0, 0, 0);`
` 513:              // *`
` 514:              // *     Compute workspace.`
` 515:              // *     (Note: Comments in the code beginning "Workspace:" describe the`
` 516:              // *     minimal amount of workspace needed at that point in the code,`
` 517:              // *     as well as the preferred amount for good performance.`
` 518:              // *     NB refers to the optimal block size for the immediately`
` 519:              // *     following subroutine, as returned by ILAENV.)`
` 520:              // *`
` 521:              MINWRK = 1;`
` 522:              MINMN = Math.Max(1, MINMN);`
` 523:              NLVL = Math.Max(Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(MINMN) / Convert.ToDouble(SMLSIZ + 1)) / Math.Log(TWO))) + 1, 0);`
` 524:              // *`
` 525:              if (INFO == 0)`
` 526:              {`
` 527:                  MAXWRK = 0;`
` 528:                  MM = M;`
` 529:                  if (M >= N && M >= MNTHR)`
` 530:                  {`
` 531:                      // *`
` 532:                      // *           Path 1a - overdetermined, with many more rows than columns.`
` 533:                      // *`
` 534:                      MM = N;`
` 535:                      MAXWRK = Math.Max(MAXWRK, N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1));`
` 536:                      MAXWRK = Math.Max(MAXWRK, N + NRHS * this._ilaenv.Run(1, "DORMQR", "LT", M, NRHS, N,  - 1));`
` 537:                  }`
` 538:                  if (M >= N)`
` 539:                  {`
` 540:                      // *`
` 541:                      // *           Path 1 - overdetermined or exactly determined.`
` 542:                      // *`
` 543:                      MAXWRK = Math.Max(MAXWRK, 3 * N + (MM + N) * this._ilaenv.Run(1, "DGEBRD", " ", MM, N,  - 1,  - 1));`
` 544:                      MAXWRK = Math.Max(MAXWRK, 3 * N + NRHS * this._ilaenv.Run(1, "DORMBR", "QLT", MM, NRHS, N,  - 1));`
` 545:                      MAXWRK = Math.Max(MAXWRK, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORMBR", "PLN", N, NRHS, N,  - 1));`
` 546:                      WLALSD = 9 * N + 2 * N * SMLSIZ + 8 * N * NLVL + N * NRHS + (int)Math.Pow(SMLSIZ + 1, 2);`
` 547:                      MAXWRK = Math.Max(MAXWRK, 3 * N + WLALSD);`
` 548:                      MINWRK = Math.Max(3 * N + MM, Math.Max(3 * N + NRHS, 3 * N + WLALSD));`
` 549:                  }`
` 550:                  if (N > M)`
` 551:                  {`
` 552:                      WLALSD = 9 * M + 2 * M * SMLSIZ + 8 * M * NLVL + M * NRHS + (int)Math.Pow(SMLSIZ + 1, 2);`
` 553:                      if (N >= MNTHR)`
` 554:                      {`
` 555:                          // *`
` 556:                          // *              Path 2a - underdetermined, with many more columns`
` 557:                          // *              than rows.`
` 558:                          // *`
` 559:                          MAXWRK = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);`
` 560:                          MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));`
` 561:                          MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + NRHS * this._ilaenv.Run(1, "DORMBR", "QLT", M, NRHS, M,  - 1));`
` 562:                          MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + (M - 1) * this._ilaenv.Run(1, "DORMBR", "PLN", M, NRHS, M,  - 1));`
` 563:                          if (NRHS > 1)`
` 564:                          {`
` 565:                              MAXWRK = Math.Max(MAXWRK, M * M + M + M * NRHS);`
` 566:                          }`
` 567:                          else`
` 568:                          {`
` 569:                              MAXWRK = Math.Max(MAXWRK, M * M + 2 * M);`
` 570:                          }`
` 571:                          MAXWRK = Math.Max(MAXWRK, M + NRHS * this._ilaenv.Run(1, "DORMLQ", "LT", N, NRHS, M,  - 1));`
` 572:                          MAXWRK = Math.Max(MAXWRK, M * M + 4 * M + WLALSD);`
` 573:                      }`
` 574:                      else`
` 575:                      {`
` 576:                          // *`
` 577:                          // *              Path 2 - remaining underdetermined cases.`
` 578:                          // *`
` 579:                          MAXWRK = 3 * M + (N + M) * this._ilaenv.Run(1, "DGEBRD", " ", M, N,  - 1,  - 1);`
` 580:                          MAXWRK = Math.Max(MAXWRK, 3 * M + NRHS * this._ilaenv.Run(1, "DORMBR", "QLT", M, NRHS, N,  - 1));`
` 581:                          MAXWRK = Math.Max(MAXWRK, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PLN", N, NRHS, M,  - 1));`
` 582:                          MAXWRK = Math.Max(MAXWRK, 3 * M + WLALSD);`
` 583:                      }`
` 584:                      MINWRK = Math.Max(3 * M + NRHS, Math.Max(3 * M + M, 3 * M + WLALSD));`
` 585:                  }`
` 586:                  MINWRK = Math.Min(MINWRK, MAXWRK);`
` 587:                  WORK[1 + o_work] = MAXWRK;`
` 588:                  if (LWORK < MINWRK && !LQUERY)`
` 589:                  {`
` 590:                      INFO =  - 12;`
` 591:                  }`
` 592:              }`
` 593:              // *`
` 594:              if (INFO != 0)`
` 595:              {`
` 596:                  this._xerbla.Run("DGELSD",  - INFO);`
` 597:                  return;`
` 598:              }`
` 599:              else`
` 600:              {`
` 601:                  if (LQUERY)`
` 602:                  {`
` 603:                      goto LABEL10;`
` 604:                  }`
` 605:              }`
` 606:              // *`
` 607:              // *     Quick return if possible.`
` 608:              // *`
` 609:              if (M == 0 || N == 0)`
` 610:              {`
` 611:                  RANK = 0;`
` 612:                  return;`
` 613:              }`
` 614:              // *`
` 615:              // *     Get machine parameters.`
` 616:              // *`
` 617:              EPS = this._dlamch.Run("P");`
` 618:              SFMIN = this._dlamch.Run("S");`
` 619:              SMLNUM = SFMIN / EPS;`
` 620:              BIGNUM = ONE / SMLNUM;`
` 621:              this._dlabad.Run(ref SMLNUM, ref BIGNUM);`
` 622:              // *`
` 623:              // *     Scale A if max entry outside range [SMLNUM,BIGNUM].`
` 624:              // *`
` 625:              ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref WORK, offset_work);`
` 626:              IASCL = 0;`
` 627:              if (ANRM > ZERO && ANRM < SMLNUM)`
` 628:              {`
` 629:                  // *`
` 630:                  // *        Scale matrix norm up to SMLNUM.`
` 631:                  // *`
` 632:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M`
` 633:                                   , N, ref A, offset_a, LDA, ref INFO);`
` 634:                  IASCL = 1;`
` 635:              }`
` 636:              else`
` 637:              {`
` 638:                  if (ANRM > BIGNUM)`
` 639:                  {`
` 640:                      // *`
` 641:                      // *        Scale matrix norm down to BIGNUM.`
` 642:                      // *`
` 643:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M`
` 644:                                       , N, ref A, offset_a, LDA, ref INFO);`
` 645:                      IASCL = 2;`
` 646:                  }`
` 647:                  else`
` 648:                  {`
` 649:                      if (ANRM == ZERO)`
` 650:                      {`
` 651:                          // *`
` 652:                          // *        Matrix all zero. Return zero solution.`
` 653:                          // *`
` 654:                          this._dlaset.Run("F", Math.Max(M, N), NRHS, ZERO, ZERO, ref B, offset_b`
` 655:                                           , LDB);`
` 656:                          this._dlaset.Run("F", MINMN, 1, ZERO, ZERO, ref S, offset_s`
` 657:                                           , 1);`
` 658:                          RANK = 0;`
` 659:                          goto LABEL10;`
` 660:                      }`
` 661:                  }`
` 662:              }`
` 663:              // *`
` 664:              // *     Scale B if max entry outside range [SMLNUM,BIGNUM].`
` 665:              // *`
` 666:              BNRM = this._dlange.Run("M", M, NRHS, B, offset_b, LDB, ref WORK, offset_work);`
` 667:              IBSCL = 0;`
` 668:              if (BNRM > ZERO && BNRM < SMLNUM)`
` 669:              {`
` 670:                  // *`
` 671:                  // *        Scale matrix norm up to SMLNUM.`
` 672:                  // *`
` 673:                  this._dlascl.Run("G", 0, 0, BNRM, SMLNUM, M`
` 674:                                   , NRHS, ref B, offset_b, LDB, ref INFO);`
` 675:                  IBSCL = 1;`
` 676:              }`
` 677:              else`
` 678:              {`
` 679:                  if (BNRM > BIGNUM)`
` 680:                  {`
` 681:                      // *`
` 682:                      // *        Scale matrix norm down to BIGNUM.`
` 683:                      // *`
` 684:                      this._dlascl.Run("G", 0, 0, BNRM, BIGNUM, M`
` 685:                                       , NRHS, ref B, offset_b, LDB, ref INFO);`
` 686:                      IBSCL = 2;`
` 687:                  }`
` 688:              }`
` 689:              // *`
` 690:              // *     If M < N make sure certain entries of B are zero.`
` 691:              // *`
` 692:              if (M < N)`
` 693:              {`
` 694:                  this._dlaset.Run("F", N - M, NRHS, ZERO, ZERO, ref B, M + 1+1 * LDB + o_b`
` 695:                                   , LDB);`
` 696:              }`
` 697:              // *`
` 698:              // *     Overdetermined case.`
` 699:              // *`
` 700:              if (M >= N)`
` 701:              {`
` 702:                  // *`
` 703:                  // *        Path 1 - overdetermined or exactly determined.`
` 704:                  // *`
` 705:                  MM = M;`
` 706:                  if (M >= MNTHR)`
` 707:                  {`
` 708:                      // *`
` 709:                      // *           Path 1a - overdetermined, with many more rows than columns.`
` 710:                      // *`
` 711:                      MM = N;`
` 712:                      ITAU = 1;`
` 713:                      NWORK = ITAU + N;`
` 714:                      // *`
` 715:                      // *           Compute A=Q*R.`
` 716:                      // *           (Workspace: need 2*N, prefer N+N*NB)`
` 717:                      // *`
` 718:                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
` 719:                                       , LWORK - NWORK + 1, ref INFO);`
` 720:                      // *`
` 721:                      // *           Multiply B by transpose(Q).`
` 722:                      // *           (Workspace: need N+NRHS, prefer N+NRHS*NB)`
` 723:                      // *`
` 724:                      this._dormqr.Run("L", "T", M, NRHS, N, ref A, offset_a`
` 725:                                       , LDA, WORK, ITAU + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work, LWORK - NWORK + 1`
` 726:                                       , ref INFO);`
` 727:                      // *`
` 728:                      // *           Zero out below R.`
` 729:                      // *`
` 730:                      if (N > 1)`
` 731:                      {`
` 732:                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a`
` 733:                                           , LDA);`
` 734:                      }`
` 735:                  }`
` 736:                  // *`
` 737:                  IE = 1;`
` 738:                  ITAUQ = IE + N;`
` 739:                  ITAUP = ITAUQ + N;`
` 740:                  NWORK = ITAUP + N;`
` 741:                  // *`
` 742:                  // *        Bidiagonalize R in A.`
` 743:                  // *        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)`
` 744:                  // *`
` 745:                  this._dgebrd.Run(MM, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
` 746:                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref INFO);`
` 747:                  // *`
` 748:                  // *        Multiply B by transpose of left bidiagonalizing vectors of R.`
` 749:                  // *        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)`
` 750:                  // *`
` 751:                  this._dormbr.Run("Q", "L", "T", MM, NRHS, N`
` 752:                                   , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work`
` 753:                                   , LWORK - NWORK + 1, ref INFO);`
` 754:                  // *`
` 755:                  // *        Solve the bidiagonal least squares problem.`
` 756:                  // *`
` 757:                  this._dlalsd.Run("U", SMLSIZ, N, NRHS, ref S, offset_s, ref WORK, IE + o_work`
` 758:                                   , ref B, offset_b, LDB, RCOND, ref RANK, ref WORK, NWORK + o_work, ref IWORK, offset_iwork`
` 759:                                   , ref INFO);`
` 760:                  if (INFO != 0)`
` 761:                  {`
` 762:                      goto LABEL10;`
` 763:                  }`
` 764:                  // *`
` 765:                  // *        Multiply B by right bidiagonalizing vectors of R.`
` 766:                  // *`
` 767:                  this._dormbr.Run("P", "L", "N", N, NRHS, N`
` 768:                                   , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work`
` 769:                                   , LWORK - NWORK + 1, ref INFO);`
` 770:                  // *`
` 771:              }`
` 772:              else`
` 773:              {`
` 774:                  if (N >= MNTHR && LWORK >= 4 * M + M * M + Math.Max(M, Math.Max(2 * M - 4, Math.Max(NRHS, Math.Max(N - 3 * M, WLALSD)))))`
` 775:                  {`
` 776:                      // *`
` 777:                      // *        Path 2a - underdetermined, with many more columns than rows`
` 778:                      // *        and sufficient workspace for an efficient algorithm.`
` 779:                      // *`
` 780:                      LDWORK = M;`
` 781:                      if (LWORK >= Math.Max(4 * M + M * LDA + Math.Max(M, Math.Max(2 * M - 4, Math.Max(NRHS, N - 3 * M))), Math.Max(M * LDA + M + M * NRHS, 4 * M + M * LDA + WLALSD))) LDWORK = LDA;`
` 782:                      ITAU = 1;`
` 783:                      NWORK = M + 1;`
` 784:                      // *`
` 785:                      // *        Compute A=L*Q.`
` 786:                      // *        (Workspace: need 2*M, prefer M+M*NB)`
` 787:                      // *`
` 788:                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
` 789:                                       , LWORK - NWORK + 1, ref INFO);`
` 790:                      IL = NWORK;`
` 791:                      // *`
` 792:                      // *        Copy L to WORK(IL), zeroing out above its diagonal.`
` 793:                      // *`
` 794:                      this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IL + o_work`
` 795:                                       , LDWORK);`
` 796:                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IL + LDWORK + o_work`
` 797:                                       , LDWORK);`
` 798:                      IE = IL + LDWORK * M;`
` 799:                      ITAUQ = IE + M;`
` 800:                      ITAUP = ITAUQ + M;`
` 801:                      NWORK = ITAUP + M;`
` 802:                      // *`
` 803:                      // *        Bidiagonalize L in WORK(IL).`
` 804:                      // *        (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)`
` 805:                      // *`
` 806:                      this._dgebrd.Run(M, M, ref WORK, IL + o_work, LDWORK, ref S, offset_s, ref WORK, IE + o_work`
` 807:                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref INFO);`
` 808:                      // *`
` 809:                      // *        Multiply B by transpose of left bidiagonalizing vectors of L.`
` 810:                      // *        (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)`
` 811:                      // *`
` 812:                      this._dormbr.Run("Q", "L", "T", M, NRHS, M`
` 813:                                       , ref WORK, IL + o_work, LDWORK, WORK, ITAUQ + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work`
` 814:                                       , LWORK - NWORK + 1, ref INFO);`
` 815:                      // *`
` 816:                      // *        Solve the bidiagonal least squares problem.`
` 817:                      // *`
` 818:                      this._dlalsd.Run("U", SMLSIZ, M, NRHS, ref S, offset_s, ref WORK, IE + o_work`
` 819:                                       , ref B, offset_b, LDB, RCOND, ref RANK, ref WORK, NWORK + o_work, ref IWORK, offset_iwork`
` 820:                                       , ref INFO);`
` 821:                      if (INFO != 0)`
` 822:                      {`
` 823:                          goto LABEL10;`
` 824:                      }`
` 825:                      // *`
` 826:                      // *        Multiply B by right bidiagonalizing vectors of L.`
` 827:                      // *`
` 828:                      this._dormbr.Run("P", "L", "N", M, NRHS, M`
` 829:                                       , ref WORK, IL + o_work, LDWORK, WORK, ITAUP + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work`
` 830:                                       , LWORK - NWORK + 1, ref INFO);`
` 831:                      // *`
` 832:                      // *        Zero out below first M rows of B.`
` 833:                      // *`
` 834:                      this._dlaset.Run("F", N - M, NRHS, ZERO, ZERO, ref B, M + 1+1 * LDB + o_b`
` 835:                                       , LDB);`
` 836:                      NWORK = ITAU + M;`
` 837:                      // *`
` 838:                      // *        Multiply transpose(Q) by B.`
` 839:                      // *        (Workspace: need M+NRHS, prefer M+NRHS*NB)`
` 840:                      // *`
` 841:                      this._dormlq.Run("L", "T", N, NRHS, M, ref A, offset_a`
` 842:                                       , LDA, WORK, ITAU + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work, LWORK - NWORK + 1`
` 843:                                       , ref INFO);`
` 844:                      // *`
` 845:                  }`
` 846:                  else`
` 847:                  {`
` 848:                      // *`
` 849:                      // *        Path 2 - remaining underdetermined cases.`
` 850:                      // *`
` 851:                      IE = 1;`
` 852:                      ITAUQ = IE + M;`
` 853:                      ITAUP = ITAUQ + M;`
` 854:                      NWORK = ITAUP + M;`
` 855:                      // *`
` 856:                      // *        Bidiagonalize A.`
` 857:                      // *        (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)`
` 858:                      // *`
` 859:                      this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
` 860:                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref INFO);`
` 861:                      // *`
` 862:                      // *        Multiply B by transpose of left bidiagonalizing vectors.`
` 863:                      // *        (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)`
` 864:                      // *`
` 865:                      this._dormbr.Run("Q", "L", "T", M, NRHS, N`
` 866:                                       , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work`
` 867:                                       , LWORK - NWORK + 1, ref INFO);`
` 868:                      // *`
` 869:                      // *        Solve the bidiagonal least squares problem.`
` 870:                      // *`
` 871:                      this._dlalsd.Run("L", SMLSIZ, M, NRHS, ref S, offset_s, ref WORK, IE + o_work`
` 872:                                       , ref B, offset_b, LDB, RCOND, ref RANK, ref WORK, NWORK + o_work, ref IWORK, offset_iwork`
` 873:                                       , ref INFO);`
` 874:                      if (INFO != 0)`
` 875:                      {`
` 876:                          goto LABEL10;`
` 877:                      }`
` 878:                      // *`
` 879:                      // *        Multiply B by right bidiagonalizing vectors of A.`
` 880:                      // *`
` 881:                      this._dormbr.Run("P", "L", "N", N, NRHS, M`
` 882:                                       , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref B, offset_b, LDB, ref WORK, NWORK + o_work`
` 883:                                       , LWORK - NWORK + 1, ref INFO);`
` 884:                      // *`
` 885:                  }`
` 886:              }`
` 887:              // *`
` 888:              // *     Undo scaling.`
` 889:              // *`
` 890:              if (IASCL == 1)`
` 891:              {`
` 892:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, N`
` 893:                                   , NRHS, ref B, offset_b, LDB, ref INFO);`
` 894:                  this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN`
` 895:                                   , 1, ref S, offset_s, MINMN, ref INFO);`
` 896:              }`
` 897:              else`
` 898:              {`
` 899:                  if (IASCL == 2)`
` 900:                  {`
` 901:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, N`
` 902:                                       , NRHS, ref B, offset_b, LDB, ref INFO);`
` 903:                      this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN`
` 904:                                       , 1, ref S, offset_s, MINMN, ref INFO);`
` 905:                  }`
` 906:              }`
` 907:              if (IBSCL == 1)`
` 908:              {`
` 909:                  this._dlascl.Run("G", 0, 0, SMLNUM, BNRM, N`
` 910:                                   , NRHS, ref B, offset_b, LDB, ref INFO);`
` 911:              }`
` 912:              else`
` 913:              {`
` 914:                  if (IBSCL == 2)`
` 915:                  {`
` 916:                      this._dlascl.Run("G", 0, 0, BIGNUM, BNRM, N`
` 917:                                       , NRHS, ref B, offset_b, LDB, ref INFO);`
` 918:                  }`
` 919:              }`
` 920:              // *`
` 921:          LABEL10:;`
` 922:              WORK[1 + o_work] = MAXWRK;`
` 923:              return;`
` 924:              // *`
` 925:              // *     End of DGELSD`
` 926:              // *`
` 927:   `
` 928:              #endregion`
` 929:   `
` 930:          }`
` 931:      }`
` 932:  }`