`   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:      /// DGESDD computes the singular value decomposition (SVD) of a real`
`  27:      /// M-by-N matrix A, optionally computing the left and right singular`
`  28:      /// vectors.  If singular vectors are desired, it uses a`
`  29:      /// divide-and-conquer algorithm.`
`  30:      /// `
`  31:      /// The SVD is written`
`  32:      /// `
`  33:      /// A = U * SIGMA * transpose(V)`
`  34:      /// `
`  35:      /// where SIGMA is an M-by-N matrix which is zero except for its`
`  36:      /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and`
`  37:      /// V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA`
`  38:      /// are the singular values of A; they are real and non-negative, and`
`  39:      /// are returned in descending order.  The first min(m,n) columns of`
`  40:      /// U and V are the left and right singular vectors of A.`
`  41:      /// `
`  42:      /// Note that the routine returns VT = V**T, not V.`
`  43:      /// `
`  44:      /// The divide and conquer algorithm makes very mild assumptions about`
`  45:      /// floating point arithmetic. It will work on machines with a guard`
`  46:      /// digit in add/subtract, or on those binary machines without guard`
`  47:      /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or`
`  48:      /// Cray-2. It could conceivably fail on hexadecimal or decimal machines`
`  49:      /// without guard digits, but we know of none.`
`  50:      /// `
`  51:      ///</summary>`
`  52:      public class DGESDD`
`  53:      {`
`  54:      `
`  55:   `
`  56:          #region Dependencies`
`  57:          `
`  58:          DBDSDC _dbdsdc; DGEBRD _dgebrd; DGELQF _dgelqf; DGEMM _dgemm; DGEQRF _dgeqrf; DLACPY _dlacpy; DLASCL _dlascl; `
`  59:          DLASET _dlaset;DORGBR _dorgbr; DORGLQ _dorglq; DORGQR _dorgqr; DORMBR _dormbr; XERBLA _xerbla; DLAMCH _dlamch; `
`  60:          DLANGE _dlange;ILAENV _ilaenv; LSAME _lsame; `
`  61:   `
`  62:          #endregion`
`  63:   `
`  64:   `
`  65:          #region Fields`
`  66:          `
`  67:          const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool WNTQA = false; bool WNTQAS = false; `
`  68:          bool WNTQN = false;bool WNTQO = false; bool WNTQS = false; int BDSPAC = 0; int BLK = 0; int CHUNK = 0; int I = 0; `
`  69:          int IE = 0;int IERR = 0; int IL = 0; int IR = 0; int ISCL = 0; int ITAU = 0; int ITAUP = 0; int ITAUQ = 0; int IU = 0; `
`  70:          int IVT = 0;int LDWKVT = 0; int LDWRKL = 0; int LDWRKR = 0; int LDWRKU = 0; int MAXWRK = 0; int MINMN = 0; int MINWRK = 0; `
`  71:          int MNTHR = 0;int NWORK = 0; int WRKBL = 0; double ANRM = 0; double BIGNUM = 0; double EPS = 0; double SMLNUM = 0; `
`  72:          int[] IDUM = new int[1]; int offset_idum = 0; double[] DUM = new double[1]; int offset_dum = 0;`
`  73:   `
`  74:          #endregion`
`  75:   `
`  76:          public DGESDD(DBDSDC dbdsdc, DGEBRD dgebrd, DGELQF dgelqf, DGEMM dgemm, DGEQRF dgeqrf, DLACPY dlacpy, DLASCL dlascl, DLASET dlaset, DORGBR dorgbr, DORGLQ dorglq`
`  77:                        , DORGQR dorgqr, DORMBR dormbr, XERBLA xerbla, DLAMCH dlamch, DLANGE dlange, ILAENV ilaenv, LSAME lsame)`
`  78:          {`
`  79:      `
`  80:   `
`  81:              #region Set Dependencies`
`  82:              `
`  83:              this._dbdsdc = dbdsdc; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf; `
`  84:              this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq; `
`  85:              this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._dlamch = dlamch; this._dlange = dlange; `
`  86:              this._ilaenv = ilaenv;this._lsame = lsame; `
`  87:   `
`  88:              #endregion`
`  89:   `
`  90:          }`
`  91:      `
`  92:          public DGESDD()`
`  93:          {`
`  94:      `
`  95:   `
`  96:              #region Dependencies (Initialization)`
`  97:              `
`  98:              LSAME lsame = new LSAME();`
`  99:              IEEECK ieeeck = new IEEECK();`
` 100:              IPARMQ iparmq = new IPARMQ();`
` 101:              DLAMC3 dlamc3 = new DLAMC3();`
` 102:              DLASSQ dlassq = new DLASSQ();`
` 103:              DCOPY dcopy = new DCOPY();`
` 104:              XERBLA xerbla = new XERBLA();`
` 105:              DLAMRG dlamrg = new DLAMRG();`
` 106:              DLAPY2 dlapy2 = new DLAPY2();`
` 107:              DROT drot = new DROT();`
` 108:              DNRM2 dnrm2 = new DNRM2();`
` 109:              DLASD5 dlasd5 = new DLASD5();`
` 110:              DLAS2 dlas2 = new DLAS2();`
` 111:              DLASQ5 dlasq5 = new DLASQ5();`
` 112:              DLAZQ4 dlazq4 = new DLAZQ4();`
` 113:              DSCAL dscal = new DSCAL();`
` 114:              DSWAP dswap = new DSWAP();`
` 115:              DLASDT dlasdt = new DLASDT();`
` 116:              DDOT ddot = new DDOT();`
` 117:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 118:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 119:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 120:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 121:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 122:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 123:              DLANST dlanst = new DLANST(lsame, dlassq);`
` 124:              DLARTG dlartg = new DLARTG(dlamch);`
` 125:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);`
` 126:              DLACPY dlacpy = new DLACPY(lsame);`
` 127:              DLASET dlaset = new DLASET(lsame);`
` 128:              DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);`
` 129:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
` 130:              DLAED6 dlaed6 = new DLAED6(dlamch);`
` 131:              DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);`
` 132:              DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);`
` 133:              DLASD1 dlasd1 = new DLASD1(dlamrg, dlascl, dlasd2, dlasd3, xerbla);`
` 134:              DLASQ6 dlasq6 = new DLASQ6(dlamch);`
` 135:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);`
` 136:              DLASRT dlasrt = new DLASRT(lsame, xerbla);`
` 137:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);`
` 138:              DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);`
` 139:              DLASR dlasr = new DLASR(lsame, xerbla);`
` 140:              DLASV2 dlasv2 = new DLASV2(dlamch);`
` 141:              DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap`
` 142:                                         , xerbla);`
` 143:              DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);`
` 144:              DLASD0 dlasd0 = new DLASD0(dlasd1, dlasdq, dlasdt, xerbla);`
` 145:              DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);`
` 146:              DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);`
` 147:              DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);`
` 148:              DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);`
` 149:              DBDSDC dbdsdc = new DBDSDC(lsame, ilaenv, dlamch, dlanst, dcopy, dlartg, dlascl, dlasd0, dlasda, dlasdq`
` 150:                                         , dlaset, dlasr, dswap, xerbla);`
` 151:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
` 152:              DGER dger = new DGER(xerbla);`
` 153:              DLARF dlarf = new DLARF(dgemv, dger, lsame);`
` 154:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
` 155:              DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);`
` 156:              DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);`
` 157:              DGEBRD dgebrd = new DGEBRD(dgebd2, dgemm, dlabrd, xerbla, ilaenv);`
` 158:              DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);`
` 159:              DTRMM dtrmm = new DTRMM(lsame, xerbla);`
` 160:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);`
` 161:              DTRMV dtrmv = new DTRMV(lsame, xerbla);`
` 162:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);`
` 163:              DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);`
` 164:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);`
` 165:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);`
` 166:              DORGL2 dorgl2 = new DORGL2(dlarf, dscal, xerbla);`
` 167:              DORGLQ dorglq = new DORGLQ(dlarfb, dlarft, dorgl2, xerbla, ilaenv);`
` 168:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);`
` 169:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);`
` 170:              DORGBR dorgbr = new DORGBR(lsame, ilaenv, dorglq, dorgqr, xerbla);`
` 171:              DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);`
` 172:              DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);`
` 173:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);`
` 174:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);`
` 175:              DORMBR dormbr = new DORMBR(lsame, ilaenv, dormlq, dormqr, xerbla);`
` 176:              DLANGE dlange = new DLANGE(dlassq, lsame);`
` 177:   `
` 178:              #endregion`
` 179:   `
` 180:   `
` 181:              #region Set Dependencies`
` 182:              `
` 183:              this._dbdsdc = dbdsdc; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf; `
` 184:              this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq; `
` 185:              this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._dlamch = dlamch; this._dlange = dlange; `
` 186:              this._ilaenv = ilaenv;this._lsame = lsame; `
` 187:   `
` 188:              #endregion`
` 189:   `
` 190:          }`
` 191:          /// <summary>`
` 192:          /// Purpose`
` 193:          /// =======`
` 194:          /// `
` 195:          /// DGESDD computes the singular value decomposition (SVD) of a real`
` 196:          /// M-by-N matrix A, optionally computing the left and right singular`
` 197:          /// vectors.  If singular vectors are desired, it uses a`
` 198:          /// divide-and-conquer algorithm.`
` 199:          /// `
` 200:          /// The SVD is written`
` 201:          /// `
` 202:          /// A = U * SIGMA * transpose(V)`
` 203:          /// `
` 204:          /// where SIGMA is an M-by-N matrix which is zero except for its`
` 205:          /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and`
` 206:          /// V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA`
` 207:          /// are the singular values of A; they are real and non-negative, and`
` 208:          /// are returned in descending order.  The first min(m,n) columns of`
` 209:          /// U and V are the left and right singular vectors of A.`
` 210:          /// `
` 211:          /// Note that the routine returns VT = V**T, not V.`
` 212:          /// `
` 213:          /// The divide and conquer algorithm makes very mild assumptions about`
` 214:          /// floating point arithmetic. It will work on machines with a guard`
` 215:          /// digit in add/subtract, or on those binary machines without guard`
` 216:          /// digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or`
` 217:          /// Cray-2. It could conceivably fail on hexadecimal or decimal machines`
` 218:          /// without guard digits, but we know of none.`
` 219:          /// `
` 220:          ///</summary>`
` 221:          /// <param name="JOBZ">`
` 222:          /// (input) CHARACTER*1`
` 223:          /// Specifies options for computing all or part of the matrix U:`
` 224:          /// = 'A':  all M columns of U and all N rows of V**T are`
` 225:          /// returned in the arrays U and VT;`
` 226:          /// = 'S':  the first min(M,N) columns of U and the first`
` 227:          /// min(M,N) rows of V**T are returned in the arrays U`
` 228:          /// and VT;`
` 229:          /// = 'O':  If M .GE. N, the first N columns of U are overwritten`
` 230:          /// on the array A and all rows of V**T are returned in`
` 231:          /// the array VT;`
` 232:          /// otherwise, all columns of U are returned in the`
` 233:          /// array U and the first M rows of V**T are overwritten`
` 234:          /// in the array A;`
` 235:          /// = 'N':  no columns of U or rows of V**T are computed.`
` 236:          ///</param>`
` 237:          /// <param name="M">`
` 238:          /// (input) INTEGER`
` 239:          /// The number of rows of the input matrix A.  M .GE. 0.`
` 240:          ///</param>`
` 241:          /// <param name="N">`
` 242:          /// (input) INTEGER`
` 243:          /// The number of columns of the input matrix A.  N .GE. 0.`
` 244:          ///</param>`
` 245:          /// <param name="A">`
` 246:          /// = U * SIGMA * transpose(V)`
` 247:          ///</param>`
` 248:          /// <param name="LDA">`
` 249:          /// (input) INTEGER`
` 250:          /// The leading dimension of the array A.  LDA .GE. max(1,M).`
` 251:          ///</param>`
` 252:          /// <param name="S">`
` 253:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))`
` 254:          /// The singular values of A, sorted so that S(i) .GE. S(i+1).`
` 255:          ///</param>`
` 256:          /// <param name="U">`
` 257:          /// (output) DOUBLE PRECISION array, dimension (LDU,UCOL)`
` 258:          /// UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M .LT. N;`
` 259:          /// UCOL = min(M,N) if JOBZ = 'S'.`
` 260:          /// If JOBZ = 'A' or JOBZ = 'O' and M .LT. N, U contains the M-by-M`
` 261:          /// orthogonal matrix U;`
` 262:          /// if JOBZ = 'S', U contains the first min(M,N) columns of U`
` 263:          /// (the left singular vectors, stored columnwise);`
` 264:          /// if JOBZ = 'O' and M .GE. N, or JOBZ = 'N', U is not referenced.`
` 265:          ///</param>`
` 266:          /// <param name="LDU">`
` 267:          /// (input) INTEGER`
` 268:          /// The leading dimension of the array U.  LDU .GE. 1; if`
` 269:          /// JOBZ = 'S' or 'A' or JOBZ = 'O' and M .LT. N, LDU .GE. M.`
` 270:          ///</param>`
` 271:          /// <param name="VT">`
` 272:          /// (output) DOUBLE PRECISION array, dimension (LDVT,N)`
` 273:          /// If JOBZ = 'A' or JOBZ = 'O' and M .GE. N, VT contains the`
` 274:          /// N-by-N orthogonal matrix V**T;`
` 275:          /// if JOBZ = 'S', VT contains the first min(M,N) rows of`
` 276:          /// V**T (the right singular vectors, stored rowwise);`
` 277:          /// if JOBZ = 'O' and M .LT. N, or JOBZ = 'N', VT is not referenced.`
` 278:          ///</param>`
` 279:          /// <param name="LDVT">`
` 280:          /// (input) INTEGER`
` 281:          /// The leading dimension of the array VT.  LDVT .GE. 1; if`
` 282:          /// JOBZ = 'A' or JOBZ = 'O' and M .GE. N, LDVT .GE. N;`
` 283:          /// if JOBZ = 'S', LDVT .GE. min(M,N).`
` 284:          ///</param>`
` 285:          /// <param name="WORK">`
` 286:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 287:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK;`
` 288:          ///</param>`
` 289:          /// <param name="LWORK">`
` 290:          /// (input) INTEGER`
` 291:          /// The dimension of the array WORK. LWORK .GE. 1.`
` 292:          /// If JOBZ = 'N',`
` 293:          /// LWORK .GE. 3*min(M,N) + max(max(M,N),7*min(M,N)).`
` 294:          /// If JOBZ = 'O',`
` 295:          /// LWORK .GE. 3*min(M,N)*min(M,N) + `
` 296:          /// max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).`
` 297:          /// If JOBZ = 'S' or 'A'`
` 298:          /// LWORK .GE. 3*min(M,N)*min(M,N) +`
` 299:          /// max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).`
` 300:          /// For good performance, LWORK should generally be larger.`
` 301:          /// If LWORK = -1 but other input arguments are legal, WORK(1)`
` 302:          /// returns the optimal LWORK.`
` 303:          ///</param>`
` 304:          /// <param name="IWORK">`
` 305:          /// (workspace) INTEGER array, dimension (8*min(M,N))`
` 306:          ///</param>`
` 307:          /// <param name="INFO">`
` 308:          /// (output) INTEGER`
` 309:          /// = 0:  successful exit.`
` 310:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 311:          /// .GT. 0:  DBDSDC did not converge, updating process failed.`
` 312:          ///</param>`
` 313:          public void Run(string JOBZ, int M, int N, ref double[] A, int offset_a, int LDA, ref double[] S, int offset_s`
` 314:                           , ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT, ref double[] WORK, int offset_work, int LWORK`
` 315:                           , ref int[] IWORK, int offset_iwork, ref int INFO)`
` 316:          {`
` 317:   `
` 318:              #region Array Index Correction`
` 319:              `
` 320:               int o_a = -1 - LDA + offset_a;  int o_s = -1 + offset_s;  int o_u = -1 - LDU + offset_u; `
` 321:               int o_vt = -1 - LDVT + offset_vt; int o_work = -1 + offset_work;  int o_iwork = -1 + offset_iwork; `
` 322:   `
` 323:              #endregion`
` 324:   `
` 325:   `
` 326:              #region Strings`
` 327:              `
` 328:              JOBZ = JOBZ.Substring(0, 1);  `
` 329:   `
` 330:              #endregion`
` 331:   `
` 332:   `
` 333:              #region Prolog`
` 334:              `
` 335:              // *`
` 336:              // *  -- LAPACK driver routine (version 3.1) --`
` 337:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 338:              // *     November 2006`
` 339:              // *`
` 340:              // *     .. Scalar Arguments ..`
` 341:              // *     ..`
` 342:              // *     .. Array Arguments ..`
` 343:              // *     ..`
` 344:              // *`
` 345:              // *  Purpose`
` 346:              // *  =======`
` 347:              // *`
` 348:              // *  DGESDD computes the singular value decomposition (SVD) of a real`
` 349:              // *  M-by-N matrix A, optionally computing the left and right singular`
` 350:              // *  vectors.  If singular vectors are desired, it uses a`
` 351:              // *  divide-and-conquer algorithm.`
` 352:              // *`
` 353:              // *  The SVD is written`
` 354:              // *`
` 355:              // *       A = U * SIGMA * transpose(V)`
` 356:              // *`
` 357:              // *  where SIGMA is an M-by-N matrix which is zero except for its`
` 358:              // *  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and`
` 359:              // *  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA`
` 360:              // *  are the singular values of A; they are real and non-negative, and`
` 361:              // *  are returned in descending order.  The first min(m,n) columns of`
` 362:              // *  U and V are the left and right singular vectors of A.`
` 363:              // *`
` 364:              // *  Note that the routine returns VT = V**T, not V.`
` 365:              // *`
` 366:              // *  The divide and conquer algorithm makes very mild assumptions about`
` 367:              // *  floating point arithmetic. It will work on machines with a guard`
` 368:              // *  digit in add/subtract, or on those binary machines without guard`
` 369:              // *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or`
` 370:              // *  Cray-2. It could conceivably fail on hexadecimal or decimal machines`
` 371:              // *  without guard digits, but we know of none.`
` 372:              // *`
` 373:              // *  Arguments`
` 374:              // *  =========`
` 375:              // *`
` 376:              // *  JOBZ    (input) CHARACTER*1`
` 377:              // *          Specifies options for computing all or part of the matrix U:`
` 378:              // *          = 'A':  all M columns of U and all N rows of V**T are`
` 379:              // *                  returned in the arrays U and VT;`
` 380:              // *          = 'S':  the first min(M,N) columns of U and the first`
` 381:              // *                  min(M,N) rows of V**T are returned in the arrays U`
` 382:              // *                  and VT;`
` 383:              // *          = 'O':  If M >= N, the first N columns of U are overwritten`
` 384:              // *                  on the array A and all rows of V**T are returned in`
` 385:              // *                  the array VT;`
` 386:              // *                  otherwise, all columns of U are returned in the`
` 387:              // *                  array U and the first M rows of V**T are overwritten`
` 388:              // *                  in the array A;`
` 389:              // *          = 'N':  no columns of U or rows of V**T are computed.`
` 390:              // *`
` 391:              // *  M       (input) INTEGER`
` 392:              // *          The number of rows of the input matrix A.  M >= 0.`
` 393:              // *`
` 394:              // *  N       (input) INTEGER`
` 395:              // *          The number of columns of the input matrix A.  N >= 0.`
` 396:              // *`
` 397:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 398:              // *          On entry, the M-by-N matrix A.`
` 399:              // *          On exit,`
` 400:              // *          if JOBZ = 'O',  A is overwritten with the first N columns`
` 401:              // *                          of U (the left singular vectors, stored`
` 402:              // *                          columnwise) if M >= N;`
` 403:              // *                          A is overwritten with the first M rows`
` 404:              // *                          of V**T (the right singular vectors, stored`
` 405:              // *                          rowwise) otherwise.`
` 406:              // *          if JOBZ .ne. 'O', the contents of A are destroyed.`
` 407:              // *`
` 408:              // *  LDA     (input) INTEGER`
` 409:              // *          The leading dimension of the array A.  LDA >= max(1,M).`
` 410:              // *`
` 411:              // *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))`
` 412:              // *          The singular values of A, sorted so that S(i) >= S(i+1).`
` 413:              // *`
` 414:              // *  U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)`
` 415:              // *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;`
` 416:              // *          UCOL = min(M,N) if JOBZ = 'S'.`
` 417:              // *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M`
` 418:              // *          orthogonal matrix U;`
` 419:              // *          if JOBZ = 'S', U contains the first min(M,N) columns of U`
` 420:              // *          (the left singular vectors, stored columnwise);`
` 421:              // *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.`
` 422:              // *`
` 423:              // *  LDU     (input) INTEGER`
` 424:              // *          The leading dimension of the array U.  LDU >= 1; if`
` 425:              // *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.`
` 426:              // *`
` 427:              // *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)`
` 428:              // *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the`
` 429:              // *          N-by-N orthogonal matrix V**T;`
` 430:              // *          if JOBZ = 'S', VT contains the first min(M,N) rows of`
` 431:              // *          V**T (the right singular vectors, stored rowwise);`
` 432:              // *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.`
` 433:              // *`
` 434:              // *  LDVT    (input) INTEGER`
` 435:              // *          The leading dimension of the array VT.  LDVT >= 1; if`
` 436:              // *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;`
` 437:              // *          if JOBZ = 'S', LDVT >= min(M,N).`
` 438:              // *`
` 439:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 440:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;`
` 441:              // *`
` 442:              // *  LWORK   (input) INTEGER`
` 443:              // *          The dimension of the array WORK. LWORK >= 1.`
` 444:              // *          If JOBZ = 'N',`
` 445:              // *            LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).`
` 446:              // *          If JOBZ = 'O',`
` 447:              // *            LWORK >= 3*min(M,N)*min(M,N) + `
` 448:              // *                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).`
` 449:              // *          If JOBZ = 'S' or 'A'`
` 450:              // *            LWORK >= 3*min(M,N)*min(M,N) +`
` 451:              // *                     max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).`
` 452:              // *          For good performance, LWORK should generally be larger.`
` 453:              // *          If LWORK = -1 but other input arguments are legal, WORK(1)`
` 454:              // *          returns the optimal LWORK.`
` 455:              // *`
` 456:              // *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))`
` 457:              // *`
` 458:              // *  INFO    (output) INTEGER`
` 459:              // *          = 0:  successful exit.`
` 460:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 461:              // *          > 0:  DBDSDC did not converge, updating process failed.`
` 462:              // *`
` 463:              // *  Further Details`
` 464:              // *  ===============`
` 465:              // *`
` 466:              // *  Based on contributions by`
` 467:              // *     Ming Gu and Huan Ren, Computer Science Division, University of`
` 468:              // *     California at Berkeley, USA`
` 469:              // *`
` 470:              // *  =====================================================================`
` 471:              // *`
` 472:              // *     .. Parameters ..`
` 473:              // *     ..`
` 474:              // *     .. Local Scalars ..`
` 475:              // *     ..`
` 476:              // *     .. Local Arrays ..`
` 477:              // *     ..`
` 478:              // *     .. External Subroutines ..`
` 479:              // *     ..`
` 480:              // *     .. External Functions ..`
` 481:              // *     ..`
` 482:              // *     .. Intrinsic Functions ..`
` 483:              //      INTRINSIC          INT, MAX, MIN, SQRT;`
` 484:              // *     ..`
` 485:              // *     .. Executable Statements ..`
` 486:              // *`
` 487:              // *     Test the input arguments`
` 488:              // *`
` 489:   `
` 490:              #endregion`
` 491:   `
` 492:   `
` 493:              #region Body`
` 494:              `
` 495:              INFO = 0;`
` 496:              MINMN = Math.Min(M, N);`
` 497:              WNTQA = this._lsame.Run(JOBZ, "A");`
` 498:              WNTQS = this._lsame.Run(JOBZ, "S");`
` 499:              WNTQAS = WNTQA || WNTQS;`
` 500:              WNTQO = this._lsame.Run(JOBZ, "O");`
` 501:              WNTQN = this._lsame.Run(JOBZ, "N");`
` 502:              LQUERY = (LWORK ==  - 1);`
` 503:              // *`
` 504:              if (!(WNTQA || WNTQS || WNTQO || WNTQN))`
` 505:              {`
` 506:                  INFO =  - 1;`
` 507:              }`
` 508:              else`
` 509:              {`
` 510:                  if (M < 0)`
` 511:                  {`
` 512:                      INFO =  - 2;`
` 513:                  }`
` 514:                  else`
` 515:                  {`
` 516:                      if (N < 0)`
` 517:                      {`
` 518:                          INFO =  - 3;`
` 519:                      }`
` 520:                      else`
` 521:                      {`
` 522:                          if (LDA < Math.Max(1, M))`
` 523:                          {`
` 524:                              INFO =  - 5;`
` 525:                          }`
` 526:                          else`
` 527:                          {`
` 528:                              if (LDU < 1 || (WNTQAS && LDU < M) || (WNTQO && M < N && LDU < M))`
` 529:                              {`
` 530:                                  INFO =  - 8;`
` 531:                              }`
` 532:                              else`
` 533:                              {`
` 534:                                  if (LDVT < 1 || (WNTQA && LDVT < N) || (WNTQS && LDVT < MINMN) || (WNTQO && M >= N && LDVT < N))`
` 535:                                  {`
` 536:                                      INFO =  - 10;`
` 537:                                  }`
` 538:                              }`
` 539:                          }`
` 540:                      }`
` 541:                  }`
` 542:              }`
` 543:              // *`
` 544:              // *     Compute workspace`
` 545:              // *      (Note: Comments in the code beginning "Workspace:" describe the`
` 546:              // *       minimal amount of workspace needed at that point in the code,`
` 547:              // *       as well as the preferred amount for good performance.`
` 548:              // *       NB refers to the optimal block size for the immediately`
` 549:              // *       following subroutine, as returned by ILAENV.)`
` 550:              // *`
` 551:              if (INFO == 0)`
` 552:              {`
` 553:                  MINWRK = 1;`
` 554:                  MAXWRK = 1;`
` 555:                  if (M >= N && MINMN > 0)`
` 556:                  {`
` 557:                      // *`
` 558:                      // *           Compute space needed for DBDSDC`
` 559:                      // *`
` 560:                      MNTHR = Convert.ToInt32(Math.Truncate(MINMN * 11.0E0 / 6.0E0));`
` 561:                      if (WNTQN)`
` 562:                      {`
` 563:                          BDSPAC = 7 * N;`
` 564:                      }`
` 565:                      else`
` 566:                      {`
` 567:                          BDSPAC = 3 * N * N + 4 * N;`
` 568:                      }`
` 569:                      if (M >= MNTHR)`
` 570:                      {`
` 571:                          if (WNTQN)`
` 572:                          {`
` 573:                              // *`
` 574:                              // *                 Path 1 (M much larger than N, JOBZ='N')`
` 575:                              // *`
` 576:                              WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);`
` 577:                              WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));`
` 578:                              MAXWRK = Math.Max(WRKBL, BDSPAC + N);`
` 579:                              MINWRK = BDSPAC + N;`
` 580:                          }`
` 581:                          else`
` 582:                          {`
` 583:                              if (WNTQO)`
` 584:                              {`
` 585:                                  // *`
` 586:                                  // *                 Path 2 (M much larger than N, JOBZ='O')`
` 587:                                  // *`
` 588:                                  WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);`
` 589:                                  WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));`
` 590:                                  WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));`
` 591:                                  WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", N, N, N,  - 1));`
` 592:                                  WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N,  - 1));`
` 593:                                  WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);`
` 594:                                  MAXWRK = WRKBL + 2 * N * N;`
` 595:                                  MINWRK = BDSPAC + 2 * N * N + 3 * N;`
` 596:                              }`
` 597:                              else`
` 598:                              {`
` 599:                                  if (WNTQS)`
` 600:                                  {`
` 601:                                      // *`
` 602:                                      // *                 Path 3 (M much larger than N, JOBZ='S')`
` 603:                                      // *`
` 604:                                      WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);`
` 605:                                      WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));`
` 606:                                      WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));`
` 607:                                      WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", N, N, N,  - 1));`
` 608:                                      WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N,  - 1));`
` 609:                                      WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);`
` 610:                                      MAXWRK = WRKBL + N * N;`
` 611:                                      MINWRK = BDSPAC + N * N + 3 * N;`
` 612:                                  }`
` 613:                                  else`
` 614:                                  {`
` 615:                                      if (WNTQA)`
` 616:                                      {`
` 617:                                          // *`
` 618:                                          // *                 Path 4 (M much larger than N, JOBZ='A')`
` 619:                                          // *`
` 620:                                          WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);`
` 621:                                          WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N,  - 1));`
` 622:                                          WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));`
` 623:                                          WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", N, N, N,  - 1));`
` 624:                                          WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N,  - 1));`
` 625:                                          WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);`
` 626:                                          MAXWRK = WRKBL + N * N;`
` 627:                                          MINWRK = BDSPAC + N * N + 3 * N;`
` 628:                                      }`
` 629:                                  }`
` 630:                              }`
` 631:                          }`
` 632:                      }`
` 633:                      else`
` 634:                      {`
` 635:                          // *`
` 636:                          // *              Path 5 (M at least N, but not much larger)`
` 637:                          // *`
` 638:                          WRKBL = 3 * N + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N,  - 1,  - 1);`
` 639:                          if (WNTQN)`
` 640:                          {`
` 641:                              MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * N);`
` 642:                              MINWRK = 3 * N + Math.Max(M, BDSPAC);`
` 643:                          }`
` 644:                          else`
` 645:                          {`
` 646:                              if (WNTQO)`
` 647:                              {`
` 648:                                  WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", M, N, N,  - 1));`
` 649:                                  WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N,  - 1));`
` 650:                                  WRKBL = Math.Max(WRKBL, BDSPAC + 3 * N);`
` 651:                                  MAXWRK = WRKBL + M * N;`
` 652:                                  MINWRK = 3 * N + Math.Max(M, N * N + BDSPAC);`
` 653:                              }`
` 654:                              else`
` 655:                              {`
` 656:                                  if (WNTQS)`
` 657:                                  {`
` 658:                                      WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "QLN", M, N, N,  - 1));`
` 659:                                      WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N,  - 1));`
` 660:                                      MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * N);`
` 661:                                      MINWRK = 3 * N + Math.Max(M, BDSPAC);`
` 662:                                  }`
` 663:                                  else`
` 664:                                  {`
` 665:                                      if (WNTQA)`
` 666:                                      {`
` 667:                                          WRKBL = Math.Max(WRKBL, 3 * N + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N,  - 1));`
` 668:                                          WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, N,  - 1));`
` 669:                                          MAXWRK = Math.Max(MAXWRK, BDSPAC + 3 * N);`
` 670:                                          MINWRK = 3 * N + Math.Max(M, BDSPAC);`
` 671:                                      }`
` 672:                                  }`
` 673:                              }`
` 674:                          }`
` 675:                      }`
` 676:                  }`
` 677:                  else`
` 678:                  {`
` 679:                      if (MINMN > 0)`
` 680:                      {`
` 681:                          // *`
` 682:                          // *           Compute space needed for DBDSDC`
` 683:                          // *`
` 684:                          MNTHR = Convert.ToInt32(Math.Truncate(MINMN * 11.0E0 / 6.0E0));`
` 685:                          if (WNTQN)`
` 686:                          {`
` 687:                              BDSPAC = 7 * M;`
` 688:                          }`
` 689:                          else`
` 690:                          {`
` 691:                              BDSPAC = 3 * M * M + 4 * M;`
` 692:                          }`
` 693:                          if (N >= MNTHR)`
` 694:                          {`
` 695:                              if (WNTQN)`
` 696:                              {`
` 697:                                  // *`
` 698:                                  // *                 Path 1t (N much larger than M, JOBZ='N')`
` 699:                                  // *`
` 700:                                  WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);`
` 701:                                  WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));`
` 702:                                  MAXWRK = Math.Max(WRKBL, BDSPAC + M);`
` 703:                                  MINWRK = BDSPAC + M;`
` 704:                              }`
` 705:                              else`
` 706:                              {`
` 707:                                  if (WNTQO)`
` 708:                                  {`
` 709:                                      // *`
` 710:                                      // *                 Path 2t (N much larger than M, JOBZ='O')`
` 711:                                      // *`
` 712:                                      WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);`
` 713:                                      WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));`
` 714:                                      WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));`
` 715:                                      WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, M,  - 1));`
` 716:                                      WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, M, M,  - 1));`
` 717:                                      WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 718:                                      MAXWRK = WRKBL + 2 * M * M;`
` 719:                                      MINWRK = BDSPAC + 2 * M * M + 3 * M;`
` 720:                                  }`
` 721:                                  else`
` 722:                                  {`
` 723:                                      if (WNTQS)`
` 724:                                      {`
` 725:                                          // *`
` 726:                                          // *                 Path 3t (N much larger than M, JOBZ='S')`
` 727:                                          // *`
` 728:                                          WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);`
` 729:                                          WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));`
` 730:                                          WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));`
` 731:                                          WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, M,  - 1));`
` 732:                                          WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, M, M,  - 1));`
` 733:                                          WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 734:                                          MAXWRK = WRKBL + M * M;`
` 735:                                          MINWRK = BDSPAC + M * M + 3 * M;`
` 736:                                      }`
` 737:                                      else`
` 738:                                      {`
` 739:                                          if (WNTQA)`
` 740:                                          {`
` 741:                                              // *`
` 742:                                              // *                 Path 4t (N much larger than M, JOBZ='A')`
` 743:                                              // *`
` 744:                                              WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);`
` 745:                                              WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M,  - 1));`
` 746:                                              WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));`
` 747:                                              WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, M,  - 1));`
` 748:                                              WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, M, M,  - 1));`
` 749:                                              WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 750:                                              MAXWRK = WRKBL + M * M;`
` 751:                                              MINWRK = BDSPAC + M * M + 3 * M;`
` 752:                                          }`
` 753:                                      }`
` 754:                                  }`
` 755:                              }`
` 756:                          }`
` 757:                          else`
` 758:                          {`
` 759:                              // *`
` 760:                              // *              Path 5t (N greater than M, but not much larger)`
` 761:                              // *`
` 762:                              WRKBL = 3 * M + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N,  - 1,  - 1);`
` 763:                              if (WNTQN)`
` 764:                              {`
` 765:                                  MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 766:                                  MINWRK = 3 * M + Math.Max(N, BDSPAC);`
` 767:                              }`
` 768:                              else`
` 769:                              {`
` 770:                                  if (WNTQO)`
` 771:                                  {`
` 772:                                      WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N,  - 1));`
` 773:                                      WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, N, M,  - 1));`
` 774:                                      WRKBL = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 775:                                      MAXWRK = WRKBL + M * N;`
` 776:                                      MINWRK = 3 * M + Math.Max(N, M * M + BDSPAC);`
` 777:                                  }`
` 778:                                  else`
` 779:                                  {`
` 780:                                      if (WNTQS)`
` 781:                                      {`
` 782:                                          WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N,  - 1));`
` 783:                                          WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", M, N, M,  - 1));`
` 784:                                          MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 785:                                          MINWRK = 3 * M + Math.Max(N, BDSPAC);`
` 786:                                      }`
` 787:                                      else`
` 788:                                      {`
` 789:                                          if (WNTQA)`
` 790:                                          {`
` 791:                                              WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "QLN", M, M, N,  - 1));`
` 792:                                              WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORMBR", "PRT", N, N, M,  - 1));`
` 793:                                              MAXWRK = Math.Max(WRKBL, BDSPAC + 3 * M);`
` 794:                                              MINWRK = 3 * M + Math.Max(N, BDSPAC);`
` 795:                                          }`
` 796:                                      }`
` 797:                                  }`
` 798:                              }`
` 799:                          }`
` 800:                      }`
` 801:                  }`
` 802:                  MAXWRK = Math.Max(MAXWRK, MINWRK);`
` 803:                  WORK[1 + o_work] = MAXWRK;`
` 804:                  // *`
` 805:                  if (LWORK < MINWRK && !LQUERY)`
` 806:                  {`
` 807:                      INFO =  - 12;`
` 808:                  }`
` 809:              }`
` 810:              // *`
` 811:              if (INFO != 0)`
` 812:              {`
` 813:                  this._xerbla.Run("DGESDD",  - INFO);`
` 814:                  return;`
` 815:              }`
` 816:              else`
` 817:              {`
` 818:                  if (LQUERY)`
` 819:                  {`
` 820:                      return;`
` 821:                  }`
` 822:              }`
` 823:              // *`
` 824:              // *     Quick return if possible`
` 825:              // *`
` 826:              if (M == 0 || N == 0)`
` 827:              {`
` 828:                  return;`
` 829:              }`
` 830:              // *`
` 831:              // *     Get machine constants`
` 832:              // *`
` 833:              EPS = this._dlamch.Run("P");`
` 834:              SMLNUM = Math.Sqrt(this._dlamch.Run("S")) / EPS;`
` 835:              BIGNUM = ONE / SMLNUM;`
` 836:              // *`
` 837:              // *     Scale A if max element outside range [SMLNUM,BIGNUM]`
` 838:              // *`
` 839:              ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref DUM, offset_dum);`
` 840:              ISCL = 0;`
` 841:              if (ANRM > ZERO && ANRM < SMLNUM)`
` 842:              {`
` 843:                  ISCL = 1;`
` 844:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M`
` 845:                                   , N, ref A, offset_a, LDA, ref IERR);`
` 846:              }`
` 847:              else`
` 848:              {`
` 849:                  if (ANRM > BIGNUM)`
` 850:                  {`
` 851:                      ISCL = 1;`
` 852:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M`
` 853:                                       , N, ref A, offset_a, LDA, ref IERR);`
` 854:                  }`
` 855:              }`
` 856:              // *`
` 857:              if (M >= N)`
` 858:              {`
` 859:                  // *`
` 860:                  // *        A has at least as many rows as columns. If A has sufficiently`
` 861:                  // *        more rows than columns, first reduce using the QR`
` 862:                  // *        decomposition (if sufficient workspace available)`
` 863:                  // *`
` 864:                  if (M >= MNTHR)`
` 865:                  {`
` 866:                      // *`
` 867:                      if (WNTQN)`
` 868:                      {`
` 869:                          // *`
` 870:                          // *              Path 1 (M much larger than N, JOBZ='N')`
` 871:                          // *              No singular vectors to be computed`
` 872:                          // *`
` 873:                          ITAU = 1;`
` 874:                          NWORK = ITAU + N;`
` 875:                          // *`
` 876:                          // *              Compute A=Q*R`
` 877:                          // *              (Workspace: need 2*N, prefer N+N*NB)`
` 878:                          // *`
` 879:                          this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
` 880:                                           , LWORK - NWORK + 1, ref IERR);`
` 881:                          // *`
` 882:                          // *              Zero out below R`
` 883:                          // *`
` 884:                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a`
` 885:                                           , LDA);`
` 886:                          IE = 1;`
` 887:                          ITAUQ = IE + N;`
` 888:                          ITAUP = ITAUQ + N;`
` 889:                          NWORK = ITAUP + N;`
` 890:                          // *`
` 891:                          // *              Bidiagonalize R in A`
` 892:                          // *              (Workspace: need 4*N, prefer 3*N+2*N*NB)`
` 893:                          // *`
` 894:                          this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
` 895:                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
` 896:                          NWORK = IE + N;`
` 897:                          // *`
` 898:                          // *              Perform bidiagonal SVD, computing singular values only`
` 899:                          // *              (Workspace: need N+BDSPAC)`
` 900:                          // *`
` 901:                          this._dbdsdc.Run("U", "N", N, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum`
` 902:                                           , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
` 903:                                           , ref IWORK, offset_iwork, ref INFO);`
` 904:                          // *`
` 905:                      }`
` 906:                      else`
` 907:                      {`
` 908:                          if (WNTQO)`
` 909:                          {`
` 910:                              // *`
` 911:                              // *              Path 2 (M much larger than N, JOBZ = 'O')`
` 912:                              // *              N left singular vectors to be overwritten on A and`
` 913:                              // *              N right singular vectors to be computed in VT`
` 914:                              // *`
` 915:                              IR = 1;`
` 916:                              // *`
` 917:                              // *              WORK(IR) is LDWRKR by N`
` 918:                              // *`
` 919:                              if (LWORK >= LDA * N + N * N + 3 * N + BDSPAC)`
` 920:                              {`
` 921:                                  LDWRKR = LDA;`
` 922:                              }`
` 923:                              else`
` 924:                              {`
` 925:                                  LDWRKR = (LWORK - N * N - 3 * N - BDSPAC) / N;`
` 926:                              }`
` 927:                              ITAU = IR + LDWRKR * N;`
` 928:                              NWORK = ITAU + N;`
` 929:                              // *`
` 930:                              // *              Compute A=Q*R`
` 931:                              // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
` 932:                              // *`
` 933:                              this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
` 934:                                               , LWORK - NWORK + 1, ref IERR);`
` 935:                              // *`
` 936:                              // *              Copy R to WORK(IR), zeroing out below it`
` 937:                              // *`
` 938:                              this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work`
` 939:                                               , LDWRKR);`
` 940:                              this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work`
` 941:                                               , LDWRKR);`
` 942:                              // *`
` 943:                              // *              Generate Q in A`
` 944:                              // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
` 945:                              // *`
` 946:                              this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work`
` 947:                                               , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
` 948:                              IE = ITAU;`
` 949:                              ITAUQ = IE + N;`
` 950:                              ITAUP = ITAUQ + N;`
` 951:                              NWORK = ITAUP + N;`
` 952:                              // *`
` 953:                              // *              Bidiagonalize R in VT, copying result to WORK(IR)`
` 954:                              // *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)`
` 955:                              // *`
` 956:                              this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work`
` 957:                                               , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
` 958:                              // *`
` 959:                              // *              WORK(IU) is N by N`
` 960:                              // *`
` 961:                              IU = NWORK;`
` 962:                              NWORK = IU + N * N;`
` 963:                              // *`
` 964:                              // *              Perform bidiagonal SVD, computing left singular vectors`
` 965:                              // *              of bidiagonal matrix in WORK(IU) and computing right`
` 966:                              // *              singular vectors of bidiagonal matrix in VT`
` 967:                              // *              (Workspace: need N+N*N+BDSPAC)`
` 968:                              // *`
` 969:                              this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref WORK, IU + o_work`
` 970:                                               , N, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
` 971:                                               , ref IWORK, offset_iwork, ref INFO);`
` 972:                              // *`
` 973:                              // *              Overwrite WORK(IU) by left singular vectors of R`
` 974:                              // *              and VT by right singular vectors of R`
` 975:                              // *              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)`
` 976:                              // *`
` 977:                              this._dormbr.Run("Q", "L", "N", N, N, N`
` 978:                                               , ref WORK, IR + o_work, LDWRKR, WORK, ITAUQ + o_work, ref WORK, IU + o_work, N, ref WORK, NWORK + o_work`
` 979:                                               , LWORK - NWORK + 1, ref IERR);`
` 980:                              this._dormbr.Run("P", "R", "T", N, N, N`
` 981:                                               , ref WORK, IR + o_work, LDWRKR, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
` 982:                                               , LWORK - NWORK + 1, ref IERR);`
` 983:                              // *`
` 984:                              // *              Multiply Q in A by left singular vectors of R in`
` 985:                              // *              WORK(IU), storing result in WORK(IR) and copying to A`
` 986:                              // *              (Workspace: need 2*N*N, prefer N*N+M*N)`
` 987:                              // *`
` 988:                              for (I = 1; (LDWRKR >= 0) ? (I <= M) : (I >= M); I += LDWRKR)`
` 989:                              {`
` 990:                                  CHUNK = Math.Min(M - I + 1, LDWRKR);`
` 991:                                  this._dgemm.Run("N", "N", CHUNK, N, N, ONE`
` 992:                                                  , A, I+1 * LDA + o_a, LDA, WORK, IU + o_work, N, ZERO, ref WORK, IR + o_work`
` 993:                                                  , LDWRKR);`
` 994:                                  this._dlacpy.Run("F", CHUNK, N, WORK, IR + o_work, LDWRKR, ref A, I+1 * LDA + o_a`
` 995:                                                   , LDA);`
` 996:                              }`
` 997:                              // *`
` 998:                          }`
` 999:                          else`
`1000:                          {`
`1001:                              if (WNTQS)`
`1002:                              {`
`1003:                                  // *`
`1004:                                  // *              Path 3 (M much larger than N, JOBZ='S')`
`1005:                                  // *              N left singular vectors to be computed in U and`
`1006:                                  // *              N right singular vectors to be computed in VT`
`1007:                                  // *`
`1008:                                  IR = 1;`
`1009:                                  // *`
`1010:                                  // *              WORK(IR) is N by N`
`1011:                                  // *`
`1012:                                  LDWRKR = N;`
`1013:                                  ITAU = IR + LDWRKR * N;`
`1014:                                  NWORK = ITAU + N;`
`1015:                                  // *`
`1016:                                  // *              Compute A=Q*R`
`1017:                                  // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1018:                                  // *`
`1019:                                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
`1020:                                                   , LWORK - NWORK + 1, ref IERR);`
`1021:                                  // *`
`1022:                                  // *              Copy R to WORK(IR), zeroing out below it`
`1023:                                  // *`
`1024:                                  this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work`
`1025:                                                   , LDWRKR);`
`1026:                                  this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work`
`1027:                                                   , LDWRKR);`
`1028:                                  // *`
`1029:                                  // *              Generate Q in A`
`1030:                                  // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1031:                                  // *`
`1032:                                  this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work`
`1033:                                                   , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1034:                                  IE = ITAU;`
`1035:                                  ITAUQ = IE + N;`
`1036:                                  ITAUP = ITAUQ + N;`
`1037:                                  NWORK = ITAUP + N;`
`1038:                                  // *`
`1039:                                  // *              Bidiagonalize R in WORK(IR)`
`1040:                                  // *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)`
`1041:                                  // *`
`1042:                                  this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work`
`1043:                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1044:                                  // *`
`1045:                                  // *              Perform bidiagonal SVD, computing left singular vectors`
`1046:                                  // *              of bidiagoal matrix in U and computing right singular`
`1047:                                  // *              vectors of bidiagonal matrix in VT`
`1048:                                  // *              (Workspace: need N+BDSPAC)`
`1049:                                  // *`
`1050:                                  this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1051:                                                   , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1052:                                                   , ref IWORK, offset_iwork, ref INFO);`
`1053:                                  // *`
`1054:                                  // *              Overwrite U by left singular vectors of R and VT`
`1055:                                  // *              by right singular vectors of R`
`1056:                                  // *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)`
`1057:                                  // *`
`1058:                                  this._dormbr.Run("Q", "L", "N", N, N, N`
`1059:                                                   , ref WORK, IR + o_work, LDWRKR, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1060:                                                   , LWORK - NWORK + 1, ref IERR);`
`1061:                                  // *`
`1062:                                  this._dormbr.Run("P", "R", "T", N, N, N`
`1063:                                                   , ref WORK, IR + o_work, LDWRKR, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1064:                                                   , LWORK - NWORK + 1, ref IERR);`
`1065:                                  // *`
`1066:                                  // *              Multiply Q in A by left singular vectors of R in`
`1067:                                  // *              WORK(IR), storing result in U`
`1068:                                  // *              (Workspace: need N*N)`
`1069:                                  // *`
`1070:                                  this._dlacpy.Run("F", N, N, U, offset_u, LDU, ref WORK, IR + o_work`
`1071:                                                   , LDWRKR);`
`1072:                                  this._dgemm.Run("N", "N", M, N, N, ONE`
`1073:                                                  , A, offset_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref U, offset_u`
`1074:                                                  , LDU);`
`1075:                                  // *`
`1076:                              }`
`1077:                              else`
`1078:                              {`
`1079:                                  if (WNTQA)`
`1080:                                  {`
`1081:                                      // *`
`1082:                                      // *              Path 4 (M much larger than N, JOBZ='A')`
`1083:                                      // *              M left singular vectors to be computed in U and`
`1084:                                      // *              N right singular vectors to be computed in VT`
`1085:                                      // *`
`1086:                                      IU = 1;`
`1087:                                      // *`
`1088:                                      // *              WORK(IU) is N by N`
`1089:                                      // *`
`1090:                                      LDWRKU = N;`
`1091:                                      ITAU = IU + LDWRKU * N;`
`1092:                                      NWORK = ITAU + N;`
`1093:                                      // *`
`1094:                                      // *              Compute A=Q*R, copying result to U`
`1095:                                      // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1096:                                      // *`
`1097:                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
`1098:                                                       , LWORK - NWORK + 1, ref IERR);`
`1099:                                      this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u`
`1100:                                                       , LDU);`
`1101:                                      // *`
`1102:                                      // *              Generate Q in U`
`1103:                                      // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1104:                                      this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work`
`1105:                                                       , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1106:                                      // *`
`1107:                                      // *              Produce R in A, zeroing out other entries`
`1108:                                      // *`
`1109:                                      this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a`
`1110:                                                       , LDA);`
`1111:                                      IE = ITAU;`
`1112:                                      ITAUQ = IE + N;`
`1113:                                      ITAUP = ITAUQ + N;`
`1114:                                      NWORK = ITAUP + N;`
`1115:                                      // *`
`1116:                                      // *              Bidiagonalize R in A`
`1117:                                      // *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)`
`1118:                                      // *`
`1119:                                      this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
`1120:                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1121:                                      // *`
`1122:                                      // *              Perform bidiagonal SVD, computing left singular vectors`
`1123:                                      // *              of bidiagonal matrix in WORK(IU) and computing right`
`1124:                                      // *              singular vectors of bidiagonal matrix in VT`
`1125:                                      // *              (Workspace: need N+N*N+BDSPAC)`
`1126:                                      // *`
`1127:                                      this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref WORK, IU + o_work`
`1128:                                                       , N, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1129:                                                       , ref IWORK, offset_iwork, ref INFO);`
`1130:                                      // *`
`1131:                                      // *              Overwrite WORK(IU) by left singular vectors of R and VT`
`1132:                                      // *              by right singular vectors of R`
`1133:                                      // *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)`
`1134:                                      // *`
`1135:                                      this._dormbr.Run("Q", "L", "N", N, N, N`
`1136:                                                       , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, NWORK + o_work`
`1137:                                                       , LWORK - NWORK + 1, ref IERR);`
`1138:                                      this._dormbr.Run("P", "R", "T", N, N, N`
`1139:                                                       , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1140:                                                       , LWORK - NWORK + 1, ref IERR);`
`1141:                                      // *`
`1142:                                      // *              Multiply Q in U by left singular vectors of R in`
`1143:                                      // *              WORK(IU), storing result in A`
`1144:                                      // *              (Workspace: need N*N)`
`1145:                                      // *`
`1146:                                      this._dgemm.Run("N", "N", M, N, N, ONE`
`1147:                                                      , U, offset_u, LDU, WORK, IU + o_work, LDWRKU, ZERO, ref A, offset_a`
`1148:                                                      , LDA);`
`1149:                                      // *`
`1150:                                      // *              Copy left singular vectors of A from A to U`
`1151:                                      // *`
`1152:                                      this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u`
`1153:                                                       , LDU);`
`1154:                                      // *`
`1155:                                  }`
`1156:                              }`
`1157:                          }`
`1158:                      }`
`1159:                      // *`
`1160:                  }`
`1161:                  else`
`1162:                  {`
`1163:                      // *`
`1164:                      // *           M .LT. MNTHR`
`1165:                      // *`
`1166:                      // *           Path 5 (M at least N, but not much larger)`
`1167:                      // *           Reduce to bidiagonal form without QR decomposition`
`1168:                      // *`
`1169:                      IE = 1;`
`1170:                      ITAUQ = IE + N;`
`1171:                      ITAUP = ITAUQ + N;`
`1172:                      NWORK = ITAUP + N;`
`1173:                      // *`
`1174:                      // *           Bidiagonalize A`
`1175:                      // *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)`
`1176:                      // *`
`1177:                      this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
`1178:                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1179:                      if (WNTQN)`
`1180:                      {`
`1181:                          // *`
`1182:                          // *              Perform bidiagonal SVD, only computing singular values`
`1183:                          // *              (Workspace: need N+BDSPAC)`
`1184:                          // *`
`1185:                          this._dbdsdc.Run("U", "N", N, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum`
`1186:                                           , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1187:                                           , ref IWORK, offset_iwork, ref INFO);`
`1188:                      }`
`1189:                      else`
`1190:                      {`
`1191:                          if (WNTQO)`
`1192:                          {`
`1193:                              IU = NWORK;`
`1194:                              if (LWORK >= M * N + 3 * N + BDSPAC)`
`1195:                              {`
`1196:                                  // *`
`1197:                                  // *                 WORK( IU ) is M by N`
`1198:                                  // *`
`1199:                                  LDWRKU = M;`
`1200:                                  NWORK = IU + LDWRKU * N;`
`1201:                                  this._dlaset.Run("F", M, N, ZERO, ZERO, ref WORK, IU + o_work`
`1202:                                                   , LDWRKU);`
`1203:                              }`
`1204:                              else`
`1205:                              {`
`1206:                                  // *`
`1207:                                  // *                 WORK( IU ) is N by N`
`1208:                                  // *`
`1209:                                  LDWRKU = N;`
`1210:                                  NWORK = IU + LDWRKU * N;`
`1211:                                  // *`
`1212:                                  // *                 WORK(IR) is LDWRKR by N`
`1213:                                  // *`
`1214:                                  IR = NWORK;`
`1215:                                  LDWRKR = (LWORK - N * N - 3 * N) / N;`
`1216:                              }`
`1217:                              NWORK = IU + LDWRKU * N;`
`1218:                              // *`
`1219:                              // *              Perform bidiagonal SVD, computing left singular vectors`
`1220:                              // *              of bidiagonal matrix in WORK(IU) and computing right`
`1221:                              // *              singular vectors of bidiagonal matrix in VT`
`1222:                              // *              (Workspace: need N+N*N+BDSPAC)`
`1223:                              // *`
`1224:                              this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref WORK, IU + o_work`
`1225:                                               , LDWRKU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1226:                                               , ref IWORK, offset_iwork, ref INFO);`
`1227:                              // *`
`1228:                              // *              Overwrite VT by right singular vectors of A`
`1229:                              // *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1230:                              // *`
`1231:                              this._dormbr.Run("P", "R", "T", N, N, N`
`1232:                                               , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1233:                                               , LWORK - NWORK + 1, ref IERR);`
`1234:                              // *`
`1235:                              if (LWORK >= M * N + 3 * N + BDSPAC)`
`1236:                              {`
`1237:                                  // *`
`1238:                                  // *                 Overwrite WORK(IU) by left singular vectors of A`
`1239:                                  // *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1240:                                  // *`
`1241:                                  this._dormbr.Run("Q", "L", "N", M, N, N`
`1242:                                                   , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, NWORK + o_work`
`1243:                                                   , LWORK - NWORK + 1, ref IERR);`
`1244:                                  // *`
`1245:                                  // *                 Copy left singular vectors of A from WORK(IU) to A`
`1246:                                  // *`
`1247:                                  this._dlacpy.Run("F", M, N, WORK, IU + o_work, LDWRKU, ref A, offset_a`
`1248:                                                   , LDA);`
`1249:                              }`
`1250:                              else`
`1251:                              {`
`1252:                                  // *`
`1253:                                  // *                 Generate Q in A`
`1254:                                  // *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)`
`1255:                                  // *`
`1256:                                  this._dorgbr.Run("Q", M, N, N, ref A, offset_a, LDA`
`1257:                                                   , WORK, ITAUQ + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1258:                                  // *`
`1259:                                  // *                 Multiply Q in A by left singular vectors of`
`1260:                                  // *                 bidiagonal matrix in WORK(IU), storing result in`
`1261:                                  // *                 WORK(IR) and copying to A`
`1262:                                  // *                 (Workspace: need 2*N*N, prefer N*N+M*N)`
`1263:                                  // *`
`1264:                                  for (I = 1; (LDWRKR >= 0) ? (I <= M) : (I >= M); I += LDWRKR)`
`1265:                                  {`
`1266:                                      CHUNK = Math.Min(M - I + 1, LDWRKR);`
`1267:                                      this._dgemm.Run("N", "N", CHUNK, N, N, ONE`
`1268:                                                      , A, I+1 * LDA + o_a, LDA, WORK, IU + o_work, LDWRKU, ZERO, ref WORK, IR + o_work`
`1269:                                                      , LDWRKR);`
`1270:                                      this._dlacpy.Run("F", CHUNK, N, WORK, IR + o_work, LDWRKR, ref A, I+1 * LDA + o_a`
`1271:                                                       , LDA);`
`1272:                                  }`
`1273:                              }`
`1274:                              // *`
`1275:                          }`
`1276:                          else`
`1277:                          {`
`1278:                              if (WNTQS)`
`1279:                              {`
`1280:                                  // *`
`1281:                                  // *              Perform bidiagonal SVD, computing left singular vectors`
`1282:                                  // *              of bidiagonal matrix in U and computing right singular`
`1283:                                  // *              vectors of bidiagonal matrix in VT`
`1284:                                  // *              (Workspace: need N+BDSPAC)`
`1285:                                  // *`
`1286:                                  this._dlaset.Run("F", M, N, ZERO, ZERO, ref U, offset_u`
`1287:                                                   , LDU);`
`1288:                                  this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1289:                                                   , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1290:                                                   , ref IWORK, offset_iwork, ref INFO);`
`1291:                                  // *`
`1292:                                  // *              Overwrite U by left singular vectors of A and VT`
`1293:                                  // *              by right singular vectors of A`
`1294:                                  // *              (Workspace: need 3*N, prefer 2*N+N*NB)`
`1295:                                  // *`
`1296:                                  this._dormbr.Run("Q", "L", "N", M, N, N`
`1297:                                                   , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1298:                                                   , LWORK - NWORK + 1, ref IERR);`
`1299:                                  this._dormbr.Run("P", "R", "T", N, N, N`
`1300:                                                   , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1301:                                                   , LWORK - NWORK + 1, ref IERR);`
`1302:                              }`
`1303:                              else`
`1304:                              {`
`1305:                                  if (WNTQA)`
`1306:                                  {`
`1307:                                      // *`
`1308:                                      // *              Perform bidiagonal SVD, computing left singular vectors`
`1309:                                      // *              of bidiagonal matrix in U and computing right singular`
`1310:                                      // *              vectors of bidiagonal matrix in VT`
`1311:                                      // *              (Workspace: need N+BDSPAC)`
`1312:                                      // *`
`1313:                                      this._dlaset.Run("F", M, M, ZERO, ZERO, ref U, offset_u`
`1314:                                                       , LDU);`
`1315:                                      this._dbdsdc.Run("U", "I", N, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1316:                                                       , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1317:                                                       , ref IWORK, offset_iwork, ref INFO);`
`1318:                                      // *`
`1319:                                      // *              Set the right corner of U to identity matrix`
`1320:                                      // *`
`1321:                                      if (M > N)`
`1322:                                      {`
`1323:                                          this._dlaset.Run("F", M - N, M - N, ZERO, ONE, ref U, N + 1+(N + 1) * LDU + o_u`
`1324:                                                           , LDU);`
`1325:                                      }`
`1326:                                      // *`
`1327:                                      // *              Overwrite U by left singular vectors of A and VT`
`1328:                                      // *              by right singular vectors of A`
`1329:                                      // *              (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)`
`1330:                                      // *`
`1331:                                      this._dormbr.Run("Q", "L", "N", M, M, N`
`1332:                                                       , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1333:                                                       , LWORK - NWORK + 1, ref IERR);`
`1334:                                      this._dormbr.Run("P", "R", "T", N, N, M`
`1335:                                                       , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1336:                                                       , LWORK - NWORK + 1, ref IERR);`
`1337:                                  }`
`1338:                              }`
`1339:                          }`
`1340:                      }`
`1341:                      // *`
`1342:                  }`
`1343:                  // *`
`1344:              }`
`1345:              else`
`1346:              {`
`1347:                  // *`
`1348:                  // *        A has more columns than rows. If A has sufficiently more`
`1349:                  // *        columns than rows, first reduce using the LQ decomposition (if`
`1350:                  // *        sufficient workspace available)`
`1351:                  // *`
`1352:                  if (N >= MNTHR)`
`1353:                  {`
`1354:                      // *`
`1355:                      if (WNTQN)`
`1356:                      {`
`1357:                          // *`
`1358:                          // *              Path 1t (N much larger than M, JOBZ='N')`
`1359:                          // *              No singular vectors to be computed`
`1360:                          // *`
`1361:                          ITAU = 1;`
`1362:                          NWORK = ITAU + M;`
`1363:                          // *`
`1364:                          // *              Compute A=L*Q`
`1365:                          // *              (Workspace: need 2*M, prefer M+M*NB)`
`1366:                          // *`
`1367:                          this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
`1368:                                           , LWORK - NWORK + 1, ref IERR);`
`1369:                          // *`
`1370:                          // *              Zero out above L`
`1371:                          // *`
`1372:                          this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a`
`1373:                                           , LDA);`
`1374:                          IE = 1;`
`1375:                          ITAUQ = IE + M;`
`1376:                          ITAUP = ITAUQ + M;`
`1377:                          NWORK = ITAUP + M;`
`1378:                          // *`
`1379:                          // *              Bidiagonalize L in A`
`1380:                          // *              (Workspace: need 4*M, prefer 3*M+2*M*NB)`
`1381:                          // *`
`1382:                          this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
`1383:                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1384:                          NWORK = IE + M;`
`1385:                          // *`
`1386:                          // *              Perform bidiagonal SVD, computing singular values only`
`1387:                          // *              (Workspace: need M+BDSPAC)`
`1388:                          // *`
`1389:                          this._dbdsdc.Run("U", "N", M, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum`
`1390:                                           , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1391:                                           , ref IWORK, offset_iwork, ref INFO);`
`1392:                          // *`
`1393:                      }`
`1394:                      else`
`1395:                      {`
`1396:                          if (WNTQO)`
`1397:                          {`
`1398:                              // *`
`1399:                              // *              Path 2t (N much larger than M, JOBZ='O')`
`1400:                              // *              M right singular vectors to be overwritten on A and`
`1401:                              // *              M left singular vectors to be computed in U`
`1402:                              // *`
`1403:                              IVT = 1;`
`1404:                              // *`
`1405:                              // *              IVT is M by M`
`1406:                              // *`
`1407:                              IL = IVT + M * M;`
`1408:                              if (LWORK >= M * N + M * M + 3 * M + BDSPAC)`
`1409:                              {`
`1410:                                  // *`
`1411:                                  // *                 WORK(IL) is M by N`
`1412:                                  // *`
`1413:                                  LDWRKL = M;`
`1414:                                  CHUNK = N;`
`1415:                              }`
`1416:                              else`
`1417:                              {`
`1418:                                  LDWRKL = M;`
`1419:                                  CHUNK = (LWORK - M * M) / M;`
`1420:                              }`
`1421:                              ITAU = IL + LDWRKL * M;`
`1422:                              NWORK = ITAU + M;`
`1423:                              // *`
`1424:                              // *              Compute A=L*Q`
`1425:                              // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1426:                              // *`
`1427:                              this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
`1428:                                               , LWORK - NWORK + 1, ref IERR);`
`1429:                              // *`
`1430:                              // *              Copy L to WORK(IL), zeroing about above it`
`1431:                              // *`
`1432:                              this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IL + o_work`
`1433:                                               , LDWRKL);`
`1434:                              this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IL + LDWRKL + o_work`
`1435:                                               , LDWRKL);`
`1436:                              // *`
`1437:                              // *              Generate Q in A`
`1438:                              // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1439:                              // *`
`1440:                              this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work`
`1441:                                               , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1442:                              IE = ITAU;`
`1443:                              ITAUQ = IE + M;`
`1444:                              ITAUP = ITAUQ + M;`
`1445:                              NWORK = ITAUP + M;`
`1446:                              // *`
`1447:                              // *              Bidiagonalize L in WORK(IL)`
`1448:                              // *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)`
`1449:                              // *`
`1450:                              this._dgebrd.Run(M, M, ref WORK, IL + o_work, LDWRKL, ref S, offset_s, ref WORK, IE + o_work`
`1451:                                               , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1452:                              // *`
`1453:                              // *              Perform bidiagonal SVD, computing left singular vectors`
`1454:                              // *              of bidiagonal matrix in U, and computing right singular`
`1455:                              // *              vectors of bidiagonal matrix in WORK(IVT)`
`1456:                              // *              (Workspace: need M+M*M+BDSPAC)`
`1457:                              // *`
`1458:                              this._dbdsdc.Run("U", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1459:                                               , LDU, ref WORK, IVT + o_work, M, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1460:                                               , ref IWORK, offset_iwork, ref INFO);`
`1461:                              // *`
`1462:                              // *              Overwrite U by left singular vectors of L and WORK(IVT)`
`1463:                              // *              by right singular vectors of L`
`1464:                              // *              (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)`
`1465:                              // *`
`1466:                              this._dormbr.Run("Q", "L", "N", M, M, M`
`1467:                                               , ref WORK, IL + o_work, LDWRKL, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1468:                                               , LWORK - NWORK + 1, ref IERR);`
`1469:                              this._dormbr.Run("P", "R", "T", M, M, M`
`1470:                                               , ref WORK, IL + o_work, LDWRKL, WORK, ITAUP + o_work, ref WORK, IVT + o_work, M, ref WORK, NWORK + o_work`
`1471:                                               , LWORK - NWORK + 1, ref IERR);`
`1472:                              // *`
`1473:                              // *              Multiply right singular vectors of L in WORK(IVT) by Q`
`1474:                              // *              in A, storing result in WORK(IL) and copying to A`
`1475:                              // *              (Workspace: need 2*M*M, prefer M*M+M*N)`
`1476:                              // *`
`1477:                              for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)`
`1478:                              {`
`1479:                                  BLK = Math.Min(N - I + 1, CHUNK);`
`1480:                                  this._dgemm.Run("N", "N", M, BLK, M, ONE`
`1481:                                                  , WORK, IVT + o_work, M, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IL + o_work`
`1482:                                                  , LDWRKL);`
`1483:                                  this._dlacpy.Run("F", M, BLK, WORK, IL + o_work, LDWRKL, ref A, 1+I * LDA + o_a`
`1484:                                                   , LDA);`
`1485:                              }`
`1486:                              // *`
`1487:                          }`
`1488:                          else`
`1489:                          {`
`1490:                              if (WNTQS)`
`1491:                              {`
`1492:                                  // *`
`1493:                                  // *              Path 3t (N much larger than M, JOBZ='S')`
`1494:                                  // *              M right singular vectors to be computed in VT and`
`1495:                                  // *              M left singular vectors to be computed in U`
`1496:                                  // *`
`1497:                                  IL = 1;`
`1498:                                  // *`
`1499:                                  // *              WORK(IL) is M by M`
`1500:                                  // *`
`1501:                                  LDWRKL = M;`
`1502:                                  ITAU = IL + LDWRKL * M;`
`1503:                                  NWORK = ITAU + M;`
`1504:                                  // *`
`1505:                                  // *              Compute A=L*Q`
`1506:                                  // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1507:                                  // *`
`1508:                                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
`1509:                                                   , LWORK - NWORK + 1, ref IERR);`
`1510:                                  // *`
`1511:                                  // *              Copy L to WORK(IL), zeroing out above it`
`1512:                                  // *`
`1513:                                  this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IL + o_work`
`1514:                                                   , LDWRKL);`
`1515:                                  this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IL + LDWRKL + o_work`
`1516:                                                   , LDWRKL);`
`1517:                                  // *`
`1518:                                  // *              Generate Q in A`
`1519:                                  // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1520:                                  // *`
`1521:                                  this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work`
`1522:                                                   , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1523:                                  IE = ITAU;`
`1524:                                  ITAUQ = IE + M;`
`1525:                                  ITAUP = ITAUQ + M;`
`1526:                                  NWORK = ITAUP + M;`
`1527:                                  // *`
`1528:                                  // *              Bidiagonalize L in WORK(IU), copying result to U`
`1529:                                  // *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)`
`1530:                                  // *`
`1531:                                  this._dgebrd.Run(M, M, ref WORK, IL + o_work, LDWRKL, ref S, offset_s, ref WORK, IE + o_work`
`1532:                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1533:                                  // *`
`1534:                                  // *              Perform bidiagonal SVD, computing left singular vectors`
`1535:                                  // *              of bidiagonal matrix in U and computing right singular`
`1536:                                  // *              vectors of bidiagonal matrix in VT`
`1537:                                  // *              (Workspace: need M+BDSPAC)`
`1538:                                  // *`
`1539:                                  this._dbdsdc.Run("U", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1540:                                                   , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1541:                                                   , ref IWORK, offset_iwork, ref INFO);`
`1542:                                  // *`
`1543:                                  // *              Overwrite U by left singular vectors of L and VT`
`1544:                                  // *              by right singular vectors of L`
`1545:                                  // *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)`
`1546:                                  // *`
`1547:                                  this._dormbr.Run("Q", "L", "N", M, M, M`
`1548:                                                   , ref WORK, IL + o_work, LDWRKL, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1549:                                                   , LWORK - NWORK + 1, ref IERR);`
`1550:                                  this._dormbr.Run("P", "R", "T", M, M, M`
`1551:                                                   , ref WORK, IL + o_work, LDWRKL, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1552:                                                   , LWORK - NWORK + 1, ref IERR);`
`1553:                                  // *`
`1554:                                  // *              Multiply right singular vectors of L in WORK(IL) by`
`1555:                                  // *              Q in A, storing result in VT`
`1556:                                  // *              (Workspace: need M*M)`
`1557:                                  // *`
`1558:                                  this._dlacpy.Run("F", M, M, VT, offset_vt, LDVT, ref WORK, IL + o_work`
`1559:                                                   , LDWRKL);`
`1560:                                  this._dgemm.Run("N", "N", M, N, M, ONE`
`1561:                                                  , WORK, IL + o_work, LDWRKL, A, offset_a, LDA, ZERO, ref VT, offset_vt`
`1562:                                                  , LDVT);`
`1563:                                  // *`
`1564:                              }`
`1565:                              else`
`1566:                              {`
`1567:                                  if (WNTQA)`
`1568:                                  {`
`1569:                                      // *`
`1570:                                      // *              Path 4t (N much larger than M, JOBZ='A')`
`1571:                                      // *              N right singular vectors to be computed in VT and`
`1572:                                      // *              M left singular vectors to be computed in U`
`1573:                                      // *`
`1574:                                      IVT = 1;`
`1575:                                      // *`
`1576:                                      // *              WORK(IVT) is M by M`
`1577:                                      // *`
`1578:                                      LDWKVT = M;`
`1579:                                      ITAU = IVT + LDWKVT * M;`
`1580:                                      NWORK = ITAU + M;`
`1581:                                      // *`
`1582:                                      // *              Compute A=L*Q, copying result to VT`
`1583:                                      // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1584:                                      // *`
`1585:                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, NWORK + o_work`
`1586:                                                       , LWORK - NWORK + 1, ref IERR);`
`1587:                                      this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt`
`1588:                                                       , LDVT);`
`1589:                                      // *`
`1590:                                      // *              Generate Q in VT`
`1591:                                      // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1592:                                      // *`
`1593:                                      this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work`
`1594:                                                       , ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1595:                                      // *`
`1596:                                      // *              Produce L in A, zeroing out other entries`
`1597:                                      // *`
`1598:                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a`
`1599:                                                       , LDA);`
`1600:                                      IE = ITAU;`
`1601:                                      ITAUQ = IE + M;`
`1602:                                      ITAUP = ITAUQ + M;`
`1603:                                      NWORK = ITAUP + M;`
`1604:                                      // *`
`1605:                                      // *              Bidiagonalize L in A`
`1606:                                      // *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)`
`1607:                                      // *`
`1608:                                      this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
`1609:                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1610:                                      // *`
`1611:                                      // *              Perform bidiagonal SVD, computing left singular vectors`
`1612:                                      // *              of bidiagonal matrix in U and computing right singular`
`1613:                                      // *              vectors of bidiagonal matrix in WORK(IVT)`
`1614:                                      // *              (Workspace: need M+M*M+BDSPAC)`
`1615:                                      // *`
`1616:                                      this._dbdsdc.Run("U", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1617:                                                       , LDU, ref WORK, IVT + o_work, LDWKVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1618:                                                       , ref IWORK, offset_iwork, ref INFO);`
`1619:                                      // *`
`1620:                                      // *              Overwrite U by left singular vectors of L and WORK(IVT)`
`1621:                                      // *              by right singular vectors of L`
`1622:                                      // *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)`
`1623:                                      // *`
`1624:                                      this._dormbr.Run("Q", "L", "N", M, M, M`
`1625:                                                       , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1626:                                                       , LWORK - NWORK + 1, ref IERR);`
`1627:                                      this._dormbr.Run("P", "R", "T", M, M, M`
`1628:                                                       , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref WORK, IVT + o_work, LDWKVT, ref WORK, NWORK + o_work`
`1629:                                                       , LWORK - NWORK + 1, ref IERR);`
`1630:                                      // *`
`1631:                                      // *              Multiply right singular vectors of L in WORK(IVT) by`
`1632:                                      // *              Q in VT, storing result in A`
`1633:                                      // *              (Workspace: need M*M)`
`1634:                                      // *`
`1635:                                      this._dgemm.Run("N", "N", M, N, M, ONE`
`1636:                                                      , WORK, IVT + o_work, LDWKVT, VT, offset_vt, LDVT, ZERO, ref A, offset_a`
`1637:                                                      , LDA);`
`1638:                                      // *`
`1639:                                      // *              Copy right singular vectors of A from A to VT`
`1640:                                      // *`
`1641:                                      this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt`
`1642:                                                       , LDVT);`
`1643:                                      // *`
`1644:                                  }`
`1645:                              }`
`1646:                          }`
`1647:                      }`
`1648:                      // *`
`1649:                  }`
`1650:                  else`
`1651:                  {`
`1652:                      // *`
`1653:                      // *           N .LT. MNTHR`
`1654:                      // *`
`1655:                      // *           Path 5t (N greater than M, but not much larger)`
`1656:                      // *           Reduce to bidiagonal form without LQ decomposition`
`1657:                      // *`
`1658:                      IE = 1;`
`1659:                      ITAUQ = IE + M;`
`1660:                      ITAUP = ITAUQ + M;`
`1661:                      NWORK = ITAUP + M;`
`1662:                      // *`
`1663:                      // *           Bidiagonalize A`
`1664:                      // *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)`
`1665:                      // *`
`1666:                      this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work`
`1667:                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1668:                      if (WNTQN)`
`1669:                      {`
`1670:                          // *`
`1671:                          // *              Perform bidiagonal SVD, only computing singular values`
`1672:                          // *              (Workspace: need M+BDSPAC)`
`1673:                          // *`
`1674:                          this._dbdsdc.Run("L", "N", M, ref S, offset_s, ref WORK, IE + o_work, ref DUM, offset_dum`
`1675:                                           , 1, ref DUM, offset_dum, 1, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1676:                                           , ref IWORK, offset_iwork, ref INFO);`
`1677:                      }`
`1678:                      else`
`1679:                      {`
`1680:                          if (WNTQO)`
`1681:                          {`
`1682:                              LDWKVT = M;`
`1683:                              IVT = NWORK;`
`1684:                              if (LWORK >= M * N + 3 * M + BDSPAC)`
`1685:                              {`
`1686:                                  // *`
`1687:                                  // *                 WORK( IVT ) is M by N`
`1688:                                  // *`
`1689:                                  this._dlaset.Run("F", M, N, ZERO, ZERO, ref WORK, IVT + o_work`
`1690:                                                   , LDWKVT);`
`1691:                                  NWORK = IVT + LDWKVT * N;`
`1692:                              }`
`1693:                              else`
`1694:                              {`
`1695:                                  // *`
`1696:                                  // *                 WORK( IVT ) is M by M`
`1697:                                  // *`
`1698:                                  NWORK = IVT + LDWKVT * M;`
`1699:                                  IL = NWORK;`
`1700:                                  // *`
`1701:                                  // *                 WORK(IL) is M by CHUNK`
`1702:                                  // *`
`1703:                                  CHUNK = (LWORK - M * M - 3 * M) / M;`
`1704:                              }`
`1705:                              // *`
`1706:                              // *              Perform bidiagonal SVD, computing left singular vectors`
`1707:                              // *              of bidiagonal matrix in U and computing right singular`
`1708:                              // *              vectors of bidiagonal matrix in WORK(IVT)`
`1709:                              // *              (Workspace: need M*M+BDSPAC)`
`1710:                              // *`
`1711:                              this._dbdsdc.Run("L", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1712:                                               , LDU, ref WORK, IVT + o_work, LDWKVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1713:                                               , ref IWORK, offset_iwork, ref INFO);`
`1714:                              // *`
`1715:                              // *              Overwrite U by left singular vectors of A`
`1716:                              // *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1717:                              // *`
`1718:                              this._dormbr.Run("Q", "L", "N", M, M, N`
`1719:                                               , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1720:                                               , LWORK - NWORK + 1, ref IERR);`
`1721:                              // *`
`1722:                              if (LWORK >= M * N + 3 * M + BDSPAC)`
`1723:                              {`
`1724:                                  // *`
`1725:                                  // *                 Overwrite WORK(IVT) by left singular vectors of A`
`1726:                                  // *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1727:                                  // *`
`1728:                                  this._dormbr.Run("P", "R", "T", M, N, M`
`1729:                                                   , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref WORK, IVT + o_work, LDWKVT, ref WORK, NWORK + o_work`
`1730:                                                   , LWORK - NWORK + 1, ref IERR);`
`1731:                                  // *`
`1732:                                  // *                 Copy right singular vectors of A from WORK(IVT) to A`
`1733:                                  // *`
`1734:                                  this._dlacpy.Run("F", M, N, WORK, IVT + o_work, LDWKVT, ref A, offset_a`
`1735:                                                   , LDA);`
`1736:                              }`
`1737:                              else`
`1738:                              {`
`1739:                                  // *`
`1740:                                  // *                 Generate P**T in A`
`1741:                                  // *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)`
`1742:                                  // *`
`1743:                                  this._dorgbr.Run("P", M, N, M, ref A, offset_a, LDA`
`1744:                                                   , WORK, ITAUP + o_work, ref WORK, NWORK + o_work, LWORK - NWORK + 1, ref IERR);`
`1745:                                  // *`
`1746:                                  // *                 Multiply Q in A by right singular vectors of`
`1747:                                  // *                 bidiagonal matrix in WORK(IVT), storing result in`
`1748:                                  // *                 WORK(IL) and copying to A`
`1749:                                  // *                 (Workspace: need 2*M*M, prefer M*M+M*N)`
`1750:                                  // *`
`1751:                                  for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)`
`1752:                                  {`
`1753:                                      BLK = Math.Min(N - I + 1, CHUNK);`
`1754:                                      this._dgemm.Run("N", "N", M, BLK, M, ONE`
`1755:                                                      , WORK, IVT + o_work, LDWKVT, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IL + o_work`
`1756:                                                      , M);`
`1757:                                      this._dlacpy.Run("F", M, BLK, WORK, IL + o_work, M, ref A, 1+I * LDA + o_a`
`1758:                                                       , LDA);`
`1759:                                  }`
`1760:                              }`
`1761:                          }`
`1762:                          else`
`1763:                          {`
`1764:                              if (WNTQS)`
`1765:                              {`
`1766:                                  // *`
`1767:                                  // *              Perform bidiagonal SVD, computing left singular vectors`
`1768:                                  // *              of bidiagonal matrix in U and computing right singular`
`1769:                                  // *              vectors of bidiagonal matrix in VT`
`1770:                                  // *              (Workspace: need M+BDSPAC)`
`1771:                                  // *`
`1772:                                  this._dlaset.Run("F", M, N, ZERO, ZERO, ref VT, offset_vt`
`1773:                                                   , LDVT);`
`1774:                                  this._dbdsdc.Run("L", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1775:                                                   , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1776:                                                   , ref IWORK, offset_iwork, ref INFO);`
`1777:                                  // *`
`1778:                                  // *              Overwrite U by left singular vectors of A and VT`
`1779:                                  // *              by right singular vectors of A`
`1780:                                  // *              (Workspace: need 3*M, prefer 2*M+M*NB)`
`1781:                                  // *`
`1782:                                  this._dormbr.Run("Q", "L", "N", M, M, N`
`1783:                                                   , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1784:                                                   , LWORK - NWORK + 1, ref IERR);`
`1785:                                  this._dormbr.Run("P", "R", "T", M, N, M`
`1786:                                                   , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1787:                                                   , LWORK - NWORK + 1, ref IERR);`
`1788:                              }`
`1789:                              else`
`1790:                              {`
`1791:                                  if (WNTQA)`
`1792:                                  {`
`1793:                                      // *`
`1794:                                      // *              Perform bidiagonal SVD, computing left singular vectors`
`1795:                                      // *              of bidiagonal matrix in U and computing right singular`
`1796:                                      // *              vectors of bidiagonal matrix in VT`
`1797:                                      // *              (Workspace: need M+BDSPAC)`
`1798:                                      // *`
`1799:                                      this._dlaset.Run("F", N, N, ZERO, ZERO, ref VT, offset_vt`
`1800:                                                       , LDVT);`
`1801:                                      this._dbdsdc.Run("L", "I", M, ref S, offset_s, ref WORK, IE + o_work, ref U, offset_u`
`1802:                                                       , LDU, ref VT, offset_vt, LDVT, ref DUM, offset_dum, ref IDUM, offset_idum, ref WORK, NWORK + o_work`
`1803:                                                       , ref IWORK, offset_iwork, ref INFO);`
`1804:                                      // *`
`1805:                                      // *              Set the right corner of VT to identity matrix`
`1806:                                      // *`
`1807:                                      if (N > M)`
`1808:                                      {`
`1809:                                          this._dlaset.Run("F", N - M, N - M, ZERO, ONE, ref VT, M + 1+(M + 1) * LDVT + o_vt`
`1810:                                                           , LDVT);`
`1811:                                      }`
`1812:                                      // *`
`1813:                                      // *              Overwrite U by left singular vectors of A and VT`
`1814:                                      // *              by right singular vectors of A`
`1815:                                      // *              (Workspace: need 2*M+N, prefer 2*M+N*NB)`
`1816:                                      // *`
`1817:                                      this._dormbr.Run("Q", "L", "N", M, M, N`
`1818:                                                       , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, NWORK + o_work`
`1819:                                                       , LWORK - NWORK + 1, ref IERR);`
`1820:                                      this._dormbr.Run("P", "R", "T", N, N, M`
`1821:                                                       , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, NWORK + o_work`
`1822:                                                       , LWORK - NWORK + 1, ref IERR);`
`1823:                                  }`
`1824:                              }`
`1825:                          }`
`1826:                      }`
`1827:                      // *`
`1828:                  }`
`1829:                  // *`
`1830:              }`
`1831:              // *`
`1832:              // *     Undo scaling if necessary`
`1833:              // *`
`1834:              if (ISCL == 1)`
`1835:              {`
`1836:                  if (ANRM > BIGNUM)`
`1837:                  {`
`1838:                      this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN`
`1839:                                       , 1, ref S, offset_s, MINMN, ref IERR);`
`1840:                  }`
`1841:                  if (ANRM < SMLNUM)`
`1842:                  {`
`1843:                      this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN`
`1844:                                       , 1, ref S, offset_s, MINMN, ref IERR);`
`1845:                  }`
`1846:              }`
`1847:              // *`
`1848:              // *     Return optimal workspace in WORK(1)`
`1849:              // *`
`1850:              WORK[1 + o_work] = MAXWRK;`
`1851:              // *`
`1852:              return;`
`1853:              // *`
`1854:              // *     End of DGESDD`
`1855:              // *`
`1856:   `
`1857:              #endregion`
`1858:   `
`1859:          }`
`1860:      }`
`1861:  }`