Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   1:  #region Translated by Jose Antonio De Santiago-Castillo.
   2:   
   3:  //Translated by Jose Antonio De Santiago-Castillo. 
   4:  //E-mail:JAntonioDeSantiago@gmail.com
   5:  //Web: www.DotNumerics.com
   6:  //
   7:  //Fortran to C# Translation.
   8:  //Translated by:
   9:  //F2CSharp Version 0.71 (November 10, 2009)
  10:  //Code Optimizations: None
  11:  //
  12:  #endregion
  13:   
  14:  using System;
  15:  using DotNumerics.FortranLibrary;
  16:   
  17:  namespace DotNumerics.CSLapack
  18:  {
  19:      /// <summary>
  20:      /// -- LAPACK driver routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DGESVD computes the singular value decomposition (SVD) of a real
  27:      /// M-by-N matrix A, optionally computing the left and/or right singular
  28:      /// vectors. The SVD is written
  29:      /// 
  30:      /// A = U * SIGMA * transpose(V)
  31:      /// 
  32:      /// where SIGMA is an M-by-N matrix which is zero except for its
  33:      /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  34:      /// V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  35:      /// are the singular values of A; they are real and non-negative, and
  36:      /// are returned in descending order.  The first min(m,n) columns of
  37:      /// U and V are the left and right singular vectors of A.
  38:      /// 
  39:      /// Note that the routine returns V**T, not V.
  40:      /// 
  41:      ///</summary>
  42:      public class DGESVD
  43:      {
  44:      
  45:   
  46:          #region Dependencies
  47:          
  48:          DBDSQR _dbdsqr; DGEBRD _dgebrd; DGELQF _dgelqf; DGEMM _dgemm; DGEQRF _dgeqrf; DLACPY _dlacpy; DLASCL _dlascl; 
  49:          DLASET _dlaset;DORGBR _dorgbr; DORGLQ _dorglq; DORGQR _dorgqr; DORMBR _dormbr; XERBLA _xerbla; LSAME _lsame; 
  50:          ILAENV _ilaenv;DLAMCH _dlamch; DLANGE _dlange; 
  51:   
  52:          #endregion
  53:   
  54:   
  55:          #region Fields
  56:          
  57:          const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool WNTUA = false; bool WNTUAS = false; 
  58:          bool WNTUN = false;bool WNTUO = false; bool WNTUS = false; bool WNTVA = false; bool WNTVAS = false; bool WNTVN = false; 
  59:          bool WNTVO = false;bool WNTVS = false; int BDSPAC = 0; int BLK = 0; int CHUNK = 0; int I = 0; int IE = 0; int IERR = 0; 
  60:          int IR = 0;int ISCL = 0; int ITAU = 0; int ITAUP = 0; int ITAUQ = 0; int IU = 0; int IWORK = 0; int LDWRKR = 0; 
  61:          int LDWRKU = 0;int MAXWRK = 0; int MINMN = 0; int MINWRK = 0; int MNTHR = 0; int NCU = 0; int NCVT = 0; int NRU = 0; 
  62:          int NRVT = 0;int WRKBL = 0; double ANRM = 0; double BIGNUM = 0; double EPS = 0; double SMLNUM = 0; 
  63:          double[] DUM = new double[1]; int offset_dum = 0;
  64:   
  65:          #endregion
  66:   
  67:          public DGESVD(DBDSQR dbdsqr, DGEBRD dgebrd, DGELQF dgelqf, DGEMM dgemm, DGEQRF dgeqrf, DLACPY dlacpy, DLASCL dlascl, DLASET dlaset, DORGBR dorgbr, DORGLQ dorglq
  68:                        , DORGQR dorgqr, DORMBR dormbr, XERBLA xerbla, LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange)
  69:          {
  70:      
  71:   
  72:              #region Set Dependencies
  73:              
  74:              this._dbdsqr = dbdsqr; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf; 
  75:              this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq; 
  76:              this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._lsame = lsame; this._ilaenv = ilaenv; 
  77:              this._dlamch = dlamch;this._dlange = dlange; 
  78:   
  79:              #endregion
  80:   
  81:          }
  82:      
  83:          public DGESVD()
  84:          {
  85:      
  86:   
  87:              #region Dependencies (Initialization)
  88:              
  89:              LSAME lsame = new LSAME();
  90:              DLAMC3 dlamc3 = new DLAMC3();
  91:              DLAS2 dlas2 = new DLAS2();
  92:              DCOPY dcopy = new DCOPY();
  93:              XERBLA xerbla = new XERBLA();
  94:              DLASQ5 dlasq5 = new DLASQ5();
  95:              DLAZQ4 dlazq4 = new DLAZQ4();
  96:              IEEECK ieeeck = new IEEECK();
  97:              IPARMQ iparmq = new IPARMQ();
  98:              DROT drot = new DROT();
  99:              DSCAL dscal = new DSCAL();
 100:              DSWAP dswap = new DSWAP();
 101:              DLAPY2 dlapy2 = new DLAPY2();
 102:              DNRM2 dnrm2 = new DNRM2();
 103:              DLASSQ dlassq = new DLASSQ();
 104:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 105:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 106:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 107:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 108:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 109:              DLARTG dlartg = new DLARTG(dlamch);
 110:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 111:              DLASQ6 dlasq6 = new DLASQ6(dlamch);
 112:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);
 113:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
 114:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 115:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);
 116:              DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);
 117:              DLASR dlasr = new DLASR(lsame, xerbla);
 118:              DLASV2 dlasv2 = new DLASV2(dlamch);
 119:              DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap
 120:                                         , xerbla);
 121:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 122:              DGER dger = new DGER(xerbla);
 123:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 124:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 125:              DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);
 126:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 127:              DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);
 128:              DGEBRD dgebrd = new DGEBRD(dgebd2, dgemm, dlabrd, xerbla, ilaenv);
 129:              DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);
 130:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
 131:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
 132:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
 133:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
 134:              DGELQF dgelqf = new DGELQF(dgelq2, dlarfb, dlarft, xerbla, ilaenv);
 135:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
 136:              DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
 137:              DLACPY dlacpy = new DLACPY(lsame);
 138:              DLASET dlaset = new DLASET(lsame);
 139:              DORGL2 dorgl2 = new DORGL2(dlarf, dscal, xerbla);
 140:              DORGLQ dorglq = new DORGLQ(dlarfb, dlarft, dorgl2, xerbla, ilaenv);
 141:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
 142:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
 143:              DORGBR dorgbr = new DORGBR(lsame, ilaenv, dorglq, dorgqr, xerbla);
 144:              DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);
 145:              DORMLQ dormlq = new DORMLQ(lsame, ilaenv, dlarfb, dlarft, dorml2, xerbla);
 146:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
 147:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
 148:              DORMBR dormbr = new DORMBR(lsame, ilaenv, dormlq, dormqr, xerbla);
 149:              DLANGE dlange = new DLANGE(dlassq, lsame);
 150:   
 151:              #endregion
 152:   
 153:   
 154:              #region Set Dependencies
 155:              
 156:              this._dbdsqr = dbdsqr; this._dgebrd = dgebrd; this._dgelqf = dgelqf; this._dgemm = dgemm; this._dgeqrf = dgeqrf; 
 157:              this._dlacpy = dlacpy;this._dlascl = dlascl; this._dlaset = dlaset; this._dorgbr = dorgbr; this._dorglq = dorglq; 
 158:              this._dorgqr = dorgqr;this._dormbr = dormbr; this._xerbla = xerbla; this._lsame = lsame; this._ilaenv = ilaenv; 
 159:              this._dlamch = dlamch;this._dlange = dlange; 
 160:   
 161:              #endregion
 162:   
 163:          }
 164:          /// <summary>
 165:          /// Purpose
 166:          /// =======
 167:          /// 
 168:          /// DGESVD computes the singular value decomposition (SVD) of a real
 169:          /// M-by-N matrix A, optionally computing the left and/or right singular
 170:          /// vectors. The SVD is written
 171:          /// 
 172:          /// A = U * SIGMA * transpose(V)
 173:          /// 
 174:          /// where SIGMA is an M-by-N matrix which is zero except for its
 175:          /// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
 176:          /// V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
 177:          /// are the singular values of A; they are real and non-negative, and
 178:          /// are returned in descending order.  The first min(m,n) columns of
 179:          /// U and V are the left and right singular vectors of A.
 180:          /// 
 181:          /// Note that the routine returns V**T, not V.
 182:          /// 
 183:          ///</summary>
 184:          /// <param name="JOBU">
 185:          /// (input) CHARACTER*1
 186:          /// Specifies options for computing all or part of the matrix U:
 187:          /// = 'A':  all M columns of U are returned in array U:
 188:          /// = 'S':  the first min(m,n) columns of U (the left singular
 189:          /// vectors) are returned in the array U;
 190:          /// = 'O':  the first min(m,n) columns of U (the left singular
 191:          /// vectors) are overwritten on the array A;
 192:          /// = 'N':  no columns of U (no left singular vectors) are
 193:          /// computed.
 194:          ///</param>
 195:          /// <param name="JOBVT">
 196:          /// (input) CHARACTER*1
 197:          /// Specifies options for computing all or part of the matrix
 198:          /// V**T:
 199:          /// = 'A':  all N rows of V**T are returned in the array VT;
 200:          /// = 'S':  the first min(m,n) rows of V**T (the right singular
 201:          /// vectors) are returned in the array VT;
 202:          /// = 'O':  the first min(m,n) rows of V**T (the right singular
 203:          /// vectors) are overwritten on the array A;
 204:          /// = 'N':  no rows of V**T (no right singular vectors) are
 205:          /// computed.
 206:          /// 
 207:          /// JOBVT and JOBU cannot both be 'O'.
 208:          ///</param>
 209:          /// <param name="M">
 210:          /// (input) INTEGER
 211:          /// The number of rows of the input matrix A.  M .GE. 0.
 212:          ///</param>
 213:          /// <param name="N">
 214:          /// (input) INTEGER
 215:          /// The number of columns of the input matrix A.  N .GE. 0.
 216:          ///</param>
 217:          /// <param name="A">
 218:          /// = U * SIGMA * transpose(V)
 219:          ///</param>
 220:          /// <param name="LDA">
 221:          /// (input) INTEGER
 222:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 223:          ///</param>
 224:          /// <param name="S">
 225:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 226:          /// The singular values of A, sorted so that S(i) .GE. S(i+1).
 227:          ///</param>
 228:          /// <param name="U">
 229:          /// (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
 230:          /// (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
 231:          /// If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
 232:          /// if JOBU = 'S', U contains the first min(m,n) columns of U
 233:          /// (the left singular vectors, stored columnwise);
 234:          /// if JOBU = 'N' or 'O', U is not referenced.
 235:          ///</param>
 236:          /// <param name="LDU">
 237:          /// (input) INTEGER
 238:          /// The leading dimension of the array U.  LDU .GE. 1; if
 239:          /// JOBU = 'S' or 'A', LDU .GE. M.
 240:          ///</param>
 241:          /// <param name="VT">
 242:          /// (output) DOUBLE PRECISION array, dimension (LDVT,N)
 243:          /// If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
 244:          /// V**T;
 245:          /// if JOBVT = 'S', VT contains the first min(m,n) rows of
 246:          /// V**T (the right singular vectors, stored rowwise);
 247:          /// if JOBVT = 'N' or 'O', VT is not referenced.
 248:          ///</param>
 249:          /// <param name="LDVT">
 250:          /// (input) INTEGER
 251:          /// The leading dimension of the array VT.  LDVT .GE. 1; if
 252:          /// JOBVT = 'A', LDVT .GE. N; if JOBVT = 'S', LDVT .GE. min(M,N).
 253:          ///</param>
 254:          /// <param name="WORK">
 255:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 256:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
 257:          /// if INFO .GT. 0, WORK(2:MIN(M,N)) contains the unconverged
 258:          /// superdiagonal elements of an upper bidiagonal matrix B
 259:          /// whose diagonal is in S (not necessarily sorted). B
 260:          /// satisfies A = U * B * VT, so it has the same singular values
 261:          /// as A, and singular vectors related by U and VT.
 262:          ///</param>
 263:          /// <param name="LWORK">
 264:          /// (input) INTEGER
 265:          /// The dimension of the array WORK.
 266:          /// LWORK .GE. MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
 267:          /// For good performance, LWORK should generally be larger.
 268:          /// 
 269:          /// If LWORK = -1, then a workspace query is assumed; the routine
 270:          /// only calculates the optimal size of the WORK array, returns
 271:          /// this value as the first entry of the WORK array, and no error
 272:          /// message related to LWORK is issued by XERBLA.
 273:          ///</param>
 274:          /// <param name="INFO">
 275:          /// (output) INTEGER
 276:          /// = 0:  successful exit.
 277:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 278:          /// .GT. 0:  if DBDSQR did not converge, INFO specifies how many
 279:          /// superdiagonals of an intermediate bidiagonal form B
 280:          /// did not converge to zero. See the description of WORK
 281:          /// above for details.
 282:          ///</param>
 283:          public void Run(string JOBU, string JOBVT, int M, int N, ref double[] A, int offset_a, int LDA
 284:                           , ref double[] S, int offset_s, ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT, ref double[] WORK, int offset_work
 285:                           , int LWORK, ref int INFO)
 286:          {
 287:   
 288:              #region Array Index Correction
 289:              
 290:               int o_a = -1 - LDA + offset_a;  int o_s = -1 + offset_s;  int o_u = -1 - LDU + offset_u; 
 291:               int o_vt = -1 - LDVT + offset_vt; int o_work = -1 + offset_work; 
 292:   
 293:              #endregion
 294:   
 295:   
 296:              #region Strings
 297:              
 298:              JOBU = JOBU.Substring(0, 1);  JOBVT = JOBVT.Substring(0, 1);  
 299:   
 300:              #endregion
 301:   
 302:   
 303:              #region Prolog
 304:              
 305:              // *
 306:              // *  -- LAPACK driver routine (version 3.1) --
 307:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 308:              // *     November 2006
 309:              // *
 310:              // *     .. Scalar Arguments ..
 311:              // *     ..
 312:              // *     .. Array Arguments ..
 313:              // *     ..
 314:              // *
 315:              // *  Purpose
 316:              // *  =======
 317:              // *
 318:              // *  DGESVD computes the singular value decomposition (SVD) of a real
 319:              // *  M-by-N matrix A, optionally computing the left and/or right singular
 320:              // *  vectors. The SVD is written
 321:              // *
 322:              // *       A = U * SIGMA * transpose(V)
 323:              // *
 324:              // *  where SIGMA is an M-by-N matrix which is zero except for its
 325:              // *  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
 326:              // *  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
 327:              // *  are the singular values of A; they are real and non-negative, and
 328:              // *  are returned in descending order.  The first min(m,n) columns of
 329:              // *  U and V are the left and right singular vectors of A.
 330:              // *
 331:              // *  Note that the routine returns V**T, not V.
 332:              // *
 333:              // *  Arguments
 334:              // *  =========
 335:              // *
 336:              // *  JOBU    (input) CHARACTER*1
 337:              // *          Specifies options for computing all or part of the matrix U:
 338:              // *          = 'A':  all M columns of U are returned in array U:
 339:              // *          = 'S':  the first min(m,n) columns of U (the left singular
 340:              // *                  vectors) are returned in the array U;
 341:              // *          = 'O':  the first min(m,n) columns of U (the left singular
 342:              // *                  vectors) are overwritten on the array A;
 343:              // *          = 'N':  no columns of U (no left singular vectors) are
 344:              // *                  computed.
 345:              // *
 346:              // *  JOBVT   (input) CHARACTER*1
 347:              // *          Specifies options for computing all or part of the matrix
 348:              // *          V**T:
 349:              // *          = 'A':  all N rows of V**T are returned in the array VT;
 350:              // *          = 'S':  the first min(m,n) rows of V**T (the right singular
 351:              // *                  vectors) are returned in the array VT;
 352:              // *          = 'O':  the first min(m,n) rows of V**T (the right singular
 353:              // *                  vectors) are overwritten on the array A;
 354:              // *          = 'N':  no rows of V**T (no right singular vectors) are
 355:              // *                  computed.
 356:              // *
 357:              // *          JOBVT and JOBU cannot both be 'O'.
 358:              // *
 359:              // *  M       (input) INTEGER
 360:              // *          The number of rows of the input matrix A.  M >= 0.
 361:              // *
 362:              // *  N       (input) INTEGER
 363:              // *          The number of columns of the input matrix A.  N >= 0.
 364:              // *
 365:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 366:              // *          On entry, the M-by-N matrix A.
 367:              // *          On exit,
 368:              // *          if JOBU = 'O',  A is overwritten with the first min(m,n)
 369:              // *                          columns of U (the left singular vectors,
 370:              // *                          stored columnwise);
 371:              // *          if JOBVT = 'O', A is overwritten with the first min(m,n)
 372:              // *                          rows of V**T (the right singular vectors,
 373:              // *                          stored rowwise);
 374:              // *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
 375:              // *                          are destroyed.
 376:              // *
 377:              // *  LDA     (input) INTEGER
 378:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 379:              // *
 380:              // *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
 381:              // *          The singular values of A, sorted so that S(i) >= S(i+1).
 382:              // *
 383:              // *  U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
 384:              // *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
 385:              // *          If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
 386:              // *          if JOBU = 'S', U contains the first min(m,n) columns of U
 387:              // *          (the left singular vectors, stored columnwise);
 388:              // *          if JOBU = 'N' or 'O', U is not referenced.
 389:              // *
 390:              // *  LDU     (input) INTEGER
 391:              // *          The leading dimension of the array U.  LDU >= 1; if
 392:              // *          JOBU = 'S' or 'A', LDU >= M.
 393:              // *
 394:              // *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
 395:              // *          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
 396:              // *          V**T;
 397:              // *          if JOBVT = 'S', VT contains the first min(m,n) rows of
 398:              // *          V**T (the right singular vectors, stored rowwise);
 399:              // *          if JOBVT = 'N' or 'O', VT is not referenced.
 400:              // *
 401:              // *  LDVT    (input) INTEGER
 402:              // *          The leading dimension of the array VT.  LDVT >= 1; if
 403:              // *          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
 404:              // *
 405:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 406:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
 407:              // *          if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged
 408:              // *          superdiagonal elements of an upper bidiagonal matrix B
 409:              // *          whose diagonal is in S (not necessarily sorted). B
 410:              // *          satisfies A = U * B * VT, so it has the same singular values
 411:              // *          as A, and singular vectors related by U and VT.
 412:              // *
 413:              // *  LWORK   (input) INTEGER
 414:              // *          The dimension of the array WORK.
 415:              // *          LWORK >= MAX(1,3*MIN(M,N)+MAX(M,N),5*MIN(M,N)).
 416:              // *          For good performance, LWORK should generally be larger.
 417:              // *
 418:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 419:              // *          only calculates the optimal size of the WORK array, returns
 420:              // *          this value as the first entry of the WORK array, and no error
 421:              // *          message related to LWORK is issued by XERBLA.
 422:              // *
 423:              // *  INFO    (output) INTEGER
 424:              // *          = 0:  successful exit.
 425:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 426:              // *          > 0:  if DBDSQR did not converge, INFO specifies how many
 427:              // *                superdiagonals of an intermediate bidiagonal form B
 428:              // *                did not converge to zero. See the description of WORK
 429:              // *                above for details.
 430:              // *
 431:              // *  =====================================================================
 432:              // *
 433:              // *     .. Parameters ..
 434:              // *     ..
 435:              // *     .. Local Scalars ..
 436:              // *     ..
 437:              // *     .. Local Arrays ..
 438:              // *     ..
 439:              // *     .. External Subroutines ..
 440:              // *     ..
 441:              // *     .. External Functions ..
 442:              // *     ..
 443:              // *     .. Intrinsic Functions ..
 444:              //      INTRINSIC          MAX, MIN, SQRT;
 445:              // *     ..
 446:              // *     .. Executable Statements ..
 447:              // *
 448:              // *     Test the input arguments
 449:              // *
 450:   
 451:              #endregion
 452:   
 453:   
 454:              #region Body
 455:              
 456:              INFO = 0;
 457:              MINMN = Math.Min(M, N);
 458:              WNTUA = this._lsame.Run(JOBU, "A");
 459:              WNTUS = this._lsame.Run(JOBU, "S");
 460:              WNTUAS = WNTUA || WNTUS;
 461:              WNTUO = this._lsame.Run(JOBU, "O");
 462:              WNTUN = this._lsame.Run(JOBU, "N");
 463:              WNTVA = this._lsame.Run(JOBVT, "A");
 464:              WNTVS = this._lsame.Run(JOBVT, "S");
 465:              WNTVAS = WNTVA || WNTVS;
 466:              WNTVO = this._lsame.Run(JOBVT, "O");
 467:              WNTVN = this._lsame.Run(JOBVT, "N");
 468:              LQUERY = (LWORK ==  - 1);
 469:              // *
 470:              if (!(WNTUA || WNTUS || WNTUO || WNTUN))
 471:              {
 472:                  INFO =  - 1;
 473:              }
 474:              else
 475:              {
 476:                  if (!(WNTVA || WNTVS || WNTVO || WNTVN) || (WNTVO && WNTUO))
 477:                  {
 478:                      INFO =  - 2;
 479:                  }
 480:                  else
 481:                  {
 482:                      if (M < 0)
 483:                      {
 484:                          INFO =  - 3;
 485:                      }
 486:                      else
 487:                      {
 488:                          if (N < 0)
 489:                          {
 490:                              INFO =  - 4;
 491:                          }
 492:                          else
 493:                          {
 494:                              if (LDA < Math.Max(1, M))
 495:                              {
 496:                                  INFO =  - 6;
 497:                              }
 498:                              else
 499:                              {
 500:                                  if (LDU < 1 || (WNTUAS && LDU < M))
 501:                                  {
 502:                                      INFO =  - 9;
 503:                                  }
 504:                                  else
 505:                                  {
 506:                                      if (LDVT < 1 || (WNTVA && LDVT < N) || (WNTVS && LDVT < MINMN))
 507:                                      {
 508:                                          INFO =  - 11;
 509:                                      }
 510:                                  }
 511:                              }
 512:                          }
 513:                      }
 514:                  }
 515:              }
 516:              // *
 517:              // *     Compute workspace
 518:              // *      (Note: Comments in the code beginning "Workspace:" describe the
 519:              // *       minimal amount of workspace needed at that point in the code,
 520:              // *       as well as the preferred amount for good performance.
 521:              // *       NB refers to the optimal block size for the immediately
 522:              // *       following subroutine, as returned by ILAENV.)
 523:              // *
 524:              if (INFO == 0)
 525:              {
 526:                  MINWRK = 1;
 527:                  MAXWRK = 1;
 528:                  if (M >= N && MINMN > 0)
 529:                  {
 530:                      // *
 531:                      // *           Compute space needed for DBDSQR
 532:                      // *
 533:                      MNTHR = this._ilaenv.Run(6, "DGESVD", JOBU + JOBVT, M, N, 0, 0);
 534:                      BDSPAC = 5 * N;
 535:                      if (M >= MNTHR)
 536:                      {
 537:                          if (WNTUN)
 538:                          {
 539:                              // *
 540:                              // *                 Path 1 (M much larger than N, JOBU='N')
 541:                              // *
 542:                              MAXWRK = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 543:                              MAXWRK = Math.Max(MAXWRK, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 544:                              if (WNTVO || WNTVAS) MAXWRK = Math.Max(MAXWRK, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 545:                              MAXWRK = Math.Max(MAXWRK, BDSPAC);
 546:                              MINWRK = Math.Max(4 * N, BDSPAC);
 547:                          }
 548:                          else
 549:                          {
 550:                              if (WNTUO && WNTVN)
 551:                              {
 552:                                  // *
 553:                                  // *                 Path 2 (M much larger than N, JOBU='O', JOBVT='N')
 554:                                  // *
 555:                                  WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 556:                                  WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));
 557:                                  WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 558:                                  WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 559:                                  WRKBL = Math.Max(WRKBL, BDSPAC);
 560:                                  MAXWRK = Math.Max(N * N + WRKBL, N * N + M * N + N);
 561:                                  MINWRK = Math.Max(3 * N + M, BDSPAC);
 562:                              }
 563:                              else
 564:                              {
 565:                                  if (WNTUO && WNTVAS)
 566:                                  {
 567:                                      // *
 568:                                      // *                 Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
 569:                                      // *                 'A')
 570:                                      // *
 571:                                      WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 572:                                      WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));
 573:                                      WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 574:                                      WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 575:                                      WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 576:                                      WRKBL = Math.Max(WRKBL, BDSPAC);
 577:                                      MAXWRK = Math.Max(N * N + WRKBL, N * N + M * N + N);
 578:                                      MINWRK = Math.Max(3 * N + M, BDSPAC);
 579:                                  }
 580:                                  else
 581:                                  {
 582:                                      if (WNTUS && WNTVN)
 583:                                      {
 584:                                          // *
 585:                                          // *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
 586:                                          // *
 587:                                          WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 588:                                          WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));
 589:                                          WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 590:                                          WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 591:                                          WRKBL = Math.Max(WRKBL, BDSPAC);
 592:                                          MAXWRK = N * N + WRKBL;
 593:                                          MINWRK = Math.Max(3 * N + M, BDSPAC);
 594:                                      }
 595:                                      else
 596:                                      {
 597:                                          if (WNTUS && WNTVO)
 598:                                          {
 599:                                              // *
 600:                                              // *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
 601:                                              // *
 602:                                              WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 603:                                              WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));
 604:                                              WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 605:                                              WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 606:                                              WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 607:                                              WRKBL = Math.Max(WRKBL, BDSPAC);
 608:                                              MAXWRK = 2 * N * N + WRKBL;
 609:                                              MINWRK = Math.Max(3 * N + M, BDSPAC);
 610:                                          }
 611:                                          else
 612:                                          {
 613:                                              if (WNTUS && WNTVAS)
 614:                                              {
 615:                                                  // *
 616:                                                  // *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
 617:                                                  // *                 'A')
 618:                                                  // *
 619:                                                  WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 620:                                                  WRKBL = Math.Max(WRKBL, N + N * this._ilaenv.Run(1, "DORGQR", " ", M, N, N,  - 1));
 621:                                                  WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 622:                                                  WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 623:                                                  WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 624:                                                  WRKBL = Math.Max(WRKBL, BDSPAC);
 625:                                                  MAXWRK = N * N + WRKBL;
 626:                                                  MINWRK = Math.Max(3 * N + M, BDSPAC);
 627:                                              }
 628:                                              else
 629:                                              {
 630:                                                  if (WNTUA && WNTVN)
 631:                                                  {
 632:                                                      // *
 633:                                                      // *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
 634:                                                      // *
 635:                                                      WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 636:                                                      WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N,  - 1));
 637:                                                      WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 638:                                                      WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 639:                                                      WRKBL = Math.Max(WRKBL, BDSPAC);
 640:                                                      MAXWRK = N * N + WRKBL;
 641:                                                      MINWRK = Math.Max(3 * N + M, BDSPAC);
 642:                                                  }
 643:                                                  else
 644:                                                  {
 645:                                                      if (WNTUA && WNTVO)
 646:                                                      {
 647:                                                          // *
 648:                                                          // *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
 649:                                                          // *
 650:                                                          WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 651:                                                          WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N,  - 1));
 652:                                                          WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 653:                                                          WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 654:                                                          WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 655:                                                          WRKBL = Math.Max(WRKBL, BDSPAC);
 656:                                                          MAXWRK = 2 * N * N + WRKBL;
 657:                                                          MINWRK = Math.Max(3 * N + M, BDSPAC);
 658:                                                      }
 659:                                                      else
 660:                                                      {
 661:                                                          if (WNTUA && WNTVAS)
 662:                                                          {
 663:                                                              // *
 664:                                                              // *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
 665:                                                              // *                 'A')
 666:                                                              // *
 667:                                                              WRKBL = N + N * this._ilaenv.Run(1, "DGEQRF", " ", M, N,  - 1,  - 1);
 668:                                                              WRKBL = Math.Max(WRKBL, N + M * this._ilaenv.Run(1, "DORGQR", " ", M, M, N,  - 1));
 669:                                                              WRKBL = Math.Max(WRKBL, 3 * N + 2 * N * this._ilaenv.Run(1, "DGEBRD", " ", N, N,  - 1,  - 1));
 670:                                                              WRKBL = Math.Max(WRKBL, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", N, N, N,  - 1));
 671:                                                              WRKBL = Math.Max(WRKBL, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 672:                                                              WRKBL = Math.Max(WRKBL, BDSPAC);
 673:                                                              MAXWRK = N * N + WRKBL;
 674:                                                              MINWRK = Math.Max(3 * N + M, BDSPAC);
 675:                                                          }
 676:                                                      }
 677:                                                  }
 678:                                              }
 679:                                          }
 680:                                      }
 681:                                  }
 682:                              }
 683:                          }
 684:                      }
 685:                      else
 686:                      {
 687:                          // *
 688:                          // *              Path 10 (M at least N, but not much larger)
 689:                          // *
 690:                          MAXWRK = 3 * N + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N,  - 1,  - 1);
 691:                          if (WNTUS || WNTUO) MAXWRK = Math.Max(MAXWRK, 3 * N + N * this._ilaenv.Run(1, "DORGBR", "Q", M, N, N,  - 1));
 692:                          if (WNTUA) MAXWRK = Math.Max(MAXWRK, 3 * N + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, N,  - 1));
 693:                          if (!WNTVN) MAXWRK = Math.Max(MAXWRK, 3 * N + (N - 1) * this._ilaenv.Run(1, "DORGBR", "P", N, N, N,  - 1));
 694:                          MAXWRK = Math.Max(MAXWRK, BDSPAC);
 695:                          MINWRK = Math.Max(3 * N + M, BDSPAC);
 696:                      }
 697:                  }
 698:                  else
 699:                  {
 700:                      if (MINMN > 0)
 701:                      {
 702:                          // *
 703:                          // *           Compute space needed for DBDSQR
 704:                          // *
 705:                          MNTHR = this._ilaenv.Run(6, "DGESVD", JOBU + JOBVT, M, N, 0, 0);
 706:                          BDSPAC = 5 * M;
 707:                          if (N >= MNTHR)
 708:                          {
 709:                              if (WNTVN)
 710:                              {
 711:                                  // *
 712:                                  // *                 Path 1t(N much larger than M, JOBVT='N')
 713:                                  // *
 714:                                  MAXWRK = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 715:                                  MAXWRK = Math.Max(MAXWRK, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 716:                                  if (WNTUO || WNTUAS) MAXWRK = Math.Max(MAXWRK, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 717:                                  MAXWRK = Math.Max(MAXWRK, BDSPAC);
 718:                                  MINWRK = Math.Max(4 * M, BDSPAC);
 719:                              }
 720:                              else
 721:                              {
 722:                                  if (WNTVO && WNTUN)
 723:                                  {
 724:                                      // *
 725:                                      // *                 Path 2t(N much larger than M, JOBU='N', JOBVT='O')
 726:                                      // *
 727:                                      WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 728:                                      WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));
 729:                                      WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 730:                                      WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 731:                                      WRKBL = Math.Max(WRKBL, BDSPAC);
 732:                                      MAXWRK = Math.Max(M * M + WRKBL, M * M + M * N + M);
 733:                                      MINWRK = Math.Max(3 * M + N, BDSPAC);
 734:                                  }
 735:                                  else
 736:                                  {
 737:                                      if (WNTVO && WNTUAS)
 738:                                      {
 739:                                          // *
 740:                                          // *                 Path 3t(N much larger than M, JOBU='S' or 'A',
 741:                                          // *                 JOBVT='O')
 742:                                          // *
 743:                                          WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 744:                                          WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));
 745:                                          WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 746:                                          WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 747:                                          WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 748:                                          WRKBL = Math.Max(WRKBL, BDSPAC);
 749:                                          MAXWRK = Math.Max(M * M + WRKBL, M * M + M * N + M);
 750:                                          MINWRK = Math.Max(3 * M + N, BDSPAC);
 751:                                      }
 752:                                      else
 753:                                      {
 754:                                          if (WNTVS && WNTUN)
 755:                                          {
 756:                                              // *
 757:                                              // *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
 758:                                              // *
 759:                                              WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 760:                                              WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));
 761:                                              WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 762:                                              WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 763:                                              WRKBL = Math.Max(WRKBL, BDSPAC);
 764:                                              MAXWRK = M * M + WRKBL;
 765:                                              MINWRK = Math.Max(3 * M + N, BDSPAC);
 766:                                          }
 767:                                          else
 768:                                          {
 769:                                              if (WNTVS && WNTUO)
 770:                                              {
 771:                                                  // *
 772:                                                  // *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
 773:                                                  // *
 774:                                                  WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 775:                                                  WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));
 776:                                                  WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 777:                                                  WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 778:                                                  WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 779:                                                  WRKBL = Math.Max(WRKBL, BDSPAC);
 780:                                                  MAXWRK = 2 * M * M + WRKBL;
 781:                                                  MINWRK = Math.Max(3 * M + N, BDSPAC);
 782:                                              }
 783:                                              else
 784:                                              {
 785:                                                  if (WNTVS && WNTUAS)
 786:                                                  {
 787:                                                      // *
 788:                                                      // *                 Path 6t(N much larger than M, JOBU='S' or 'A',
 789:                                                      // *                 JOBVT='S')
 790:                                                      // *
 791:                                                      WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 792:                                                      WRKBL = Math.Max(WRKBL, M + M * this._ilaenv.Run(1, "DORGLQ", " ", M, N, M,  - 1));
 793:                                                      WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 794:                                                      WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 795:                                                      WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 796:                                                      WRKBL = Math.Max(WRKBL, BDSPAC);
 797:                                                      MAXWRK = M * M + WRKBL;
 798:                                                      MINWRK = Math.Max(3 * M + N, BDSPAC);
 799:                                                  }
 800:                                                  else
 801:                                                  {
 802:                                                      if (WNTVA && WNTUN)
 803:                                                      {
 804:                                                          // *
 805:                                                          // *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
 806:                                                          // *
 807:                                                          WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 808:                                                          WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M,  - 1));
 809:                                                          WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 810:                                                          WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 811:                                                          WRKBL = Math.Max(WRKBL, BDSPAC);
 812:                                                          MAXWRK = M * M + WRKBL;
 813:                                                          MINWRK = Math.Max(3 * M + N, BDSPAC);
 814:                                                      }
 815:                                                      else
 816:                                                      {
 817:                                                          if (WNTVA && WNTUO)
 818:                                                          {
 819:                                                              // *
 820:                                                              // *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
 821:                                                              // *
 822:                                                              WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 823:                                                              WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M,  - 1));
 824:                                                              WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 825:                                                              WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 826:                                                              WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 827:                                                              WRKBL = Math.Max(WRKBL, BDSPAC);
 828:                                                              MAXWRK = 2 * M * M + WRKBL;
 829:                                                              MINWRK = Math.Max(3 * M + N, BDSPAC);
 830:                                                          }
 831:                                                          else
 832:                                                          {
 833:                                                              if (WNTVA && WNTUAS)
 834:                                                              {
 835:                                                                  // *
 836:                                                                  // *                 Path 9t(N much larger than M, JOBU='S' or 'A',
 837:                                                                  // *                 JOBVT='A')
 838:                                                                  // *
 839:                                                                  WRKBL = M + M * this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 840:                                                                  WRKBL = Math.Max(WRKBL, M + N * this._ilaenv.Run(1, "DORGLQ", " ", N, N, M,  - 1));
 841:                                                                  WRKBL = Math.Max(WRKBL, 3 * M + 2 * M * this._ilaenv.Run(1, "DGEBRD", " ", M, M,  - 1,  - 1));
 842:                                                                  WRKBL = Math.Max(WRKBL, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "P", M, M, M,  - 1));
 843:                                                                  WRKBL = Math.Max(WRKBL, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 844:                                                                  WRKBL = Math.Max(WRKBL, BDSPAC);
 845:                                                                  MAXWRK = M * M + WRKBL;
 846:                                                                  MINWRK = Math.Max(3 * M + N, BDSPAC);
 847:                                                              }
 848:                                                          }
 849:                                                      }
 850:                                                  }
 851:                                              }
 852:                                          }
 853:                                      }
 854:                                  }
 855:                              }
 856:                          }
 857:                          else
 858:                          {
 859:                              // *
 860:                              // *              Path 10t(N greater than M, but not much larger)
 861:                              // *
 862:                              MAXWRK = 3 * M + (M + N) * this._ilaenv.Run(1, "DGEBRD", " ", M, N,  - 1,  - 1);
 863:                              if (WNTVS || WNTVO) MAXWRK = Math.Max(MAXWRK, 3 * M + M * this._ilaenv.Run(1, "DORGBR", "P", M, N, M,  - 1));
 864:                              if (WNTVA) MAXWRK = Math.Max(MAXWRK, 3 * M + N * this._ilaenv.Run(1, "DORGBR", "P", N, N, M,  - 1));
 865:                              if (!WNTUN) MAXWRK = Math.Max(MAXWRK, 3 * M + (M - 1) * this._ilaenv.Run(1, "DORGBR", "Q", M, M, M,  - 1));
 866:                              MAXWRK = Math.Max(MAXWRK, BDSPAC);
 867:                              MINWRK = Math.Max(3 * M + N, BDSPAC);
 868:                          }
 869:                      }
 870:                  }
 871:                  MAXWRK = Math.Max(MAXWRK, MINWRK);
 872:                  WORK[1 + o_work] = MAXWRK;
 873:                  // *
 874:                  if (LWORK < MINWRK && !LQUERY)
 875:                  {
 876:                      INFO =  - 13;
 877:                  }
 878:              }
 879:              // *
 880:              if (INFO != 0)
 881:              {
 882:                  this._xerbla.Run("DGESVD",  - INFO);
 883:                  return;
 884:              }
 885:              else
 886:              {
 887:                  if (LQUERY)
 888:                  {
 889:                      return;
 890:                  }
 891:              }
 892:              // *
 893:              // *     Quick return if possible
 894:              // *
 895:              if (M == 0 || N == 0)
 896:              {
 897:                  return;
 898:              }
 899:              // *
 900:              // *     Get machine constants
 901:              // *
 902:              EPS = this._dlamch.Run("P");
 903:              SMLNUM = Math.Sqrt(this._dlamch.Run("S")) / EPS;
 904:              BIGNUM = ONE / SMLNUM;
 905:              // *
 906:              // *     Scale A if max element outside range [SMLNUM,BIGNUM]
 907:              // *
 908:              ANRM = this._dlange.Run("M", M, N, A, offset_a, LDA, ref DUM, offset_dum);
 909:              ISCL = 0;
 910:              if (ANRM > ZERO && ANRM < SMLNUM)
 911:              {
 912:                  ISCL = 1;
 913:                  this._dlascl.Run("G", 0, 0, ANRM, SMLNUM, M
 914:                                   , N, ref A, offset_a, LDA, ref IERR);
 915:              }
 916:              else
 917:              {
 918:                  if (ANRM > BIGNUM)
 919:                  {
 920:                      ISCL = 1;
 921:                      this._dlascl.Run("G", 0, 0, ANRM, BIGNUM, M
 922:                                       , N, ref A, offset_a, LDA, ref IERR);
 923:                  }
 924:              }
 925:              // *
 926:              if (M >= N)
 927:              {
 928:                  // *
 929:                  // *        A has at least as many rows as columns. If A has sufficiently
 930:                  // *        more rows than columns, first reduce using the QR
 931:                  // *        decomposition (if sufficient workspace available)
 932:                  // *
 933:                  if (M >= MNTHR)
 934:                  {
 935:                      // *
 936:                      if (WNTUN)
 937:                      {
 938:                          // *
 939:                          // *              Path 1 (M much larger than N, JOBU='N')
 940:                          // *              No left singular vectors to be computed
 941:                          // *
 942:                          ITAU = 1;
 943:                          IWORK = ITAU + N;
 944:                          // *
 945:                          // *              Compute A=Q*R
 946:                          // *              (Workspace: need 2*N, prefer N+N*NB)
 947:                          // *
 948:                          this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
 949:                                           , LWORK - IWORK + 1, ref IERR);
 950:                          // *
 951:                          // *              Zero out below R
 952:                          // *
 953:                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
 954:                                           , LDA);
 955:                          IE = 1;
 956:                          ITAUQ = IE + N;
 957:                          ITAUP = ITAUQ + N;
 958:                          IWORK = ITAUP + N;
 959:                          // *
 960:                          // *              Bidiagonalize R in A
 961:                          // *              (Workspace: need 4*N, prefer 3*N+2*N*NB)
 962:                          // *
 963:                          this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
 964:                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
 965:                          NCVT = 0;
 966:                          if (WNTVO || WNTVAS)
 967:                          {
 968:                              // *
 969:                              // *                 If right singular vectors desired, generate P'.
 970:                              // *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
 971:                              // *
 972:                              this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
 973:                                               , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
 974:                              NCVT = N;
 975:                          }
 976:                          IWORK = IE + N;
 977:                          // *
 978:                          // *              Perform bidiagonal QR iteration, computing right
 979:                          // *              singular vectors of A in A if desired
 980:                          // *              (Workspace: need BDSPAC)
 981:                          // *
 982:                          this._dbdsqr.Run("U", N, NCVT, 0, 0, ref S, offset_s
 983:                                           , ref WORK, IE + o_work, ref A, offset_a, LDA, ref DUM, offset_dum, 1, ref DUM, offset_dum
 984:                                           , 1, ref WORK, IWORK + o_work, ref INFO);
 985:                          // *
 986:                          // *              If right singular vectors desired in VT, copy them there
 987:                          // *
 988:                          if (WNTVAS)
 989:                          {
 990:                              this._dlacpy.Run("F", N, N, A, offset_a, LDA, ref VT, offset_vt
 991:                                               , LDVT);
 992:                          }
 993:                          // *
 994:                      }
 995:                      else
 996:                      {
 997:                          if (WNTUO && WNTVN)
 998:                          {
 999:                              // *
1000:                              // *              Path 2 (M much larger than N, JOBU='O', JOBVT='N')
1001:                              // *              N left singular vectors to be overwritten on A and
1002:                              // *              no right singular vectors to be computed
1003:                              // *
1004:                              if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1005:                              {
1006:                                  // *
1007:                                  // *                 Sufficient workspace for a fast algorithm
1008:                                  // *
1009:                                  IR = 1;
1010:                                  if (LWORK >= Math.Max(WRKBL, LDA * N + N) + LDA * N)
1011:                                  {
1012:                                      // *
1013:                                      // *                    WORK(IU) is LDA by N, WORK(IR) is LDA by N
1014:                                      // *
1015:                                      LDWRKU = LDA;
1016:                                      LDWRKR = LDA;
1017:                                  }
1018:                                  else
1019:                                  {
1020:                                      if (LWORK >= Math.Max(WRKBL, LDA * N + N) + N * N)
1021:                                      {
1022:                                          // *
1023:                                          // *                    WORK(IU) is LDA by N, WORK(IR) is N by N
1024:                                          // *
1025:                                          LDWRKU = LDA;
1026:                                          LDWRKR = N;
1027:                                      }
1028:                                      else
1029:                                      {
1030:                                          // *
1031:                                          // *                    WORK(IU) is LDWRKU by N, WORK(IR) is N by N
1032:                                          // *
1033:                                          LDWRKU = (LWORK - N * N - N) / N;
1034:                                          LDWRKR = N;
1035:                                      }
1036:                                  }
1037:                                  ITAU = IR + LDWRKR * N;
1038:                                  IWORK = ITAU + N;
1039:                                  // *
1040:                                  // *                 Compute A=Q*R
1041:                                  // *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1042:                                  // *
1043:                                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1044:                                                   , LWORK - IWORK + 1, ref IERR);
1045:                                  // *
1046:                                  // *                 Copy R to WORK(IR) and zero out below it
1047:                                  // *
1048:                                  this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1049:                                                   , LDWRKR);
1050:                                  this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1051:                                                   , LDWRKR);
1052:                                  // *
1053:                                  // *                 Generate Q in A
1054:                                  // *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1055:                                  // *
1056:                                  this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1057:                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1058:                                  IE = ITAU;
1059:                                  ITAUQ = IE + N;
1060:                                  ITAUP = ITAUQ + N;
1061:                                  IWORK = ITAUP + N;
1062:                                  // *
1063:                                  // *                 Bidiagonalize R in WORK(IR)
1064:                                  // *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1065:                                  // *
1066:                                  this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1067:                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1068:                                  // *
1069:                                  // *                 Generate left vectors bidiagonalizing R
1070:                                  // *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1071:                                  // *
1072:                                  this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1073:                                                   , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1074:                                  IWORK = IE + N;
1075:                                  // *
1076:                                  // *                 Perform bidiagonal QR iteration, computing left
1077:                                  // *                 singular vectors of R in WORK(IR)
1078:                                  // *                 (Workspace: need N*N+BDSPAC)
1079:                                  // *
1080:                                  this._dbdsqr.Run("U", N, 0, N, 0, ref S, offset_s
1081:                                                   , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1082:                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
1083:                                  IU = IE + N;
1084:                                  // *
1085:                                  // *                 Multiply Q in A by left singular vectors of R in
1086:                                  // *                 WORK(IR), storing result in WORK(IU) and copying to A
1087:                                  // *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
1088:                                  // *
1089:                                  for (I = 1; (LDWRKU >= 0) ? (I <= M) : (I >= M); I += LDWRKU)
1090:                                  {
1091:                                      CHUNK = Math.Min(M - I + 1, LDWRKU);
1092:                                      this._dgemm.Run("N", "N", CHUNK, N, N, ONE
1093:                                                      , A, I+1 * LDA + o_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref WORK, IU + o_work
1094:                                                      , LDWRKU);
1095:                                      this._dlacpy.Run("F", CHUNK, N, WORK, IU + o_work, LDWRKU, ref A, I+1 * LDA + o_a
1096:                                                       , LDA);
1097:                                  }
1098:                                  // *
1099:                              }
1100:                              else
1101:                              {
1102:                                  // *
1103:                                  // *                 Insufficient workspace for a fast algorithm
1104:                                  // *
1105:                                  IE = 1;
1106:                                  ITAUQ = IE + N;
1107:                                  ITAUP = ITAUQ + N;
1108:                                  IWORK = ITAUP + N;
1109:                                  // *
1110:                                  // *                 Bidiagonalize A
1111:                                  // *                 (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
1112:                                  // *
1113:                                  this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1114:                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1115:                                  // *
1116:                                  // *                 Generate left vectors bidiagonalizing A
1117:                                  // *                 (Workspace: need 4*N, prefer 3*N+N*NB)
1118:                                  // *
1119:                                  this._dorgbr.Run("Q", M, N, N, ref A, offset_a, LDA
1120:                                                   , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1121:                                  IWORK = IE + N;
1122:                                  // *
1123:                                  // *                 Perform bidiagonal QR iteration, computing left
1124:                                  // *                 singular vectors of A in A
1125:                                  // *                 (Workspace: need BDSPAC)
1126:                                  // *
1127:                                  this._dbdsqr.Run("U", N, 0, M, 0, ref S, offset_s
1128:                                                   , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref A, offset_a, LDA, ref DUM, offset_dum
1129:                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
1130:                                  // *
1131:                              }
1132:                              // *
1133:                          }
1134:                          else
1135:                          {
1136:                              if (WNTUO && WNTVAS)
1137:                              {
1138:                                  // *
1139:                                  // *              Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
1140:                                  // *              N left singular vectors to be overwritten on A and
1141:                                  // *              N right singular vectors to be computed in VT
1142:                                  // *
1143:                                  if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1144:                                  {
1145:                                      // *
1146:                                      // *                 Sufficient workspace for a fast algorithm
1147:                                      // *
1148:                                      IR = 1;
1149:                                      if (LWORK >= Math.Max(WRKBL, LDA * N + N) + LDA * N)
1150:                                      {
1151:                                          // *
1152:                                          // *                    WORK(IU) is LDA by N and WORK(IR) is LDA by N
1153:                                          // *
1154:                                          LDWRKU = LDA;
1155:                                          LDWRKR = LDA;
1156:                                      }
1157:                                      else
1158:                                      {
1159:                                          if (LWORK >= Math.Max(WRKBL, LDA * N + N) + N * N)
1160:                                          {
1161:                                              // *
1162:                                              // *                    WORK(IU) is LDA by N and WORK(IR) is N by N
1163:                                              // *
1164:                                              LDWRKU = LDA;
1165:                                              LDWRKR = N;
1166:                                          }
1167:                                          else
1168:                                          {
1169:                                              // *
1170:                                              // *                    WORK(IU) is LDWRKU by N and WORK(IR) is N by N
1171:                                              // *
1172:                                              LDWRKU = (LWORK - N * N - N) / N;
1173:                                              LDWRKR = N;
1174:                                          }
1175:                                      }
1176:                                      ITAU = IR + LDWRKR * N;
1177:                                      IWORK = ITAU + N;
1178:                                      // *
1179:                                      // *                 Compute A=Q*R
1180:                                      // *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1181:                                      // *
1182:                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1183:                                                       , LWORK - IWORK + 1, ref IERR);
1184:                                      // *
1185:                                      // *                 Copy R to VT, zeroing out below it
1186:                                      // *
1187:                                      this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
1188:                                                       , LDVT);
1189:                                      if (N > 1)
1190:                                      {
1191:                                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
1192:                                                           , LDVT);
1193:                                      }
1194:                                      // *
1195:                                      // *                 Generate Q in A
1196:                                      // *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1197:                                      // *
1198:                                      this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1199:                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1200:                                      IE = ITAU;
1201:                                      ITAUQ = IE + N;
1202:                                      ITAUP = ITAUQ + N;
1203:                                      IWORK = ITAUP + N;
1204:                                      // *
1205:                                      // *                 Bidiagonalize R in VT, copying result to WORK(IR)
1206:                                      // *                 (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1207:                                      // *
1208:                                      this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
1209:                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1210:                                      this._dlacpy.Run("L", N, N, VT, offset_vt, LDVT, ref WORK, IR + o_work
1211:                                                       , LDWRKR);
1212:                                      // *
1213:                                      // *                 Generate left vectors bidiagonalizing R in WORK(IR)
1214:                                      // *                 (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1215:                                      // *
1216:                                      this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1217:                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1218:                                      // *
1219:                                      // *                 Generate right vectors bidiagonalizing R in VT
1220:                                      // *                 (Workspace: need N*N+4*N-1, prefer N*N+3*N+(N-1)*NB)
1221:                                      // *
1222:                                      this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1223:                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1224:                                      IWORK = IE + N;
1225:                                      // *
1226:                                      // *                 Perform bidiagonal QR iteration, computing left
1227:                                      // *                 singular vectors of R in WORK(IR) and computing right
1228:                                      // *                 singular vectors of R in VT
1229:                                      // *                 (Workspace: need N*N+BDSPAC)
1230:                                      // *
1231:                                      this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
1232:                                                       , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1233:                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
1234:                                      IU = IE + N;
1235:                                      // *
1236:                                      // *                 Multiply Q in A by left singular vectors of R in
1237:                                      // *                 WORK(IR), storing result in WORK(IU) and copying to A
1238:                                      // *                 (Workspace: need N*N+2*N, prefer N*N+M*N+N)
1239:                                      // *
1240:                                      for (I = 1; (LDWRKU >= 0) ? (I <= M) : (I >= M); I += LDWRKU)
1241:                                      {
1242:                                          CHUNK = Math.Min(M - I + 1, LDWRKU);
1243:                                          this._dgemm.Run("N", "N", CHUNK, N, N, ONE
1244:                                                          , A, I+1 * LDA + o_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref WORK, IU + o_work
1245:                                                          , LDWRKU);
1246:                                          this._dlacpy.Run("F", CHUNK, N, WORK, IU + o_work, LDWRKU, ref A, I+1 * LDA + o_a
1247:                                                           , LDA);
1248:                                      }
1249:                                      // *
1250:                                  }
1251:                                  else
1252:                                  {
1253:                                      // *
1254:                                      // *                 Insufficient workspace for a fast algorithm
1255:                                      // *
1256:                                      ITAU = 1;
1257:                                      IWORK = ITAU + N;
1258:                                      // *
1259:                                      // *                 Compute A=Q*R
1260:                                      // *                 (Workspace: need 2*N, prefer N+N*NB)
1261:                                      // *
1262:                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1263:                                                       , LWORK - IWORK + 1, ref IERR);
1264:                                      // *
1265:                                      // *                 Copy R to VT, zeroing out below it
1266:                                      // *
1267:                                      this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
1268:                                                       , LDVT);
1269:                                      if (N > 1)
1270:                                      {
1271:                                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
1272:                                                           , LDVT);
1273:                                      }
1274:                                      // *
1275:                                      // *                 Generate Q in A
1276:                                      // *                 (Workspace: need 2*N, prefer N+N*NB)
1277:                                      // *
1278:                                      this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1279:                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1280:                                      IE = ITAU;
1281:                                      ITAUQ = IE + N;
1282:                                      ITAUP = ITAUQ + N;
1283:                                      IWORK = ITAUP + N;
1284:                                      // *
1285:                                      // *                 Bidiagonalize R in VT
1286:                                      // *                 (Workspace: need 4*N, prefer 3*N+2*N*NB)
1287:                                      // *
1288:                                      this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
1289:                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1290:                                      // *
1291:                                      // *                 Multiply Q in A by left vectors bidiagonalizing R
1292:                                      // *                 (Workspace: need 3*N+M, prefer 3*N+M*NB)
1293:                                      // *
1294:                                      this._dormbr.Run("Q", "R", "N", M, N, N
1295:                                                       , ref VT, offset_vt, LDVT, WORK, ITAUQ + o_work, ref A, offset_a, LDA, ref WORK, IWORK + o_work
1296:                                                       , LWORK - IWORK + 1, ref IERR);
1297:                                      // *
1298:                                      // *                 Generate right vectors bidiagonalizing R in VT
1299:                                      // *                 (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1300:                                      // *
1301:                                      this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1302:                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1303:                                      IWORK = IE + N;
1304:                                      // *
1305:                                      // *                 Perform bidiagonal QR iteration, computing left
1306:                                      // *                 singular vectors of A in A and computing right
1307:                                      // *                 singular vectors of A in VT
1308:                                      // *                 (Workspace: need BDSPAC)
1309:                                      // *
1310:                                      this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
1311:                                                       , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
1312:                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
1313:                                      // *
1314:                                  }
1315:                                  // *
1316:                              }
1317:                              else
1318:                              {
1319:                                  if (WNTUS)
1320:                                  {
1321:                                      // *
1322:                                      if (WNTVN)
1323:                                      {
1324:                                          // *
1325:                                          // *                 Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1326:                                          // *                 N left singular vectors to be computed in U and
1327:                                          // *                 no right singular vectors to be computed
1328:                                          // *
1329:                                          if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1330:                                          {
1331:                                              // *
1332:                                              // *                    Sufficient workspace for a fast algorithm
1333:                                              // *
1334:                                              IR = 1;
1335:                                              if (LWORK >= WRKBL + LDA * N)
1336:                                              {
1337:                                                  // *
1338:                                                  // *                       WORK(IR) is LDA by N
1339:                                                  // *
1340:                                                  LDWRKR = LDA;
1341:                                              }
1342:                                              else
1343:                                              {
1344:                                                  // *
1345:                                                  // *                       WORK(IR) is N by N
1346:                                                  // *
1347:                                                  LDWRKR = N;
1348:                                              }
1349:                                              ITAU = IR + LDWRKR * N;
1350:                                              IWORK = ITAU + N;
1351:                                              // *
1352:                                              // *                    Compute A=Q*R
1353:                                              // *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1354:                                              // *
1355:                                              this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1356:                                                               , LWORK - IWORK + 1, ref IERR);
1357:                                              // *
1358:                                              // *                    Copy R to WORK(IR), zeroing out below it
1359:                                              // *
1360:                                              this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1361:                                                               , LDWRKR);
1362:                                              this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1363:                                                               , LDWRKR);
1364:                                              // *
1365:                                              // *                    Generate Q in A
1366:                                              // *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1367:                                              // *
1368:                                              this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1369:                                                               , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1370:                                              IE = ITAU;
1371:                                              ITAUQ = IE + N;
1372:                                              ITAUP = ITAUQ + N;
1373:                                              IWORK = ITAUP + N;
1374:                                              // *
1375:                                              // *                    Bidiagonalize R in WORK(IR)
1376:                                              // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1377:                                              // *
1378:                                              this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1379:                                                               , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1380:                                              // *
1381:                                              // *                    Generate left vectors bidiagonalizing R in WORK(IR)
1382:                                              // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1383:                                              // *
1384:                                              this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1385:                                                               , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1386:                                              IWORK = IE + N;
1387:                                              // *
1388:                                              // *                    Perform bidiagonal QR iteration, computing left
1389:                                              // *                    singular vectors of R in WORK(IR)
1390:                                              // *                    (Workspace: need N*N+BDSPAC)
1391:                                              // *
1392:                                              this._dbdsqr.Run("U", N, 0, N, 0, ref S, offset_s
1393:                                                               , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1394:                                                               , 1, ref WORK, IWORK + o_work, ref INFO);
1395:                                              // *
1396:                                              // *                    Multiply Q in A by left singular vectors of R in
1397:                                              // *                    WORK(IR), storing result in U
1398:                                              // *                    (Workspace: need N*N)
1399:                                              // *
1400:                                              this._dgemm.Run("N", "N", M, N, N, ONE
1401:                                                              , A, offset_a, LDA, WORK, IR + o_work, LDWRKR, ZERO, ref U, offset_u
1402:                                                              , LDU);
1403:                                              // *
1404:                                          }
1405:                                          else
1406:                                          {
1407:                                              // *
1408:                                              // *                    Insufficient workspace for a fast algorithm
1409:                                              // *
1410:                                              ITAU = 1;
1411:                                              IWORK = ITAU + N;
1412:                                              // *
1413:                                              // *                    Compute A=Q*R, copying result to U
1414:                                              // *                    (Workspace: need 2*N, prefer N+N*NB)
1415:                                              // *
1416:                                              this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1417:                                                               , LWORK - IWORK + 1, ref IERR);
1418:                                              this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1419:                                                               , LDU);
1420:                                              // *
1421:                                              // *                    Generate Q in U
1422:                                              // *                    (Workspace: need 2*N, prefer N+N*NB)
1423:                                              // *
1424:                                              this._dorgqr.Run(M, N, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1425:                                                               , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1426:                                              IE = ITAU;
1427:                                              ITAUQ = IE + N;
1428:                                              ITAUP = ITAUQ + N;
1429:                                              IWORK = ITAUP + N;
1430:                                              // *
1431:                                              // *                    Zero out below R in A
1432:                                              // *
1433:                                              this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1434:                                                               , LDA);
1435:                                              // *
1436:                                              // *                    Bidiagonalize R in A
1437:                                              // *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
1438:                                              // *
1439:                                              this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1440:                                                               , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1441:                                              // *
1442:                                              // *                    Multiply Q in U by left vectors bidiagonalizing R
1443:                                              // *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
1444:                                              // *
1445:                                              this._dormbr.Run("Q", "R", "N", M, N, N
1446:                                                               , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1447:                                                               , LWORK - IWORK + 1, ref IERR);
1448:                                              IWORK = IE + N;
1449:                                              // *
1450:                                              // *                    Perform bidiagonal QR iteration, computing left
1451:                                              // *                    singular vectors of A in U
1452:                                              // *                    (Workspace: need BDSPAC)
1453:                                              // *
1454:                                              this._dbdsqr.Run("U", N, 0, M, 0, ref S, offset_s
1455:                                                               , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref U, offset_u, LDU, ref DUM, offset_dum
1456:                                                               , 1, ref WORK, IWORK + o_work, ref INFO);
1457:                                              // *
1458:                                          }
1459:                                          // *
1460:                                      }
1461:                                      else
1462:                                      {
1463:                                          if (WNTVO)
1464:                                          {
1465:                                              // *
1466:                                              // *                 Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1467:                                              // *                 N left singular vectors to be computed in U and
1468:                                              // *                 N right singular vectors to be overwritten on A
1469:                                              // *
1470:                                              if (LWORK >= 2 * N * N + Math.Max(4 * N, BDSPAC))
1471:                                              {
1472:                                                  // *
1473:                                                  // *                    Sufficient workspace for a fast algorithm
1474:                                                  // *
1475:                                                  IU = 1;
1476:                                                  if (LWORK >= WRKBL + 2 * LDA * N)
1477:                                                  {
1478:                                                      // *
1479:                                                      // *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
1480:                                                      // *
1481:                                                      LDWRKU = LDA;
1482:                                                      IR = IU + LDWRKU * N;
1483:                                                      LDWRKR = LDA;
1484:                                                  }
1485:                                                  else
1486:                                                  {
1487:                                                      if (LWORK >= WRKBL + (LDA + N) * N)
1488:                                                      {
1489:                                                          // *
1490:                                                          // *                       WORK(IU) is LDA by N and WORK(IR) is N by N
1491:                                                          // *
1492:                                                          LDWRKU = LDA;
1493:                                                          IR = IU + LDWRKU * N;
1494:                                                          LDWRKR = N;
1495:                                                      }
1496:                                                      else
1497:                                                      {
1498:                                                          // *
1499:                                                          // *                       WORK(IU) is N by N and WORK(IR) is N by N
1500:                                                          // *
1501:                                                          LDWRKU = N;
1502:                                                          IR = IU + LDWRKU * N;
1503:                                                          LDWRKR = N;
1504:                                                      }
1505:                                                  }
1506:                                                  ITAU = IR + LDWRKR * N;
1507:                                                  IWORK = ITAU + N;
1508:                                                  // *
1509:                                                  // *                    Compute A=Q*R
1510:                                                  // *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1511:                                                  // *
1512:                                                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1513:                                                                   , LWORK - IWORK + 1, ref IERR);
1514:                                                  // *
1515:                                                  // *                    Copy R to WORK(IU), zeroing out below it
1516:                                                  // *
1517:                                                  this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
1518:                                                                   , LDWRKU);
1519:                                                  this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
1520:                                                                   , LDWRKU);
1521:                                                  // *
1522:                                                  // *                    Generate Q in A
1523:                                                  // *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1524:                                                  // *
1525:                                                  this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1526:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1527:                                                  IE = ITAU;
1528:                                                  ITAUQ = IE + N;
1529:                                                  ITAUP = ITAUQ + N;
1530:                                                  IWORK = ITAUP + N;
1531:                                                  // *
1532:                                                  // *                    Bidiagonalize R in WORK(IU), copying result to
1533:                                                  // *                    WORK(IR)
1534:                                                  // *                    (Workspace: need 2*N*N+4*N,
1535:                                                  // *                                prefer 2*N*N+3*N+2*N*NB)
1536:                                                  // *
1537:                                                  this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
1538:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1539:                                                  this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
1540:                                                                   , LDWRKR);
1541:                                                  // *
1542:                                                  // *                    Generate left bidiagonalizing vectors in WORK(IU)
1543:                                                  // *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
1544:                                                  // *
1545:                                                  this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
1546:                                                                   , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1547:                                                  // *
1548:                                                  // *                    Generate right bidiagonalizing vectors in WORK(IR)
1549:                                                  // *                    (Workspace: need 2*N*N+4*N-1,
1550:                                                  // *                                prefer 2*N*N+3*N+(N-1)*NB)
1551:                                                  // *
1552:                                                  this._dorgbr.Run("P", N, N, N, ref WORK, IR + o_work, LDWRKR
1553:                                                                   , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1554:                                                  IWORK = IE + N;
1555:                                                  // *
1556:                                                  // *                    Perform bidiagonal QR iteration, computing left
1557:                                                  // *                    singular vectors of R in WORK(IU) and computing
1558:                                                  // *                    right singular vectors of R in WORK(IR)
1559:                                                  // *                    (Workspace: need 2*N*N+BDSPAC)
1560:                                                  // *
1561:                                                  this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
1562:                                                                   , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
1563:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
1564:                                                  // *
1565:                                                  // *                    Multiply Q in A by left singular vectors of R in
1566:                                                  // *                    WORK(IU), storing result in U
1567:                                                  // *                    (Workspace: need N*N)
1568:                                                  // *
1569:                                                  this._dgemm.Run("N", "N", M, N, N, ONE
1570:                                                                  , A, offset_a, LDA, WORK, IU + o_work, LDWRKU, ZERO, ref U, offset_u
1571:                                                                  , LDU);
1572:                                                  // *
1573:                                                  // *                    Copy right singular vectors of R to A
1574:                                                  // *                    (Workspace: need N*N)
1575:                                                  // *
1576:                                                  this._dlacpy.Run("F", N, N, WORK, IR + o_work, LDWRKR, ref A, offset_a
1577:                                                                   , LDA);
1578:                                                  // *
1579:                                              }
1580:                                              else
1581:                                              {
1582:                                                  // *
1583:                                                  // *                    Insufficient workspace for a fast algorithm
1584:                                                  // *
1585:                                                  ITAU = 1;
1586:                                                  IWORK = ITAU + N;
1587:                                                  // *
1588:                                                  // *                    Compute A=Q*R, copying result to U
1589:                                                  // *                    (Workspace: need 2*N, prefer N+N*NB)
1590:                                                  // *
1591:                                                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1592:                                                                   , LWORK - IWORK + 1, ref IERR);
1593:                                                  this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1594:                                                                   , LDU);
1595:                                                  // *
1596:                                                  // *                    Generate Q in U
1597:                                                  // *                    (Workspace: need 2*N, prefer N+N*NB)
1598:                                                  // *
1599:                                                  this._dorgqr.Run(M, N, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1600:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1601:                                                  IE = ITAU;
1602:                                                  ITAUQ = IE + N;
1603:                                                  ITAUP = ITAUQ + N;
1604:                                                  IWORK = ITAUP + N;
1605:                                                  // *
1606:                                                  // *                    Zero out below R in A
1607:                                                  // *
1608:                                                  this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1609:                                                                   , LDA);
1610:                                                  // *
1611:                                                  // *                    Bidiagonalize R in A
1612:                                                  // *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
1613:                                                  // *
1614:                                                  this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1615:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1616:                                                  // *
1617:                                                  // *                    Multiply Q in U by left vectors bidiagonalizing R
1618:                                                  // *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
1619:                                                  // *
1620:                                                  this._dormbr.Run("Q", "R", "N", M, N, N
1621:                                                                   , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1622:                                                                   , LWORK - IWORK + 1, ref IERR);
1623:                                                  // *
1624:                                                  // *                    Generate right vectors bidiagonalizing R in A
1625:                                                  // *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1626:                                                  // *
1627:                                                  this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
1628:                                                                   , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1629:                                                  IWORK = IE + N;
1630:                                                  // *
1631:                                                  // *                    Perform bidiagonal QR iteration, computing left
1632:                                                  // *                    singular vectors of A in U and computing right
1633:                                                  // *                    singular vectors of A in A
1634:                                                  // *                    (Workspace: need BDSPAC)
1635:                                                  // *
1636:                                                  this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
1637:                                                                   , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
1638:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
1639:                                                  // *
1640:                                              }
1641:                                              // *
1642:                                          }
1643:                                          else
1644:                                          {
1645:                                              if (WNTVAS)
1646:                                              {
1647:                                                  // *
1648:                                                  // *                 Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1649:                                                  // *                         or 'A')
1650:                                                  // *                 N left singular vectors to be computed in U and
1651:                                                  // *                 N right singular vectors to be computed in VT
1652:                                                  // *
1653:                                                  if (LWORK >= N * N + Math.Max(4 * N, BDSPAC))
1654:                                                  {
1655:                                                      // *
1656:                                                      // *                    Sufficient workspace for a fast algorithm
1657:                                                      // *
1658:                                                      IU = 1;
1659:                                                      if (LWORK >= WRKBL + LDA * N)
1660:                                                      {
1661:                                                          // *
1662:                                                          // *                       WORK(IU) is LDA by N
1663:                                                          // *
1664:                                                          LDWRKU = LDA;
1665:                                                      }
1666:                                                      else
1667:                                                      {
1668:                                                          // *
1669:                                                          // *                       WORK(IU) is N by N
1670:                                                          // *
1671:                                                          LDWRKU = N;
1672:                                                      }
1673:                                                      ITAU = IU + LDWRKU * N;
1674:                                                      IWORK = ITAU + N;
1675:                                                      // *
1676:                                                      // *                    Compute A=Q*R
1677:                                                      // *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1678:                                                      // *
1679:                                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1680:                                                                       , LWORK - IWORK + 1, ref IERR);
1681:                                                      // *
1682:                                                      // *                    Copy R to WORK(IU), zeroing out below it
1683:                                                      // *
1684:                                                      this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
1685:                                                                       , LDWRKU);
1686:                                                      this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
1687:                                                                       , LDWRKU);
1688:                                                      // *
1689:                                                      // *                    Generate Q in A
1690:                                                      // *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1691:                                                      // *
1692:                                                      this._dorgqr.Run(M, N, N, ref A, offset_a, LDA, WORK, ITAU + o_work
1693:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1694:                                                      IE = ITAU;
1695:                                                      ITAUQ = IE + N;
1696:                                                      ITAUP = ITAUQ + N;
1697:                                                      IWORK = ITAUP + N;
1698:                                                      // *
1699:                                                      // *                    Bidiagonalize R in WORK(IU), copying result to VT
1700:                                                      // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1701:                                                      // *
1702:                                                      this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
1703:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1704:                                                      this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref VT, offset_vt
1705:                                                                       , LDVT);
1706:                                                      // *
1707:                                                      // *                    Generate left bidiagonalizing vectors in WORK(IU)
1708:                                                      // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1709:                                                      // *
1710:                                                      this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
1711:                                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1712:                                                      // *
1713:                                                      // *                    Generate right bidiagonalizing vectors in VT
1714:                                                      // *                    (Workspace: need N*N+4*N-1,
1715:                                                      // *                                prefer N*N+3*N+(N-1)*NB)
1716:                                                      // *
1717:                                                      this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1718:                                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1719:                                                      IWORK = IE + N;
1720:                                                      // *
1721:                                                      // *                    Perform bidiagonal QR iteration, computing left
1722:                                                      // *                    singular vectors of R in WORK(IU) and computing
1723:                                                      // *                    right singular vectors of R in VT
1724:                                                      // *                    (Workspace: need N*N+BDSPAC)
1725:                                                      // *
1726:                                                      this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
1727:                                                                       , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
1728:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
1729:                                                      // *
1730:                                                      // *                    Multiply Q in A by left singular vectors of R in
1731:                                                      // *                    WORK(IU), storing result in U
1732:                                                      // *                    (Workspace: need N*N)
1733:                                                      // *
1734:                                                      this._dgemm.Run("N", "N", M, N, N, ONE
1735:                                                                      , A, offset_a, LDA, WORK, IU + o_work, LDWRKU, ZERO, ref U, offset_u
1736:                                                                      , LDU);
1737:                                                      // *
1738:                                                  }
1739:                                                  else
1740:                                                  {
1741:                                                      // *
1742:                                                      // *                    Insufficient workspace for a fast algorithm
1743:                                                      // *
1744:                                                      ITAU = 1;
1745:                                                      IWORK = ITAU + N;
1746:                                                      // *
1747:                                                      // *                    Compute A=Q*R, copying result to U
1748:                                                      // *                    (Workspace: need 2*N, prefer N+N*NB)
1749:                                                      // *
1750:                                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1751:                                                                       , LWORK - IWORK + 1, ref IERR);
1752:                                                      this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1753:                                                                       , LDU);
1754:                                                      // *
1755:                                                      // *                    Generate Q in U
1756:                                                      // *                    (Workspace: need 2*N, prefer N+N*NB)
1757:                                                      // *
1758:                                                      this._dorgqr.Run(M, N, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1759:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1760:                                                      // *
1761:                                                      // *                    Copy R to VT, zeroing out below it
1762:                                                      // *
1763:                                                      this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
1764:                                                                       , LDVT);
1765:                                                      if (N > 1)
1766:                                                      {
1767:                                                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
1768:                                                                           , LDVT);
1769:                                                      }
1770:                                                      IE = ITAU;
1771:                                                      ITAUQ = IE + N;
1772:                                                      ITAUP = ITAUQ + N;
1773:                                                      IWORK = ITAUP + N;
1774:                                                      // *
1775:                                                      // *                    Bidiagonalize R in VT
1776:                                                      // *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
1777:                                                      // *
1778:                                                      this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
1779:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1780:                                                      // *
1781:                                                      // *                    Multiply Q in U by left bidiagonalizing vectors
1782:                                                      // *                    in VT
1783:                                                      // *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
1784:                                                      // *
1785:                                                      this._dormbr.Run("Q", "R", "N", M, N, N
1786:                                                                       , ref VT, offset_vt, LDVT, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1787:                                                                       , LWORK - IWORK + 1, ref IERR);
1788:                                                      // *
1789:                                                      // *                    Generate right bidiagonalizing vectors in VT
1790:                                                      // *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
1791:                                                      // *
1792:                                                      this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
1793:                                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1794:                                                      IWORK = IE + N;
1795:                                                      // *
1796:                                                      // *                    Perform bidiagonal QR iteration, computing left
1797:                                                      // *                    singular vectors of A in U and computing right
1798:                                                      // *                    singular vectors of A in VT
1799:                                                      // *                    (Workspace: need BDSPAC)
1800:                                                      // *
1801:                                                      this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
1802:                                                                       , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
1803:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
1804:                                                      // *
1805:                                                  }
1806:                                                  // *
1807:                                              }
1808:                                          }
1809:                                      }
1810:                                      // *
1811:                                  }
1812:                                  else
1813:                                  {
1814:                                      if (WNTUA)
1815:                                      {
1816:                                          // *
1817:                                          if (WNTVN)
1818:                                          {
1819:                                              // *
1820:                                              // *                 Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1821:                                              // *                 M left singular vectors to be computed in U and
1822:                                              // *                 no right singular vectors to be computed
1823:                                              // *
1824:                                              if (LWORK >= N * N + Math.Max(N + M, Math.Max(4 * N, BDSPAC)))
1825:                                              {
1826:                                                  // *
1827:                                                  // *                    Sufficient workspace for a fast algorithm
1828:                                                  // *
1829:                                                  IR = 1;
1830:                                                  if (LWORK >= WRKBL + LDA * N)
1831:                                                  {
1832:                                                      // *
1833:                                                      // *                       WORK(IR) is LDA by N
1834:                                                      // *
1835:                                                      LDWRKR = LDA;
1836:                                                  }
1837:                                                  else
1838:                                                  {
1839:                                                      // *
1840:                                                      // *                       WORK(IR) is N by N
1841:                                                      // *
1842:                                                      LDWRKR = N;
1843:                                                  }
1844:                                                  ITAU = IR + LDWRKR * N;
1845:                                                  IWORK = ITAU + N;
1846:                                                  // *
1847:                                                  // *                    Compute A=Q*R, copying result to U
1848:                                                  // *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
1849:                                                  // *
1850:                                                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1851:                                                                   , LWORK - IWORK + 1, ref IERR);
1852:                                                  this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1853:                                                                   , LDU);
1854:                                                  // *
1855:                                                  // *                    Copy R to WORK(IR), zeroing out below it
1856:                                                  // *
1857:                                                  this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IR + o_work
1858:                                                                   , LDWRKR);
1859:                                                  this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IR + 1 + o_work
1860:                                                                   , LDWRKR);
1861:                                                  // *
1862:                                                  // *                    Generate Q in U
1863:                                                  // *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
1864:                                                  // *
1865:                                                  this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1866:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1867:                                                  IE = ITAU;
1868:                                                  ITAUQ = IE + N;
1869:                                                  ITAUP = ITAUQ + N;
1870:                                                  IWORK = ITAUP + N;
1871:                                                  // *
1872:                                                  // *                    Bidiagonalize R in WORK(IR)
1873:                                                  // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
1874:                                                  // *
1875:                                                  this._dgebrd.Run(N, N, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
1876:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1877:                                                  // *
1878:                                                  // *                    Generate left bidiagonalizing vectors in WORK(IR)
1879:                                                  // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
1880:                                                  // *
1881:                                                  this._dorgbr.Run("Q", N, N, N, ref WORK, IR + o_work, LDWRKR
1882:                                                                   , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1883:                                                  IWORK = IE + N;
1884:                                                  // *
1885:                                                  // *                    Perform bidiagonal QR iteration, computing left
1886:                                                  // *                    singular vectors of R in WORK(IR)
1887:                                                  // *                    (Workspace: need N*N+BDSPAC)
1888:                                                  // *
1889:                                                  this._dbdsqr.Run("U", N, 0, N, 0, ref S, offset_s
1890:                                                                   , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
1891:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
1892:                                                  // *
1893:                                                  // *                    Multiply Q in U by left singular vectors of R in
1894:                                                  // *                    WORK(IR), storing result in A
1895:                                                  // *                    (Workspace: need N*N)
1896:                                                  // *
1897:                                                  this._dgemm.Run("N", "N", M, N, N, ONE
1898:                                                                  , U, offset_u, LDU, WORK, IR + o_work, LDWRKR, ZERO, ref A, offset_a
1899:                                                                  , LDA);
1900:                                                  // *
1901:                                                  // *                    Copy left singular vectors of A from A to U
1902:                                                  // *
1903:                                                  this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
1904:                                                                   , LDU);
1905:                                                  // *
1906:                                              }
1907:                                              else
1908:                                              {
1909:                                                  // *
1910:                                                  // *                    Insufficient workspace for a fast algorithm
1911:                                                  // *
1912:                                                  ITAU = 1;
1913:                                                  IWORK = ITAU + N;
1914:                                                  // *
1915:                                                  // *                    Compute A=Q*R, copying result to U
1916:                                                  // *                    (Workspace: need 2*N, prefer N+N*NB)
1917:                                                  // *
1918:                                                  this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
1919:                                                                   , LWORK - IWORK + 1, ref IERR);
1920:                                                  this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
1921:                                                                   , LDU);
1922:                                                  // *
1923:                                                  // *                    Generate Q in U
1924:                                                  // *                    (Workspace: need N+M, prefer N+M*NB)
1925:                                                  // *
1926:                                                  this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
1927:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1928:                                                  IE = ITAU;
1929:                                                  ITAUQ = IE + N;
1930:                                                  ITAUP = ITAUQ + N;
1931:                                                  IWORK = ITAUP + N;
1932:                                                  // *
1933:                                                  // *                    Zero out below R in A
1934:                                                  // *
1935:                                                  this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
1936:                                                                   , LDA);
1937:                                                  // *
1938:                                                  // *                    Bidiagonalize R in A
1939:                                                  // *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
1940:                                                  // *
1941:                                                  this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
1942:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
1943:                                                  // *
1944:                                                  // *                    Multiply Q in U by left bidiagonalizing vectors
1945:                                                  // *                    in A
1946:                                                  // *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
1947:                                                  // *
1948:                                                  this._dormbr.Run("Q", "R", "N", M, N, N
1949:                                                                   , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
1950:                                                                   , LWORK - IWORK + 1, ref IERR);
1951:                                                  IWORK = IE + N;
1952:                                                  // *
1953:                                                  // *                    Perform bidiagonal QR iteration, computing left
1954:                                                  // *                    singular vectors of A in U
1955:                                                  // *                    (Workspace: need BDSPAC)
1956:                                                  // *
1957:                                                  this._dbdsqr.Run("U", N, 0, M, 0, ref S, offset_s
1958:                                                                   , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref U, offset_u, LDU, ref DUM, offset_dum
1959:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
1960:                                                  // *
1961:                                              }
1962:                                              // *
1963:                                          }
1964:                                          else
1965:                                          {
1966:                                              if (WNTVO)
1967:                                              {
1968:                                                  // *
1969:                                                  // *                 Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1970:                                                  // *                 M left singular vectors to be computed in U and
1971:                                                  // *                 N right singular vectors to be overwritten on A
1972:                                                  // *
1973:                                                  if (LWORK >= 2 * N * N + Math.Max(N + M, Math.Max(4 * N, BDSPAC)))
1974:                                                  {
1975:                                                      // *
1976:                                                      // *                    Sufficient workspace for a fast algorithm
1977:                                                      // *
1978:                                                      IU = 1;
1979:                                                      if (LWORK >= WRKBL + 2 * LDA * N)
1980:                                                      {
1981:                                                          // *
1982:                                                          // *                       WORK(IU) is LDA by N and WORK(IR) is LDA by N
1983:                                                          // *
1984:                                                          LDWRKU = LDA;
1985:                                                          IR = IU + LDWRKU * N;
1986:                                                          LDWRKR = LDA;
1987:                                                      }
1988:                                                      else
1989:                                                      {
1990:                                                          if (LWORK >= WRKBL + (LDA + N) * N)
1991:                                                          {
1992:                                                              // *
1993:                                                              // *                       WORK(IU) is LDA by N and WORK(IR) is N by N
1994:                                                              // *
1995:                                                              LDWRKU = LDA;
1996:                                                              IR = IU + LDWRKU * N;
1997:                                                              LDWRKR = N;
1998:                                                          }
1999:                                                          else
2000:                                                          {
2001:                                                              // *
2002:                                                              // *                       WORK(IU) is N by N and WORK(IR) is N by N
2003:                                                              // *
2004:                                                              LDWRKU = N;
2005:                                                              IR = IU + LDWRKU * N;
2006:                                                              LDWRKR = N;
2007:                                                          }
2008:                                                      }
2009:                                                      ITAU = IR + LDWRKR * N;
2010:                                                      IWORK = ITAU + N;
2011:                                                      // *
2012:                                                      // *                    Compute A=Q*R, copying result to U
2013:                                                      // *                    (Workspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
2014:                                                      // *
2015:                                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2016:                                                                       , LWORK - IWORK + 1, ref IERR);
2017:                                                      this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2018:                                                                       , LDU);
2019:                                                      // *
2020:                                                      // *                    Generate Q in U
2021:                                                      // *                    (Workspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
2022:                                                      // *
2023:                                                      this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2024:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2025:                                                      // *
2026:                                                      // *                    Copy R to WORK(IU), zeroing out below it
2027:                                                      // *
2028:                                                      this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
2029:                                                                       , LDWRKU);
2030:                                                      this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
2031:                                                                       , LDWRKU);
2032:                                                      IE = ITAU;
2033:                                                      ITAUQ = IE + N;
2034:                                                      ITAUP = ITAUQ + N;
2035:                                                      IWORK = ITAUP + N;
2036:                                                      // *
2037:                                                      // *                    Bidiagonalize R in WORK(IU), copying result to
2038:                                                      // *                    WORK(IR)
2039:                                                      // *                    (Workspace: need 2*N*N+4*N,
2040:                                                      // *                                prefer 2*N*N+3*N+2*N*NB)
2041:                                                      // *
2042:                                                      this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
2043:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2044:                                                      this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
2045:                                                                       , LDWRKR);
2046:                                                      // *
2047:                                                      // *                    Generate left bidiagonalizing vectors in WORK(IU)
2048:                                                      // *                    (Workspace: need 2*N*N+4*N, prefer 2*N*N+3*N+N*NB)
2049:                                                      // *
2050:                                                      this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
2051:                                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2052:                                                      // *
2053:                                                      // *                    Generate right bidiagonalizing vectors in WORK(IR)
2054:                                                      // *                    (Workspace: need 2*N*N+4*N-1,
2055:                                                      // *                                prefer 2*N*N+3*N+(N-1)*NB)
2056:                                                      // *
2057:                                                      this._dorgbr.Run("P", N, N, N, ref WORK, IR + o_work, LDWRKR
2058:                                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2059:                                                      IWORK = IE + N;
2060:                                                      // *
2061:                                                      // *                    Perform bidiagonal QR iteration, computing left
2062:                                                      // *                    singular vectors of R in WORK(IU) and computing
2063:                                                      // *                    right singular vectors of R in WORK(IR)
2064:                                                      // *                    (Workspace: need 2*N*N+BDSPAC)
2065:                                                      // *
2066:                                                      this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
2067:                                                                       , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
2068:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
2069:                                                      // *
2070:                                                      // *                    Multiply Q in U by left singular vectors of R in
2071:                                                      // *                    WORK(IU), storing result in A
2072:                                                      // *                    (Workspace: need N*N)
2073:                                                      // *
2074:                                                      this._dgemm.Run("N", "N", M, N, N, ONE
2075:                                                                      , U, offset_u, LDU, WORK, IU + o_work, LDWRKU, ZERO, ref A, offset_a
2076:                                                                      , LDA);
2077:                                                      // *
2078:                                                      // *                    Copy left singular vectors of A from A to U
2079:                                                      // *
2080:                                                      this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
2081:                                                                       , LDU);
2082:                                                      // *
2083:                                                      // *                    Copy right singular vectors of R from WORK(IR) to A
2084:                                                      // *
2085:                                                      this._dlacpy.Run("F", N, N, WORK, IR + o_work, LDWRKR, ref A, offset_a
2086:                                                                       , LDA);
2087:                                                      // *
2088:                                                  }
2089:                                                  else
2090:                                                  {
2091:                                                      // *
2092:                                                      // *                    Insufficient workspace for a fast algorithm
2093:                                                      // *
2094:                                                      ITAU = 1;
2095:                                                      IWORK = ITAU + N;
2096:                                                      // *
2097:                                                      // *                    Compute A=Q*R, copying result to U
2098:                                                      // *                    (Workspace: need 2*N, prefer N+N*NB)
2099:                                                      // *
2100:                                                      this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2101:                                                                       , LWORK - IWORK + 1, ref IERR);
2102:                                                      this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2103:                                                                       , LDU);
2104:                                                      // *
2105:                                                      // *                    Generate Q in U
2106:                                                      // *                    (Workspace: need N+M, prefer N+M*NB)
2107:                                                      // *
2108:                                                      this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2109:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2110:                                                      IE = ITAU;
2111:                                                      ITAUQ = IE + N;
2112:                                                      ITAUP = ITAUQ + N;
2113:                                                      IWORK = ITAUP + N;
2114:                                                      // *
2115:                                                      // *                    Zero out below R in A
2116:                                                      // *
2117:                                                      this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref A, 2+1 * LDA + o_a
2118:                                                                       , LDA);
2119:                                                      // *
2120:                                                      // *                    Bidiagonalize R in A
2121:                                                      // *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
2122:                                                      // *
2123:                                                      this._dgebrd.Run(N, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2124:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2125:                                                      // *
2126:                                                      // *                    Multiply Q in U by left bidiagonalizing vectors
2127:                                                      // *                    in A
2128:                                                      // *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
2129:                                                      // *
2130:                                                      this._dormbr.Run("Q", "R", "N", M, N, N
2131:                                                                       , ref A, offset_a, LDA, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
2132:                                                                       , LWORK - IWORK + 1, ref IERR);
2133:                                                      // *
2134:                                                      // *                    Generate right bidiagonalizing vectors in A
2135:                                                      // *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2136:                                                      // *
2137:                                                      this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
2138:                                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2139:                                                      IWORK = IE + N;
2140:                                                      // *
2141:                                                      // *                    Perform bidiagonal QR iteration, computing left
2142:                                                      // *                    singular vectors of A in U and computing right
2143:                                                      // *                    singular vectors of A in A
2144:                                                      // *                    (Workspace: need BDSPAC)
2145:                                                      // *
2146:                                                      this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
2147:                                                                       , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
2148:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
2149:                                                      // *
2150:                                                  }
2151:                                                  // *
2152:                                              }
2153:                                              else
2154:                                              {
2155:                                                  if (WNTVAS)
2156:                                                  {
2157:                                                      // *
2158:                                                      // *                 Path 9 (M much larger than N, JOBU='A', JOBVT='S'
2159:                                                      // *                         or 'A')
2160:                                                      // *                 M left singular vectors to be computed in U and
2161:                                                      // *                 N right singular vectors to be computed in VT
2162:                                                      // *
2163:                                                      if (LWORK >= N * N + Math.Max(N + M, Math.Max(4 * N, BDSPAC)))
2164:                                                      {
2165:                                                          // *
2166:                                                          // *                    Sufficient workspace for a fast algorithm
2167:                                                          // *
2168:                                                          IU = 1;
2169:                                                          if (LWORK >= WRKBL + LDA * N)
2170:                                                          {
2171:                                                              // *
2172:                                                              // *                       WORK(IU) is LDA by N
2173:                                                              // *
2174:                                                              LDWRKU = LDA;
2175:                                                          }
2176:                                                          else
2177:                                                          {
2178:                                                              // *
2179:                                                              // *                       WORK(IU) is N by N
2180:                                                              // *
2181:                                                              LDWRKU = N;
2182:                                                          }
2183:                                                          ITAU = IU + LDWRKU * N;
2184:                                                          IWORK = ITAU + N;
2185:                                                          // *
2186:                                                          // *                    Compute A=Q*R, copying result to U
2187:                                                          // *                    (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
2188:                                                          // *
2189:                                                          this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2190:                                                                           , LWORK - IWORK + 1, ref IERR);
2191:                                                          this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2192:                                                                           , LDU);
2193:                                                          // *
2194:                                                          // *                    Generate Q in U
2195:                                                          // *                    (Workspace: need N*N+N+M, prefer N*N+N+M*NB)
2196:                                                          // *
2197:                                                          this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2198:                                                                           , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2199:                                                          // *
2200:                                                          // *                    Copy R to WORK(IU), zeroing out below it
2201:                                                          // *
2202:                                                          this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref WORK, IU + o_work
2203:                                                                           , LDWRKU);
2204:                                                          this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref WORK, IU + 1 + o_work
2205:                                                                           , LDWRKU);
2206:                                                          IE = ITAU;
2207:                                                          ITAUQ = IE + N;
2208:                                                          ITAUP = ITAUQ + N;
2209:                                                          IWORK = ITAUP + N;
2210:                                                          // *
2211:                                                          // *                    Bidiagonalize R in WORK(IU), copying result to VT
2212:                                                          // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
2213:                                                          // *
2214:                                                          this._dgebrd.Run(N, N, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
2215:                                                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2216:                                                          this._dlacpy.Run("U", N, N, WORK, IU + o_work, LDWRKU, ref VT, offset_vt
2217:                                                                           , LDVT);
2218:                                                          // *
2219:                                                          // *                    Generate left bidiagonalizing vectors in WORK(IU)
2220:                                                          // *                    (Workspace: need N*N+4*N, prefer N*N+3*N+N*NB)
2221:                                                          // *
2222:                                                          this._dorgbr.Run("Q", N, N, N, ref WORK, IU + o_work, LDWRKU
2223:                                                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2224:                                                          // *
2225:                                                          // *                    Generate right bidiagonalizing vectors in VT
2226:                                                          // *                    (Workspace: need N*N+4*N-1,
2227:                                                          // *                                prefer N*N+3*N+(N-1)*NB)
2228:                                                          // *
2229:                                                          this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
2230:                                                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2231:                                                          IWORK = IE + N;
2232:                                                          // *
2233:                                                          // *                    Perform bidiagonal QR iteration, computing left
2234:                                                          // *                    singular vectors of R in WORK(IU) and computing
2235:                                                          // *                    right singular vectors of R in VT
2236:                                                          // *                    (Workspace: need N*N+BDSPAC)
2237:                                                          // *
2238:                                                          this._dbdsqr.Run("U", N, N, N, 0, ref S, offset_s
2239:                                                                           , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref WORK, IU + o_work, LDWRKU, ref DUM, offset_dum
2240:                                                                           , 1, ref WORK, IWORK + o_work, ref INFO);
2241:                                                          // *
2242:                                                          // *                    Multiply Q in U by left singular vectors of R in
2243:                                                          // *                    WORK(IU), storing result in A
2244:                                                          // *                    (Workspace: need N*N)
2245:                                                          // *
2246:                                                          this._dgemm.Run("N", "N", M, N, N, ONE
2247:                                                                          , U, offset_u, LDU, WORK, IU + o_work, LDWRKU, ZERO, ref A, offset_a
2248:                                                                          , LDA);
2249:                                                          // *
2250:                                                          // *                    Copy left singular vectors of A from A to U
2251:                                                          // *
2252:                                                          this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref U, offset_u
2253:                                                                           , LDU);
2254:                                                          // *
2255:                                                      }
2256:                                                      else
2257:                                                      {
2258:                                                          // *
2259:                                                          // *                    Insufficient workspace for a fast algorithm
2260:                                                          // *
2261:                                                          ITAU = 1;
2262:                                                          IWORK = ITAU + N;
2263:                                                          // *
2264:                                                          // *                    Compute A=Q*R, copying result to U
2265:                                                          // *                    (Workspace: need 2*N, prefer N+N*NB)
2266:                                                          // *
2267:                                                          this._dgeqrf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2268:                                                                           , LWORK - IWORK + 1, ref IERR);
2269:                                                          this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2270:                                                                           , LDU);
2271:                                                          // *
2272:                                                          // *                    Generate Q in U
2273:                                                          // *                    (Workspace: need N+M, prefer N+M*NB)
2274:                                                          // *
2275:                                                          this._dorgqr.Run(M, M, N, ref U, offset_u, LDU, WORK, ITAU + o_work
2276:                                                                           , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2277:                                                          // *
2278:                                                          // *                    Copy R from A to VT, zeroing out below it
2279:                                                          // *
2280:                                                          this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
2281:                                                                           , LDVT);
2282:                                                          if (N > 1)
2283:                                                          {
2284:                                                              this._dlaset.Run("L", N - 1, N - 1, ZERO, ZERO, ref VT, 2+1 * LDVT + o_vt
2285:                                                                               , LDVT);
2286:                                                          }
2287:                                                          IE = ITAU;
2288:                                                          ITAUQ = IE + N;
2289:                                                          ITAUP = ITAUQ + N;
2290:                                                          IWORK = ITAUP + N;
2291:                                                          // *
2292:                                                          // *                    Bidiagonalize R in VT
2293:                                                          // *                    (Workspace: need 4*N, prefer 3*N+2*N*NB)
2294:                                                          // *
2295:                                                          this._dgebrd.Run(N, N, ref VT, offset_vt, LDVT, ref S, offset_s, ref WORK, IE + o_work
2296:                                                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2297:                                                          // *
2298:                                                          // *                    Multiply Q in U by left bidiagonalizing vectors
2299:                                                          // *                    in VT
2300:                                                          // *                    (Workspace: need 3*N+M, prefer 3*N+M*NB)
2301:                                                          // *
2302:                                                          this._dormbr.Run("Q", "R", "N", M, N, N
2303:                                                                           , ref VT, offset_vt, LDVT, WORK, ITAUQ + o_work, ref U, offset_u, LDU, ref WORK, IWORK + o_work
2304:                                                                           , LWORK - IWORK + 1, ref IERR);
2305:                                                          // *
2306:                                                          // *                    Generate right bidiagonalizing vectors in VT
2307:                                                          // *                    (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2308:                                                          // *
2309:                                                          this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
2310:                                                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2311:                                                          IWORK = IE + N;
2312:                                                          // *
2313:                                                          // *                    Perform bidiagonal QR iteration, computing left
2314:                                                          // *                    singular vectors of A in U and computing right
2315:                                                          // *                    singular vectors of A in VT
2316:                                                          // *                    (Workspace: need BDSPAC)
2317:                                                          // *
2318:                                                          this._dbdsqr.Run("U", N, N, M, 0, ref S, offset_s
2319:                                                                           , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
2320:                                                                           , 1, ref WORK, IWORK + o_work, ref INFO);
2321:                                                          // *
2322:                                                      }
2323:                                                      // *
2324:                                                  }
2325:                                              }
2326:                                          }
2327:                                          // *
2328:                                      }
2329:                                  }
2330:                              }
2331:                          }
2332:                      }
2333:                      // *
2334:                  }
2335:                  else
2336:                  {
2337:                      // *
2338:                      // *           M .LT. MNTHR
2339:                      // *
2340:                      // *           Path 10 (M at least N, but not much larger)
2341:                      // *           Reduce to bidiagonal form without QR decomposition
2342:                      // *
2343:                      IE = 1;
2344:                      ITAUQ = IE + N;
2345:                      ITAUP = ITAUQ + N;
2346:                      IWORK = ITAUP + N;
2347:                      // *
2348:                      // *           Bidiagonalize A
2349:                      // *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
2350:                      // *
2351:                      this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2352:                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2353:                      if (WNTUAS)
2354:                      {
2355:                          // *
2356:                          // *              If left singular vectors desired in U, copy result to U
2357:                          // *              and generate left bidiagonalizing vectors in U
2358:                          // *              (Workspace: need 3*N+NCU, prefer 3*N+NCU*NB)
2359:                          // *
2360:                          this._dlacpy.Run("L", M, N, A, offset_a, LDA, ref U, offset_u
2361:                                           , LDU);
2362:                          if (WNTUS) NCU = N;
2363:                          if (WNTUA) NCU = M;
2364:                          this._dorgbr.Run("Q", M, NCU, N, ref U, offset_u, LDU
2365:                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2366:                      }
2367:                      if (WNTVAS)
2368:                      {
2369:                          // *
2370:                          // *              If right singular vectors desired in VT, copy result to
2371:                          // *              VT and generate right bidiagonalizing vectors in VT
2372:                          // *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2373:                          // *
2374:                          this._dlacpy.Run("U", N, N, A, offset_a, LDA, ref VT, offset_vt
2375:                                           , LDVT);
2376:                          this._dorgbr.Run("P", N, N, N, ref VT, offset_vt, LDVT
2377:                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2378:                      }
2379:                      if (WNTUO)
2380:                      {
2381:                          // *
2382:                          // *              If left singular vectors desired in A, generate left
2383:                          // *              bidiagonalizing vectors in A
2384:                          // *              (Workspace: need 4*N, prefer 3*N+N*NB)
2385:                          // *
2386:                          this._dorgbr.Run("Q", M, N, N, ref A, offset_a, LDA
2387:                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2388:                      }
2389:                      if (WNTVO)
2390:                      {
2391:                          // *
2392:                          // *              If right singular vectors desired in A, generate right
2393:                          // *              bidiagonalizing vectors in A
2394:                          // *              (Workspace: need 4*N-1, prefer 3*N+(N-1)*NB)
2395:                          // *
2396:                          this._dorgbr.Run("P", N, N, N, ref A, offset_a, LDA
2397:                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2398:                      }
2399:                      IWORK = IE + N;
2400:                      if (WNTUAS || WNTUO) NRU = M;
2401:                      if (WNTUN) NRU = 0;
2402:                      if (WNTVAS || WNTVO) NCVT = N;
2403:                      if (WNTVN) NCVT = 0;
2404:                      if ((!WNTUO) && (!WNTVO))
2405:                      {
2406:                          // *
2407:                          // *              Perform bidiagonal QR iteration, if desired, computing
2408:                          // *              left singular vectors in U and computing right singular
2409:                          // *              vectors in VT
2410:                          // *              (Workspace: need BDSPAC)
2411:                          // *
2412:                          this._dbdsqr.Run("U", N, NCVT, NRU, 0, ref S, offset_s
2413:                                           , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
2414:                                           , 1, ref WORK, IWORK + o_work, ref INFO);
2415:                      }
2416:                      else
2417:                      {
2418:                          if ((!WNTUO) && WNTVO)
2419:                          {
2420:                              // *
2421:                              // *              Perform bidiagonal QR iteration, if desired, computing
2422:                              // *              left singular vectors in U and computing right singular
2423:                              // *              vectors in A
2424:                              // *              (Workspace: need BDSPAC)
2425:                              // *
2426:                              this._dbdsqr.Run("U", N, NCVT, NRU, 0, ref S, offset_s
2427:                                               , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
2428:                                               , 1, ref WORK, IWORK + o_work, ref INFO);
2429:                          }
2430:                          else
2431:                          {
2432:                              // *
2433:                              // *              Perform bidiagonal QR iteration, if desired, computing
2434:                              // *              left singular vectors in A and computing right singular
2435:                              // *              vectors in VT
2436:                              // *              (Workspace: need BDSPAC)
2437:                              // *
2438:                              this._dbdsqr.Run("U", N, NCVT, NRU, 0, ref S, offset_s
2439:                                               , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
2440:                                               , 1, ref WORK, IWORK + o_work, ref INFO);
2441:                          }
2442:                      }
2443:                      // *
2444:                  }
2445:                  // *
2446:              }
2447:              else
2448:              {
2449:                  // *
2450:                  // *        A has more columns than rows. If A has sufficiently more
2451:                  // *        columns than rows, first reduce using the LQ decomposition (if
2452:                  // *        sufficient workspace available)
2453:                  // *
2454:                  if (N >= MNTHR)
2455:                  {
2456:                      // *
2457:                      if (WNTVN)
2458:                      {
2459:                          // *
2460:                          // *              Path 1t(N much larger than M, JOBVT='N')
2461:                          // *              No right singular vectors to be computed
2462:                          // *
2463:                          ITAU = 1;
2464:                          IWORK = ITAU + M;
2465:                          // *
2466:                          // *              Compute A=L*Q
2467:                          // *              (Workspace: need 2*M, prefer M+M*NB)
2468:                          // *
2469:                          this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2470:                                           , LWORK - IWORK + 1, ref IERR);
2471:                          // *
2472:                          // *              Zero out above L
2473:                          // *
2474:                          this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
2475:                                           , LDA);
2476:                          IE = 1;
2477:                          ITAUQ = IE + M;
2478:                          ITAUP = ITAUQ + M;
2479:                          IWORK = ITAUP + M;
2480:                          // *
2481:                          // *              Bidiagonalize L in A
2482:                          // *              (Workspace: need 4*M, prefer 3*M+2*M*NB)
2483:                          // *
2484:                          this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2485:                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2486:                          if (WNTUO || WNTUAS)
2487:                          {
2488:                              // *
2489:                              // *                 If left singular vectors desired, generate Q
2490:                              // *                 (Workspace: need 4*M, prefer 3*M+M*NB)
2491:                              // *
2492:                              this._dorgbr.Run("Q", M, M, M, ref A, offset_a, LDA
2493:                                               , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2494:                          }
2495:                          IWORK = IE + M;
2496:                          NRU = 0;
2497:                          if (WNTUO || WNTUAS) NRU = M;
2498:                          // *
2499:                          // *              Perform bidiagonal QR iteration, computing left singular
2500:                          // *              vectors of A in A if desired
2501:                          // *              (Workspace: need BDSPAC)
2502:                          // *
2503:                          this._dbdsqr.Run("U", M, 0, NRU, 0, ref S, offset_s
2504:                                           , ref WORK, IE + o_work, ref DUM, offset_dum, 1, ref A, offset_a, LDA, ref DUM, offset_dum
2505:                                           , 1, ref WORK, IWORK + o_work, ref INFO);
2506:                          // *
2507:                          // *              If left singular vectors desired in U, copy them there
2508:                          // *
2509:                          if (WNTUAS)
2510:                          {
2511:                              this._dlacpy.Run("F", M, M, A, offset_a, LDA, ref U, offset_u
2512:                                               , LDU);
2513:                          }
2514:                          // *
2515:                      }
2516:                      else
2517:                      {
2518:                          if (WNTVO && WNTUN)
2519:                          {
2520:                              // *
2521:                              // *              Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2522:                              // *              M right singular vectors to be overwritten on A and
2523:                              // *              no left singular vectors to be computed
2524:                              // *
2525:                              if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
2526:                              {
2527:                                  // *
2528:                                  // *                 Sufficient workspace for a fast algorithm
2529:                                  // *
2530:                                  IR = 1;
2531:                                  if (LWORK >= Math.Max(WRKBL, LDA * N + M) + LDA * M)
2532:                                  {
2533:                                      // *
2534:                                      // *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
2535:                                      // *
2536:                                      LDWRKU = LDA;
2537:                                      CHUNK = N;
2538:                                      LDWRKR = LDA;
2539:                                  }
2540:                                  else
2541:                                  {
2542:                                      if (LWORK >= Math.Max(WRKBL, LDA * N + M) + M * M)
2543:                                      {
2544:                                          // *
2545:                                          // *                    WORK(IU) is LDA by N and WORK(IR) is M by M
2546:                                          // *
2547:                                          LDWRKU = LDA;
2548:                                          CHUNK = N;
2549:                                          LDWRKR = M;
2550:                                      }
2551:                                      else
2552:                                      {
2553:                                          // *
2554:                                          // *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
2555:                                          // *
2556:                                          LDWRKU = M;
2557:                                          CHUNK = (LWORK - M * M - M) / M;
2558:                                          LDWRKR = M;
2559:                                      }
2560:                                  }
2561:                                  ITAU = IR + LDWRKR * M;
2562:                                  IWORK = ITAU + M;
2563:                                  // *
2564:                                  // *                 Compute A=L*Q
2565:                                  // *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2566:                                  // *
2567:                                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2568:                                                   , LWORK - IWORK + 1, ref IERR);
2569:                                  // *
2570:                                  // *                 Copy L to WORK(IR) and zero out above it
2571:                                  // *
2572:                                  this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IR + o_work
2573:                                                   , LDWRKR);
2574:                                  this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IR + LDWRKR + o_work
2575:                                                   , LDWRKR);
2576:                                  // *
2577:                                  // *                 Generate Q in A
2578:                                  // *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2579:                                  // *
2580:                                  this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2581:                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2582:                                  IE = ITAU;
2583:                                  ITAUQ = IE + M;
2584:                                  ITAUP = ITAUQ + M;
2585:                                  IWORK = ITAUP + M;
2586:                                  // *
2587:                                  // *                 Bidiagonalize L in WORK(IR)
2588:                                  // *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2589:                                  // *
2590:                                  this._dgebrd.Run(M, M, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
2591:                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2592:                                  // *
2593:                                  // *                 Generate right vectors bidiagonalizing L
2594:                                  // *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
2595:                                  // *
2596:                                  this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
2597:                                                   , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2598:                                  IWORK = IE + M;
2599:                                  // *
2600:                                  // *                 Perform bidiagonal QR iteration, computing right
2601:                                  // *                 singular vectors of L in WORK(IR)
2602:                                  // *                 (Workspace: need M*M+BDSPAC)
2603:                                  // *
2604:                                  this._dbdsqr.Run("U", M, M, 0, 0, ref S, offset_s
2605:                                                   , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum, 1, ref DUM, offset_dum
2606:                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
2607:                                  IU = IE + M;
2608:                                  // *
2609:                                  // *                 Multiply right singular vectors of L in WORK(IR) by Q
2610:                                  // *                 in A, storing result in WORK(IU) and copying to A
2611:                                  // *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M)
2612:                                  // *
2613:                                  for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)
2614:                                  {
2615:                                      BLK = Math.Min(N - I + 1, CHUNK);
2616:                                      this._dgemm.Run("N", "N", M, BLK, M, ONE
2617:                                                      , WORK, IR + o_work, LDWRKR, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IU + o_work
2618:                                                      , LDWRKU);
2619:                                      this._dlacpy.Run("F", M, BLK, WORK, IU + o_work, LDWRKU, ref A, 1+I * LDA + o_a
2620:                                                       , LDA);
2621:                                  }
2622:                                  // *
2623:                              }
2624:                              else
2625:                              {
2626:                                  // *
2627:                                  // *                 Insufficient workspace for a fast algorithm
2628:                                  // *
2629:                                  IE = 1;
2630:                                  ITAUQ = IE + M;
2631:                                  ITAUP = ITAUQ + M;
2632:                                  IWORK = ITAUP + M;
2633:                                  // *
2634:                                  // *                 Bidiagonalize A
2635:                                  // *                 (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
2636:                                  // *
2637:                                  this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2638:                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2639:                                  // *
2640:                                  // *                 Generate right vectors bidiagonalizing A
2641:                                  // *                 (Workspace: need 4*M, prefer 3*M+M*NB)
2642:                                  // *
2643:                                  this._dorgbr.Run("P", M, N, M, ref A, offset_a, LDA
2644:                                                   , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2645:                                  IWORK = IE + M;
2646:                                  // *
2647:                                  // *                 Perform bidiagonal QR iteration, computing right
2648:                                  // *                 singular vectors of A in A
2649:                                  // *                 (Workspace: need BDSPAC)
2650:                                  // *
2651:                                  this._dbdsqr.Run("L", M, N, 0, 0, ref S, offset_s
2652:                                                   , ref WORK, IE + o_work, ref A, offset_a, LDA, ref DUM, offset_dum, 1, ref DUM, offset_dum
2653:                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
2654:                                  // *
2655:                              }
2656:                              // *
2657:                          }
2658:                          else
2659:                          {
2660:                              if (WNTVO && WNTUAS)
2661:                              {
2662:                                  // *
2663:                                  // *              Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2664:                                  // *              M right singular vectors to be overwritten on A and
2665:                                  // *              M left singular vectors to be computed in U
2666:                                  // *
2667:                                  if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
2668:                                  {
2669:                                      // *
2670:                                      // *                 Sufficient workspace for a fast algorithm
2671:                                      // *
2672:                                      IR = 1;
2673:                                      if (LWORK >= Math.Max(WRKBL, LDA * N + M) + LDA * M)
2674:                                      {
2675:                                          // *
2676:                                          // *                    WORK(IU) is LDA by N and WORK(IR) is LDA by M
2677:                                          // *
2678:                                          LDWRKU = LDA;
2679:                                          CHUNK = N;
2680:                                          LDWRKR = LDA;
2681:                                      }
2682:                                      else
2683:                                      {
2684:                                          if (LWORK >= Math.Max(WRKBL, LDA * N + M) + M * M)
2685:                                          {
2686:                                              // *
2687:                                              // *                    WORK(IU) is LDA by N and WORK(IR) is M by M
2688:                                              // *
2689:                                              LDWRKU = LDA;
2690:                                              CHUNK = N;
2691:                                              LDWRKR = M;
2692:                                          }
2693:                                          else
2694:                                          {
2695:                                              // *
2696:                                              // *                    WORK(IU) is M by CHUNK and WORK(IR) is M by M
2697:                                              // *
2698:                                              LDWRKU = M;
2699:                                              CHUNK = (LWORK - M * M - M) / M;
2700:                                              LDWRKR = M;
2701:                                          }
2702:                                      }
2703:                                      ITAU = IR + LDWRKR * M;
2704:                                      IWORK = ITAU + M;
2705:                                      // *
2706:                                      // *                 Compute A=L*Q
2707:                                      // *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2708:                                      // *
2709:                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2710:                                                       , LWORK - IWORK + 1, ref IERR);
2711:                                      // *
2712:                                      // *                 Copy L to U, zeroing about above it
2713:                                      // *
2714:                                      this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
2715:                                                       , LDU);
2716:                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
2717:                                                       , LDU);
2718:                                      // *
2719:                                      // *                 Generate Q in A
2720:                                      // *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2721:                                      // *
2722:                                      this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2723:                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2724:                                      IE = ITAU;
2725:                                      ITAUQ = IE + M;
2726:                                      ITAUP = ITAUQ + M;
2727:                                      IWORK = ITAUP + M;
2728:                                      // *
2729:                                      // *                 Bidiagonalize L in U, copying result to WORK(IR)
2730:                                      // *                 (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2731:                                      // *
2732:                                      this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
2733:                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2734:                                      this._dlacpy.Run("U", M, M, U, offset_u, LDU, ref WORK, IR + o_work
2735:                                                       , LDWRKR);
2736:                                      // *
2737:                                      // *                 Generate right vectors bidiagonalizing L in WORK(IR)
2738:                                      // *                 (Workspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
2739:                                      // *
2740:                                      this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
2741:                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2742:                                      // *
2743:                                      // *                 Generate left vectors bidiagonalizing L in U
2744:                                      // *                 (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
2745:                                      // *
2746:                                      this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
2747:                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2748:                                      IWORK = IE + M;
2749:                                      // *
2750:                                      // *                 Perform bidiagonal QR iteration, computing left
2751:                                      // *                 singular vectors of L in U, and computing right
2752:                                      // *                 singular vectors of L in WORK(IR)
2753:                                      // *                 (Workspace: need M*M+BDSPAC)
2754:                                      // *
2755:                                      this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
2756:                                                       , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref U, offset_u, LDU, ref DUM, offset_dum
2757:                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
2758:                                      IU = IE + M;
2759:                                      // *
2760:                                      // *                 Multiply right singular vectors of L in WORK(IR) by Q
2761:                                      // *                 in A, storing result in WORK(IU) and copying to A
2762:                                      // *                 (Workspace: need M*M+2*M, prefer M*M+M*N+M))
2763:                                      // *
2764:                                      for (I = 1; (CHUNK >= 0) ? (I <= N) : (I >= N); I += CHUNK)
2765:                                      {
2766:                                          BLK = Math.Min(N - I + 1, CHUNK);
2767:                                          this._dgemm.Run("N", "N", M, BLK, M, ONE
2768:                                                          , WORK, IR + o_work, LDWRKR, A, 1+I * LDA + o_a, LDA, ZERO, ref WORK, IU + o_work
2769:                                                          , LDWRKU);
2770:                                          this._dlacpy.Run("F", M, BLK, WORK, IU + o_work, LDWRKU, ref A, 1+I * LDA + o_a
2771:                                                           , LDA);
2772:                                      }
2773:                                      // *
2774:                                  }
2775:                                  else
2776:                                  {
2777:                                      // *
2778:                                      // *                 Insufficient workspace for a fast algorithm
2779:                                      // *
2780:                                      ITAU = 1;
2781:                                      IWORK = ITAU + M;
2782:                                      // *
2783:                                      // *                 Compute A=L*Q
2784:                                      // *                 (Workspace: need 2*M, prefer M+M*NB)
2785:                                      // *
2786:                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2787:                                                       , LWORK - IWORK + 1, ref IERR);
2788:                                      // *
2789:                                      // *                 Copy L to U, zeroing out above it
2790:                                      // *
2791:                                      this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
2792:                                                       , LDU);
2793:                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
2794:                                                       , LDU);
2795:                                      // *
2796:                                      // *                 Generate Q in A
2797:                                      // *                 (Workspace: need 2*M, prefer M+M*NB)
2798:                                      // *
2799:                                      this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2800:                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2801:                                      IE = ITAU;
2802:                                      ITAUQ = IE + M;
2803:                                      ITAUP = ITAUQ + M;
2804:                                      IWORK = ITAUP + M;
2805:                                      // *
2806:                                      // *                 Bidiagonalize L in U
2807:                                      // *                 (Workspace: need 4*M, prefer 3*M+2*M*NB)
2808:                                      // *
2809:                                      this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
2810:                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2811:                                      // *
2812:                                      // *                 Multiply right vectors bidiagonalizing L by Q in A
2813:                                      // *                 (Workspace: need 3*M+N, prefer 3*M+N*NB)
2814:                                      // *
2815:                                      this._dormbr.Run("P", "L", "T", M, N, M
2816:                                                       , ref U, offset_u, LDU, WORK, ITAUP + o_work, ref A, offset_a, LDA, ref WORK, IWORK + o_work
2817:                                                       , LWORK - IWORK + 1, ref IERR);
2818:                                      // *
2819:                                      // *                 Generate left vectors bidiagonalizing L in U
2820:                                      // *                 (Workspace: need 4*M, prefer 3*M+M*NB)
2821:                                      // *
2822:                                      this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
2823:                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2824:                                      IWORK = IE + M;
2825:                                      // *
2826:                                      // *                 Perform bidiagonal QR iteration, computing left
2827:                                      // *                 singular vectors of A in U and computing right
2828:                                      // *                 singular vectors of A in A
2829:                                      // *                 (Workspace: need BDSPAC)
2830:                                      // *
2831:                                      this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
2832:                                                       , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
2833:                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
2834:                                      // *
2835:                                  }
2836:                                  // *
2837:                              }
2838:                              else
2839:                              {
2840:                                  if (WNTVS)
2841:                                  {
2842:                                      // *
2843:                                      if (WNTUN)
2844:                                      {
2845:                                          // *
2846:                                          // *                 Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2847:                                          // *                 M right singular vectors to be computed in VT and
2848:                                          // *                 no left singular vectors to be computed
2849:                                          // *
2850:                                          if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
2851:                                          {
2852:                                              // *
2853:                                              // *                    Sufficient workspace for a fast algorithm
2854:                                              // *
2855:                                              IR = 1;
2856:                                              if (LWORK >= WRKBL + LDA * M)
2857:                                              {
2858:                                                  // *
2859:                                                  // *                       WORK(IR) is LDA by M
2860:                                                  // *
2861:                                                  LDWRKR = LDA;
2862:                                              }
2863:                                              else
2864:                                              {
2865:                                                  // *
2866:                                                  // *                       WORK(IR) is M by M
2867:                                                  // *
2868:                                                  LDWRKR = M;
2869:                                              }
2870:                                              ITAU = IR + LDWRKR * M;
2871:                                              IWORK = ITAU + M;
2872:                                              // *
2873:                                              // *                    Compute A=L*Q
2874:                                              // *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2875:                                              // *
2876:                                              this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2877:                                                               , LWORK - IWORK + 1, ref IERR);
2878:                                              // *
2879:                                              // *                    Copy L to WORK(IR), zeroing out above it
2880:                                              // *
2881:                                              this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IR + o_work
2882:                                                               , LDWRKR);
2883:                                              this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IR + LDWRKR + o_work
2884:                                                               , LDWRKR);
2885:                                              // *
2886:                                              // *                    Generate Q in A
2887:                                              // *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
2888:                                              // *
2889:                                              this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
2890:                                                               , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2891:                                              IE = ITAU;
2892:                                              ITAUQ = IE + M;
2893:                                              ITAUP = ITAUQ + M;
2894:                                              IWORK = ITAUP + M;
2895:                                              // *
2896:                                              // *                    Bidiagonalize L in WORK(IR)
2897:                                              // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
2898:                                              // *
2899:                                              this._dgebrd.Run(M, M, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
2900:                                                               , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2901:                                              // *
2902:                                              // *                    Generate right vectors bidiagonalizing L in
2903:                                              // *                    WORK(IR)
2904:                                              // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
2905:                                              // *
2906:                                              this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
2907:                                                               , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2908:                                              IWORK = IE + M;
2909:                                              // *
2910:                                              // *                    Perform bidiagonal QR iteration, computing right
2911:                                              // *                    singular vectors of L in WORK(IR)
2912:                                              // *                    (Workspace: need M*M+BDSPAC)
2913:                                              // *
2914:                                              this._dbdsqr.Run("U", M, M, 0, 0, ref S, offset_s
2915:                                                               , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum, 1, ref DUM, offset_dum
2916:                                                               , 1, ref WORK, IWORK + o_work, ref INFO);
2917:                                              // *
2918:                                              // *                    Multiply right singular vectors of L in WORK(IR) by
2919:                                              // *                    Q in A, storing result in VT
2920:                                              // *                    (Workspace: need M*M)
2921:                                              // *
2922:                                              this._dgemm.Run("N", "N", M, N, M, ONE
2923:                                                              , WORK, IR + o_work, LDWRKR, A, offset_a, LDA, ZERO, ref VT, offset_vt
2924:                                                              , LDVT);
2925:                                              // *
2926:                                          }
2927:                                          else
2928:                                          {
2929:                                              // *
2930:                                              // *                    Insufficient workspace for a fast algorithm
2931:                                              // *
2932:                                              ITAU = 1;
2933:                                              IWORK = ITAU + M;
2934:                                              // *
2935:                                              // *                    Compute A=L*Q
2936:                                              // *                    (Workspace: need 2*M, prefer M+M*NB)
2937:                                              // *
2938:                                              this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
2939:                                                               , LWORK - IWORK + 1, ref IERR);
2940:                                              // *
2941:                                              // *                    Copy result to VT
2942:                                              // *
2943:                                              this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
2944:                                                               , LDVT);
2945:                                              // *
2946:                                              // *                    Generate Q in VT
2947:                                              // *                    (Workspace: need 2*M, prefer M+M*NB)
2948:                                              // *
2949:                                              this._dorglq.Run(M, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
2950:                                                               , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2951:                                              IE = ITAU;
2952:                                              ITAUQ = IE + M;
2953:                                              ITAUP = ITAUQ + M;
2954:                                              IWORK = ITAUP + M;
2955:                                              // *
2956:                                              // *                    Zero out above L in A
2957:                                              // *
2958:                                              this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
2959:                                                               , LDA);
2960:                                              // *
2961:                                              // *                    Bidiagonalize L in A
2962:                                              // *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
2963:                                              // *
2964:                                              this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
2965:                                                               , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
2966:                                              // *
2967:                                              // *                    Multiply right vectors bidiagonalizing L by Q in VT
2968:                                              // *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
2969:                                              // *
2970:                                              this._dormbr.Run("P", "L", "T", M, N, M
2971:                                                               , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
2972:                                                               , LWORK - IWORK + 1, ref IERR);
2973:                                              IWORK = IE + M;
2974:                                              // *
2975:                                              // *                    Perform bidiagonal QR iteration, computing right
2976:                                              // *                    singular vectors of A in VT
2977:                                              // *                    (Workspace: need BDSPAC)
2978:                                              // *
2979:                                              this._dbdsqr.Run("U", M, N, 0, 0, ref S, offset_s
2980:                                                               , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref DUM, offset_dum, 1, ref DUM, offset_dum
2981:                                                               , 1, ref WORK, IWORK + o_work, ref INFO);
2982:                                              // *
2983:                                          }
2984:                                          // *
2985:                                      }
2986:                                      else
2987:                                      {
2988:                                          if (WNTUO)
2989:                                          {
2990:                                              // *
2991:                                              // *                 Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2992:                                              // *                 M right singular vectors to be computed in VT and
2993:                                              // *                 M left singular vectors to be overwritten on A
2994:                                              // *
2995:                                              if (LWORK >= 2 * M * M + Math.Max(4 * M, BDSPAC))
2996:                                              {
2997:                                                  // *
2998:                                                  // *                    Sufficient workspace for a fast algorithm
2999:                                                  // *
3000:                                                  IU = 1;
3001:                                                  if (LWORK >= WRKBL + 2 * LDA * M)
3002:                                                  {
3003:                                                      // *
3004:                                                      // *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
3005:                                                      // *
3006:                                                      LDWRKU = LDA;
3007:                                                      IR = IU + LDWRKU * M;
3008:                                                      LDWRKR = LDA;
3009:                                                  }
3010:                                                  else
3011:                                                  {
3012:                                                      if (LWORK >= WRKBL + (LDA + M) * M)
3013:                                                      {
3014:                                                          // *
3015:                                                          // *                       WORK(IU) is LDA by M and WORK(IR) is M by M
3016:                                                          // *
3017:                                                          LDWRKU = LDA;
3018:                                                          IR = IU + LDWRKU * M;
3019:                                                          LDWRKR = M;
3020:                                                      }
3021:                                                      else
3022:                                                      {
3023:                                                          // *
3024:                                                          // *                       WORK(IU) is M by M and WORK(IR) is M by M
3025:                                                          // *
3026:                                                          LDWRKU = M;
3027:                                                          IR = IU + LDWRKU * M;
3028:                                                          LDWRKR = M;
3029:                                                      }
3030:                                                  }
3031:                                                  ITAU = IR + LDWRKR * M;
3032:                                                  IWORK = ITAU + M;
3033:                                                  // *
3034:                                                  // *                    Compute A=L*Q
3035:                                                  // *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3036:                                                  // *
3037:                                                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3038:                                                                   , LWORK - IWORK + 1, ref IERR);
3039:                                                  // *
3040:                                                  // *                    Copy L to WORK(IU), zeroing out below it
3041:                                                  // *
3042:                                                  this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3043:                                                                   , LDWRKU);
3044:                                                  this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3045:                                                                   , LDWRKU);
3046:                                                  // *
3047:                                                  // *                    Generate Q in A
3048:                                                  // *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3049:                                                  // *
3050:                                                  this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
3051:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3052:                                                  IE = ITAU;
3053:                                                  ITAUQ = IE + M;
3054:                                                  ITAUP = ITAUQ + M;
3055:                                                  IWORK = ITAUP + M;
3056:                                                  // *
3057:                                                  // *                    Bidiagonalize L in WORK(IU), copying result to
3058:                                                  // *                    WORK(IR)
3059:                                                  // *                    (Workspace: need 2*M*M+4*M,
3060:                                                  // *                                prefer 2*M*M+3*M+2*M*NB)
3061:                                                  // *
3062:                                                  this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3063:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3064:                                                  this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
3065:                                                                   , LDWRKR);
3066:                                                  // *
3067:                                                  // *                    Generate right bidiagonalizing vectors in WORK(IU)
3068:                                                  // *                    (Workspace: need 2*M*M+4*M-1,
3069:                                                  // *                                prefer 2*M*M+3*M+(M-1)*NB)
3070:                                                  // *
3071:                                                  this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3072:                                                                   , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3073:                                                  // *
3074:                                                  // *                    Generate left bidiagonalizing vectors in WORK(IR)
3075:                                                  // *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
3076:                                                  // *
3077:                                                  this._dorgbr.Run("Q", M, M, M, ref WORK, IR + o_work, LDWRKR
3078:                                                                   , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3079:                                                  IWORK = IE + M;
3080:                                                  // *
3081:                                                  // *                    Perform bidiagonal QR iteration, computing left
3082:                                                  // *                    singular vectors of L in WORK(IR) and computing
3083:                                                  // *                    right singular vectors of L in WORK(IU)
3084:                                                  // *                    (Workspace: need 2*M*M+BDSPAC)
3085:                                                  // *
3086:                                                  this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3087:                                                                   , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
3088:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
3089:                                                  // *
3090:                                                  // *                    Multiply right singular vectors of L in WORK(IU) by
3091:                                                  // *                    Q in A, storing result in VT
3092:                                                  // *                    (Workspace: need M*M)
3093:                                                  // *
3094:                                                  this._dgemm.Run("N", "N", M, N, M, ONE
3095:                                                                  , WORK, IU + o_work, LDWRKU, A, offset_a, LDA, ZERO, ref VT, offset_vt
3096:                                                                  , LDVT);
3097:                                                  // *
3098:                                                  // *                    Copy left singular vectors of L to A
3099:                                                  // *                    (Workspace: need M*M)
3100:                                                  // *
3101:                                                  this._dlacpy.Run("F", M, M, WORK, IR + o_work, LDWRKR, ref A, offset_a
3102:                                                                   , LDA);
3103:                                                  // *
3104:                                              }
3105:                                              else
3106:                                              {
3107:                                                  // *
3108:                                                  // *                    Insufficient workspace for a fast algorithm
3109:                                                  // *
3110:                                                  ITAU = 1;
3111:                                                  IWORK = ITAU + M;
3112:                                                  // *
3113:                                                  // *                    Compute A=L*Q, copying result to VT
3114:                                                  // *                    (Workspace: need 2*M, prefer M+M*NB)
3115:                                                  // *
3116:                                                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3117:                                                                   , LWORK - IWORK + 1, ref IERR);
3118:                                                  this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3119:                                                                   , LDVT);
3120:                                                  // *
3121:                                                  // *                    Generate Q in VT
3122:                                                  // *                    (Workspace: need 2*M, prefer M+M*NB)
3123:                                                  // *
3124:                                                  this._dorglq.Run(M, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3125:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3126:                                                  IE = ITAU;
3127:                                                  ITAUQ = IE + M;
3128:                                                  ITAUP = ITAUQ + M;
3129:                                                  IWORK = ITAUP + M;
3130:                                                  // *
3131:                                                  // *                    Zero out above L in A
3132:                                                  // *
3133:                                                  this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
3134:                                                                   , LDA);
3135:                                                  // *
3136:                                                  // *                    Bidiagonalize L in A
3137:                                                  // *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
3138:                                                  // *
3139:                                                  this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3140:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3141:                                                  // *
3142:                                                  // *                    Multiply right vectors bidiagonalizing L by Q in VT
3143:                                                  // *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
3144:                                                  // *
3145:                                                  this._dormbr.Run("P", "L", "T", M, N, M
3146:                                                                   , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3147:                                                                   , LWORK - IWORK + 1, ref IERR);
3148:                                                  // *
3149:                                                  // *                    Generate left bidiagonalizing vectors of L in A
3150:                                                  // *                    (Workspace: need 4*M, prefer 3*M+M*NB)
3151:                                                  // *
3152:                                                  this._dorgbr.Run("Q", M, M, M, ref A, offset_a, LDA
3153:                                                                   , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3154:                                                  IWORK = IE + M;
3155:                                                  // *
3156:                                                  // *                    Perform bidiagonal QR iteration, compute left
3157:                                                  // *                    singular vectors of A in A and compute right
3158:                                                  // *                    singular vectors of A in VT
3159:                                                  // *                    (Workspace: need BDSPAC)
3160:                                                  // *
3161:                                                  this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3162:                                                                   , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
3163:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
3164:                                                  // *
3165:                                              }
3166:                                              // *
3167:                                          }
3168:                                          else
3169:                                          {
3170:                                              if (WNTUAS)
3171:                                              {
3172:                                                  // *
3173:                                                  // *                 Path 6t(N much larger than M, JOBU='S' or 'A',
3174:                                                  // *                         JOBVT='S')
3175:                                                  // *                 M right singular vectors to be computed in VT and
3176:                                                  // *                 M left singular vectors to be computed in U
3177:                                                  // *
3178:                                                  if (LWORK >= M * M + Math.Max(4 * M, BDSPAC))
3179:                                                  {
3180:                                                      // *
3181:                                                      // *                    Sufficient workspace for a fast algorithm
3182:                                                      // *
3183:                                                      IU = 1;
3184:                                                      if (LWORK >= WRKBL + LDA * M)
3185:                                                      {
3186:                                                          // *
3187:                                                          // *                       WORK(IU) is LDA by N
3188:                                                          // *
3189:                                                          LDWRKU = LDA;
3190:                                                      }
3191:                                                      else
3192:                                                      {
3193:                                                          // *
3194:                                                          // *                       WORK(IU) is LDA by M
3195:                                                          // *
3196:                                                          LDWRKU = M;
3197:                                                      }
3198:                                                      ITAU = IU + LDWRKU * M;
3199:                                                      IWORK = ITAU + M;
3200:                                                      // *
3201:                                                      // *                    Compute A=L*Q
3202:                                                      // *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3203:                                                      // *
3204:                                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3205:                                                                       , LWORK - IWORK + 1, ref IERR);
3206:                                                      // *
3207:                                                      // *                    Copy L to WORK(IU), zeroing out above it
3208:                                                      // *
3209:                                                      this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3210:                                                                       , LDWRKU);
3211:                                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3212:                                                                       , LDWRKU);
3213:                                                      // *
3214:                                                      // *                    Generate Q in A
3215:                                                      // *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3216:                                                      // *
3217:                                                      this._dorglq.Run(M, N, M, ref A, offset_a, LDA, WORK, ITAU + o_work
3218:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3219:                                                      IE = ITAU;
3220:                                                      ITAUQ = IE + M;
3221:                                                      ITAUP = ITAUQ + M;
3222:                                                      IWORK = ITAUP + M;
3223:                                                      // *
3224:                                                      // *                    Bidiagonalize L in WORK(IU), copying result to U
3225:                                                      // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3226:                                                      // *
3227:                                                      this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3228:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3229:                                                      this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref U, offset_u
3230:                                                                       , LDU);
3231:                                                      // *
3232:                                                      // *                    Generate right bidiagonalizing vectors in WORK(IU)
3233:                                                      // *                    (Workspace: need M*M+4*M-1,
3234:                                                      // *                                prefer M*M+3*M+(M-1)*NB)
3235:                                                      // *
3236:                                                      this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3237:                                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3238:                                                      // *
3239:                                                      // *                    Generate left bidiagonalizing vectors in U
3240:                                                      // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
3241:                                                      // *
3242:                                                      this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3243:                                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3244:                                                      IWORK = IE + M;
3245:                                                      // *
3246:                                                      // *                    Perform bidiagonal QR iteration, computing left
3247:                                                      // *                    singular vectors of L in U and computing right
3248:                                                      // *                    singular vectors of L in WORK(IU)
3249:                                                      // *                    (Workspace: need M*M+BDSPAC)
3250:                                                      // *
3251:                                                      this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3252:                                                                       , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref U, offset_u, LDU, ref DUM, offset_dum
3253:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
3254:                                                      // *
3255:                                                      // *                    Multiply right singular vectors of L in WORK(IU) by
3256:                                                      // *                    Q in A, storing result in VT
3257:                                                      // *                    (Workspace: need M*M)
3258:                                                      // *
3259:                                                      this._dgemm.Run("N", "N", M, N, M, ONE
3260:                                                                      , WORK, IU + o_work, LDWRKU, A, offset_a, LDA, ZERO, ref VT, offset_vt
3261:                                                                      , LDVT);
3262:                                                      // *
3263:                                                  }
3264:                                                  else
3265:                                                  {
3266:                                                      // *
3267:                                                      // *                    Insufficient workspace for a fast algorithm
3268:                                                      // *
3269:                                                      ITAU = 1;
3270:                                                      IWORK = ITAU + M;
3271:                                                      // *
3272:                                                      // *                    Compute A=L*Q, copying result to VT
3273:                                                      // *                    (Workspace: need 2*M, prefer M+M*NB)
3274:                                                      // *
3275:                                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3276:                                                                       , LWORK - IWORK + 1, ref IERR);
3277:                                                      this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3278:                                                                       , LDVT);
3279:                                                      // *
3280:                                                      // *                    Generate Q in VT
3281:                                                      // *                    (Workspace: need 2*M, prefer M+M*NB)
3282:                                                      // *
3283:                                                      this._dorglq.Run(M, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3284:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3285:                                                      // *
3286:                                                      // *                    Copy L to U, zeroing out above it
3287:                                                      // *
3288:                                                      this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
3289:                                                                       , LDU);
3290:                                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
3291:                                                                       , LDU);
3292:                                                      IE = ITAU;
3293:                                                      ITAUQ = IE + M;
3294:                                                      ITAUP = ITAUQ + M;
3295:                                                      IWORK = ITAUP + M;
3296:                                                      // *
3297:                                                      // *                    Bidiagonalize L in U
3298:                                                      // *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
3299:                                                      // *
3300:                                                      this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
3301:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3302:                                                      // *
3303:                                                      // *                    Multiply right bidiagonalizing vectors in U by Q
3304:                                                      // *                    in VT
3305:                                                      // *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
3306:                                                      // *
3307:                                                      this._dormbr.Run("P", "L", "T", M, N, M
3308:                                                                       , ref U, offset_u, LDU, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3309:                                                                       , LWORK - IWORK + 1, ref IERR);
3310:                                                      // *
3311:                                                      // *                    Generate left bidiagonalizing vectors in U
3312:                                                      // *                    (Workspace: need 4*M, prefer 3*M+M*NB)
3313:                                                      // *
3314:                                                      this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3315:                                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3316:                                                      IWORK = IE + M;
3317:                                                      // *
3318:                                                      // *                    Perform bidiagonal QR iteration, computing left
3319:                                                      // *                    singular vectors of A in U and computing right
3320:                                                      // *                    singular vectors of A in VT
3321:                                                      // *                    (Workspace: need BDSPAC)
3322:                                                      // *
3323:                                                      this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3324:                                                                       , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
3325:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
3326:                                                      // *
3327:                                                  }
3328:                                                  // *
3329:                                              }
3330:                                          }
3331:                                      }
3332:                                      // *
3333:                                  }
3334:                                  else
3335:                                  {
3336:                                      if (WNTVA)
3337:                                      {
3338:                                          // *
3339:                                          if (WNTUN)
3340:                                          {
3341:                                              // *
3342:                                              // *                 Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3343:                                              // *                 N right singular vectors to be computed in VT and
3344:                                              // *                 no left singular vectors to be computed
3345:                                              // *
3346:                                              if (LWORK >= M * M + Math.Max(N + M, Math.Max(4 * M, BDSPAC)))
3347:                                              {
3348:                                                  // *
3349:                                                  // *                    Sufficient workspace for a fast algorithm
3350:                                                  // *
3351:                                                  IR = 1;
3352:                                                  if (LWORK >= WRKBL + LDA * M)
3353:                                                  {
3354:                                                      // *
3355:                                                      // *                       WORK(IR) is LDA by M
3356:                                                      // *
3357:                                                      LDWRKR = LDA;
3358:                                                  }
3359:                                                  else
3360:                                                  {
3361:                                                      // *
3362:                                                      // *                       WORK(IR) is M by M
3363:                                                      // *
3364:                                                      LDWRKR = M;
3365:                                                  }
3366:                                                  ITAU = IR + LDWRKR * M;
3367:                                                  IWORK = ITAU + M;
3368:                                                  // *
3369:                                                  // *                    Compute A=L*Q, copying result to VT
3370:                                                  // *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3371:                                                  // *
3372:                                                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3373:                                                                   , LWORK - IWORK + 1, ref IERR);
3374:                                                  this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3375:                                                                   , LDVT);
3376:                                                  // *
3377:                                                  // *                    Copy L to WORK(IR), zeroing out above it
3378:                                                  // *
3379:                                                  this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IR + o_work
3380:                                                                   , LDWRKR);
3381:                                                  this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IR + LDWRKR + o_work
3382:                                                                   , LDWRKR);
3383:                                                  // *
3384:                                                  // *                    Generate Q in VT
3385:                                                  // *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
3386:                                                  // *
3387:                                                  this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3388:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3389:                                                  IE = ITAU;
3390:                                                  ITAUQ = IE + M;
3391:                                                  ITAUP = ITAUQ + M;
3392:                                                  IWORK = ITAUP + M;
3393:                                                  // *
3394:                                                  // *                    Bidiagonalize L in WORK(IR)
3395:                                                  // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3396:                                                  // *
3397:                                                  this._dgebrd.Run(M, M, ref WORK, IR + o_work, LDWRKR, ref S, offset_s, ref WORK, IE + o_work
3398:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3399:                                                  // *
3400:                                                  // *                    Generate right bidiagonalizing vectors in WORK(IR)
3401:                                                  // *                    (Workspace: need M*M+4*M-1,
3402:                                                  // *                                prefer M*M+3*M+(M-1)*NB)
3403:                                                  // *
3404:                                                  this._dorgbr.Run("P", M, M, M, ref WORK, IR + o_work, LDWRKR
3405:                                                                   , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3406:                                                  IWORK = IE + M;
3407:                                                  // *
3408:                                                  // *                    Perform bidiagonal QR iteration, computing right
3409:                                                  // *                    singular vectors of L in WORK(IR)
3410:                                                  // *                    (Workspace: need M*M+BDSPAC)
3411:                                                  // *
3412:                                                  this._dbdsqr.Run("U", M, M, 0, 0, ref S, offset_s
3413:                                                                   , ref WORK, IE + o_work, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum, 1, ref DUM, offset_dum
3414:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
3415:                                                  // *
3416:                                                  // *                    Multiply right singular vectors of L in WORK(IR) by
3417:                                                  // *                    Q in VT, storing result in A
3418:                                                  // *                    (Workspace: need M*M)
3419:                                                  // *
3420:                                                  this._dgemm.Run("N", "N", M, N, M, ONE
3421:                                                                  , WORK, IR + o_work, LDWRKR, VT, offset_vt, LDVT, ZERO, ref A, offset_a
3422:                                                                  , LDA);
3423:                                                  // *
3424:                                                  // *                    Copy right singular vectors of A from A to VT
3425:                                                  // *
3426:                                                  this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
3427:                                                                   , LDVT);
3428:                                                  // *
3429:                                              }
3430:                                              else
3431:                                              {
3432:                                                  // *
3433:                                                  // *                    Insufficient workspace for a fast algorithm
3434:                                                  // *
3435:                                                  ITAU = 1;
3436:                                                  IWORK = ITAU + M;
3437:                                                  // *
3438:                                                  // *                    Compute A=L*Q, copying result to VT
3439:                                                  // *                    (Workspace: need 2*M, prefer M+M*NB)
3440:                                                  // *
3441:                                                  this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3442:                                                                   , LWORK - IWORK + 1, ref IERR);
3443:                                                  this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3444:                                                                   , LDVT);
3445:                                                  // *
3446:                                                  // *                    Generate Q in VT
3447:                                                  // *                    (Workspace: need M+N, prefer M+N*NB)
3448:                                                  // *
3449:                                                  this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3450:                                                                   , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3451:                                                  IE = ITAU;
3452:                                                  ITAUQ = IE + M;
3453:                                                  ITAUP = ITAUQ + M;
3454:                                                  IWORK = ITAUP + M;
3455:                                                  // *
3456:                                                  // *                    Zero out above L in A
3457:                                                  // *
3458:                                                  this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
3459:                                                                   , LDA);
3460:                                                  // *
3461:                                                  // *                    Bidiagonalize L in A
3462:                                                  // *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
3463:                                                  // *
3464:                                                  this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3465:                                                                   , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3466:                                                  // *
3467:                                                  // *                    Multiply right bidiagonalizing vectors in A by Q
3468:                                                  // *                    in VT
3469:                                                  // *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
3470:                                                  // *
3471:                                                  this._dormbr.Run("P", "L", "T", M, N, M
3472:                                                                   , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3473:                                                                   , LWORK - IWORK + 1, ref IERR);
3474:                                                  IWORK = IE + M;
3475:                                                  // *
3476:                                                  // *                    Perform bidiagonal QR iteration, computing right
3477:                                                  // *                    singular vectors of A in VT
3478:                                                  // *                    (Workspace: need BDSPAC)
3479:                                                  // *
3480:                                                  this._dbdsqr.Run("U", M, N, 0, 0, ref S, offset_s
3481:                                                                   , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref DUM, offset_dum, 1, ref DUM, offset_dum
3482:                                                                   , 1, ref WORK, IWORK + o_work, ref INFO);
3483:                                                  // *
3484:                                              }
3485:                                              // *
3486:                                          }
3487:                                          else
3488:                                          {
3489:                                              if (WNTUO)
3490:                                              {
3491:                                                  // *
3492:                                                  // *                 Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3493:                                                  // *                 N right singular vectors to be computed in VT and
3494:                                                  // *                 M left singular vectors to be overwritten on A
3495:                                                  // *
3496:                                                  if (LWORK >= 2 * M * M + Math.Max(N + M, Math.Max(4 * M, BDSPAC)))
3497:                                                  {
3498:                                                      // *
3499:                                                      // *                    Sufficient workspace for a fast algorithm
3500:                                                      // *
3501:                                                      IU = 1;
3502:                                                      if (LWORK >= WRKBL + 2 * LDA * M)
3503:                                                      {
3504:                                                          // *
3505:                                                          // *                       WORK(IU) is LDA by M and WORK(IR) is LDA by M
3506:                                                          // *
3507:                                                          LDWRKU = LDA;
3508:                                                          IR = IU + LDWRKU * M;
3509:                                                          LDWRKR = LDA;
3510:                                                      }
3511:                                                      else
3512:                                                      {
3513:                                                          if (LWORK >= WRKBL + (LDA + M) * M)
3514:                                                          {
3515:                                                              // *
3516:                                                              // *                       WORK(IU) is LDA by M and WORK(IR) is M by M
3517:                                                              // *
3518:                                                              LDWRKU = LDA;
3519:                                                              IR = IU + LDWRKU * M;
3520:                                                              LDWRKR = M;
3521:                                                          }
3522:                                                          else
3523:                                                          {
3524:                                                              // *
3525:                                                              // *                       WORK(IU) is M by M and WORK(IR) is M by M
3526:                                                              // *
3527:                                                              LDWRKU = M;
3528:                                                              IR = IU + LDWRKU * M;
3529:                                                              LDWRKR = M;
3530:                                                          }
3531:                                                      }
3532:                                                      ITAU = IR + LDWRKR * M;
3533:                                                      IWORK = ITAU + M;
3534:                                                      // *
3535:                                                      // *                    Compute A=L*Q, copying result to VT
3536:                                                      // *                    (Workspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3537:                                                      // *
3538:                                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3539:                                                                       , LWORK - IWORK + 1, ref IERR);
3540:                                                      this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3541:                                                                       , LDVT);
3542:                                                      // *
3543:                                                      // *                    Generate Q in VT
3544:                                                      // *                    (Workspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3545:                                                      // *
3546:                                                      this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3547:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3548:                                                      // *
3549:                                                      // *                    Copy L to WORK(IU), zeroing out above it
3550:                                                      // *
3551:                                                      this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3552:                                                                       , LDWRKU);
3553:                                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3554:                                                                       , LDWRKU);
3555:                                                      IE = ITAU;
3556:                                                      ITAUQ = IE + M;
3557:                                                      ITAUP = ITAUQ + M;
3558:                                                      IWORK = ITAUP + M;
3559:                                                      // *
3560:                                                      // *                    Bidiagonalize L in WORK(IU), copying result to
3561:                                                      // *                    WORK(IR)
3562:                                                      // *                    (Workspace: need 2*M*M+4*M,
3563:                                                      // *                                prefer 2*M*M+3*M+2*M*NB)
3564:                                                      // *
3565:                                                      this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3566:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3567:                                                      this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work
3568:                                                                       , LDWRKR);
3569:                                                      // *
3570:                                                      // *                    Generate right bidiagonalizing vectors in WORK(IU)
3571:                                                      // *                    (Workspace: need 2*M*M+4*M-1,
3572:                                                      // *                                prefer 2*M*M+3*M+(M-1)*NB)
3573:                                                      // *
3574:                                                      this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3575:                                                                       , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3576:                                                      // *
3577:                                                      // *                    Generate left bidiagonalizing vectors in WORK(IR)
3578:                                                      // *                    (Workspace: need 2*M*M+4*M, prefer 2*M*M+3*M+M*NB)
3579:                                                      // *
3580:                                                      this._dorgbr.Run("Q", M, M, M, ref WORK, IR + o_work, LDWRKR
3581:                                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3582:                                                      IWORK = IE + M;
3583:                                                      // *
3584:                                                      // *                    Perform bidiagonal QR iteration, computing left
3585:                                                      // *                    singular vectors of L in WORK(IR) and computing
3586:                                                      // *                    right singular vectors of L in WORK(IU)
3587:                                                      // *                    (Workspace: need 2*M*M+BDSPAC)
3588:                                                      // *
3589:                                                      this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3590:                                                                       , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref WORK, IR + o_work, LDWRKR, ref DUM, offset_dum
3591:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
3592:                                                      // *
3593:                                                      // *                    Multiply right singular vectors of L in WORK(IU) by
3594:                                                      // *                    Q in VT, storing result in A
3595:                                                      // *                    (Workspace: need M*M)
3596:                                                      // *
3597:                                                      this._dgemm.Run("N", "N", M, N, M, ONE
3598:                                                                      , WORK, IU + o_work, LDWRKU, VT, offset_vt, LDVT, ZERO, ref A, offset_a
3599:                                                                      , LDA);
3600:                                                      // *
3601:                                                      // *                    Copy right singular vectors of A from A to VT
3602:                                                      // *
3603:                                                      this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
3604:                                                                       , LDVT);
3605:                                                      // *
3606:                                                      // *                    Copy left singular vectors of A from WORK(IR) to A
3607:                                                      // *
3608:                                                      this._dlacpy.Run("F", M, M, WORK, IR + o_work, LDWRKR, ref A, offset_a
3609:                                                                       , LDA);
3610:                                                      // *
3611:                                                  }
3612:                                                  else
3613:                                                  {
3614:                                                      // *
3615:                                                      // *                    Insufficient workspace for a fast algorithm
3616:                                                      // *
3617:                                                      ITAU = 1;
3618:                                                      IWORK = ITAU + M;
3619:                                                      // *
3620:                                                      // *                    Compute A=L*Q, copying result to VT
3621:                                                      // *                    (Workspace: need 2*M, prefer M+M*NB)
3622:                                                      // *
3623:                                                      this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3624:                                                                       , LWORK - IWORK + 1, ref IERR);
3625:                                                      this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3626:                                                                       , LDVT);
3627:                                                      // *
3628:                                                      // *                    Generate Q in VT
3629:                                                      // *                    (Workspace: need M+N, prefer M+N*NB)
3630:                                                      // *
3631:                                                      this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3632:                                                                       , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3633:                                                      IE = ITAU;
3634:                                                      ITAUQ = IE + M;
3635:                                                      ITAUP = ITAUQ + M;
3636:                                                      IWORK = ITAUP + M;
3637:                                                      // *
3638:                                                      // *                    Zero out above L in A
3639:                                                      // *
3640:                                                      this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref A, 1+2 * LDA + o_a
3641:                                                                       , LDA);
3642:                                                      // *
3643:                                                      // *                    Bidiagonalize L in A
3644:                                                      // *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
3645:                                                      // *
3646:                                                      this._dgebrd.Run(M, M, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3647:                                                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3648:                                                      // *
3649:                                                      // *                    Multiply right bidiagonalizing vectors in A by Q
3650:                                                      // *                    in VT
3651:                                                      // *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
3652:                                                      // *
3653:                                                      this._dormbr.Run("P", "L", "T", M, N, M
3654:                                                                       , ref A, offset_a, LDA, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3655:                                                                       , LWORK - IWORK + 1, ref IERR);
3656:                                                      // *
3657:                                                      // *                    Generate left bidiagonalizing vectors in A
3658:                                                      // *                    (Workspace: need 4*M, prefer 3*M+M*NB)
3659:                                                      // *
3660:                                                      this._dorgbr.Run("Q", M, M, M, ref A, offset_a, LDA
3661:                                                                       , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3662:                                                      IWORK = IE + M;
3663:                                                      // *
3664:                                                      // *                    Perform bidiagonal QR iteration, computing left
3665:                                                      // *                    singular vectors of A in A and computing right
3666:                                                      // *                    singular vectors of A in VT
3667:                                                      // *                    (Workspace: need BDSPAC)
3668:                                                      // *
3669:                                                      this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3670:                                                                       , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
3671:                                                                       , 1, ref WORK, IWORK + o_work, ref INFO);
3672:                                                      // *
3673:                                                  }
3674:                                                  // *
3675:                                              }
3676:                                              else
3677:                                              {
3678:                                                  if (WNTUAS)
3679:                                                  {
3680:                                                      // *
3681:                                                      // *                 Path 9t(N much larger than M, JOBU='S' or 'A',
3682:                                                      // *                         JOBVT='A')
3683:                                                      // *                 N right singular vectors to be computed in VT and
3684:                                                      // *                 M left singular vectors to be computed in U
3685:                                                      // *
3686:                                                      if (LWORK >= M * M + Math.Max(N + M, Math.Max(4 * M, BDSPAC)))
3687:                                                      {
3688:                                                          // *
3689:                                                          // *                    Sufficient workspace for a fast algorithm
3690:                                                          // *
3691:                                                          IU = 1;
3692:                                                          if (LWORK >= WRKBL + LDA * M)
3693:                                                          {
3694:                                                              // *
3695:                                                              // *                       WORK(IU) is LDA by M
3696:                                                              // *
3697:                                                              LDWRKU = LDA;
3698:                                                          }
3699:                                                          else
3700:                                                          {
3701:                                                              // *
3702:                                                              // *                       WORK(IU) is M by M
3703:                                                              // *
3704:                                                              LDWRKU = M;
3705:                                                          }
3706:                                                          ITAU = IU + LDWRKU * M;
3707:                                                          IWORK = ITAU + M;
3708:                                                          // *
3709:                                                          // *                    Compute A=L*Q, copying result to VT
3710:                                                          // *                    (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
3711:                                                          // *
3712:                                                          this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3713:                                                                           , LWORK - IWORK + 1, ref IERR);
3714:                                                          this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3715:                                                                           , LDVT);
3716:                                                          // *
3717:                                                          // *                    Generate Q in VT
3718:                                                          // *                    (Workspace: need M*M+M+N, prefer M*M+M+N*NB)
3719:                                                          // *
3720:                                                          this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3721:                                                                           , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3722:                                                          // *
3723:                                                          // *                    Copy L to WORK(IU), zeroing out above it
3724:                                                          // *
3725:                                                          this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref WORK, IU + o_work
3726:                                                                           , LDWRKU);
3727:                                                          this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref WORK, IU + LDWRKU + o_work
3728:                                                                           , LDWRKU);
3729:                                                          IE = ITAU;
3730:                                                          ITAUQ = IE + M;
3731:                                                          ITAUP = ITAUQ + M;
3732:                                                          IWORK = ITAUP + M;
3733:                                                          // *
3734:                                                          // *                    Bidiagonalize L in WORK(IU), copying result to U
3735:                                                          // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
3736:                                                          // *
3737:                                                          this._dgebrd.Run(M, M, ref WORK, IU + o_work, LDWRKU, ref S, offset_s, ref WORK, IE + o_work
3738:                                                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3739:                                                          this._dlacpy.Run("L", M, M, WORK, IU + o_work, LDWRKU, ref U, offset_u
3740:                                                                           , LDU);
3741:                                                          // *
3742:                                                          // *                    Generate right bidiagonalizing vectors in WORK(IU)
3743:                                                          // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+(M-1)*NB)
3744:                                                          // *
3745:                                                          this._dorgbr.Run("P", M, M, M, ref WORK, IU + o_work, LDWRKU
3746:                                                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3747:                                                          // *
3748:                                                          // *                    Generate left bidiagonalizing vectors in U
3749:                                                          // *                    (Workspace: need M*M+4*M, prefer M*M+3*M+M*NB)
3750:                                                          // *
3751:                                                          this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3752:                                                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3753:                                                          IWORK = IE + M;
3754:                                                          // *
3755:                                                          // *                    Perform bidiagonal QR iteration, computing left
3756:                                                          // *                    singular vectors of L in U and computing right
3757:                                                          // *                    singular vectors of L in WORK(IU)
3758:                                                          // *                    (Workspace: need M*M+BDSPAC)
3759:                                                          // *
3760:                                                          this._dbdsqr.Run("U", M, M, M, 0, ref S, offset_s
3761:                                                                           , ref WORK, IE + o_work, ref WORK, IU + o_work, LDWRKU, ref U, offset_u, LDU, ref DUM, offset_dum
3762:                                                                           , 1, ref WORK, IWORK + o_work, ref INFO);
3763:                                                          // *
3764:                                                          // *                    Multiply right singular vectors of L in WORK(IU) by
3765:                                                          // *                    Q in VT, storing result in A
3766:                                                          // *                    (Workspace: need M*M)
3767:                                                          // *
3768:                                                          this._dgemm.Run("N", "N", M, N, M, ONE
3769:                                                                          , WORK, IU + o_work, LDWRKU, VT, offset_vt, LDVT, ZERO, ref A, offset_a
3770:                                                                          , LDA);
3771:                                                          // *
3772:                                                          // *                    Copy right singular vectors of A from A to VT
3773:                                                          // *
3774:                                                          this._dlacpy.Run("F", M, N, A, offset_a, LDA, ref VT, offset_vt
3775:                                                                           , LDVT);
3776:                                                          // *
3777:                                                      }
3778:                                                      else
3779:                                                      {
3780:                                                          // *
3781:                                                          // *                    Insufficient workspace for a fast algorithm
3782:                                                          // *
3783:                                                          ITAU = 1;
3784:                                                          IWORK = ITAU + M;
3785:                                                          // *
3786:                                                          // *                    Compute A=L*Q, copying result to VT
3787:                                                          // *                    (Workspace: need 2*M, prefer M+M*NB)
3788:                                                          // *
3789:                                                          this._dgelqf.Run(M, N, ref A, offset_a, LDA, ref WORK, ITAU + o_work, ref WORK, IWORK + o_work
3790:                                                                           , LWORK - IWORK + 1, ref IERR);
3791:                                                          this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3792:                                                                           , LDVT);
3793:                                                          // *
3794:                                                          // *                    Generate Q in VT
3795:                                                          // *                    (Workspace: need M+N, prefer M+N*NB)
3796:                                                          // *
3797:                                                          this._dorglq.Run(N, N, M, ref VT, offset_vt, LDVT, WORK, ITAU + o_work
3798:                                                                           , ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3799:                                                          // *
3800:                                                          // *                    Copy L to U, zeroing out above it
3801:                                                          // *
3802:                                                          this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
3803:                                                                           , LDU);
3804:                                                          this._dlaset.Run("U", M - 1, M - 1, ZERO, ZERO, ref U, 1+2 * LDU + o_u
3805:                                                                           , LDU);
3806:                                                          IE = ITAU;
3807:                                                          ITAUQ = IE + M;
3808:                                                          ITAUP = ITAUQ + M;
3809:                                                          IWORK = ITAUP + M;
3810:                                                          // *
3811:                                                          // *                    Bidiagonalize L in U
3812:                                                          // *                    (Workspace: need 4*M, prefer 3*M+2*M*NB)
3813:                                                          // *
3814:                                                          this._dgebrd.Run(M, M, ref U, offset_u, LDU, ref S, offset_s, ref WORK, IE + o_work
3815:                                                                           , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3816:                                                          // *
3817:                                                          // *                    Multiply right bidiagonalizing vectors in U by Q
3818:                                                          // *                    in VT
3819:                                                          // *                    (Workspace: need 3*M+N, prefer 3*M+N*NB)
3820:                                                          // *
3821:                                                          this._dormbr.Run("P", "L", "T", M, N, M
3822:                                                                           , ref U, offset_u, LDU, WORK, ITAUP + o_work, ref VT, offset_vt, LDVT, ref WORK, IWORK + o_work
3823:                                                                           , LWORK - IWORK + 1, ref IERR);
3824:                                                          // *
3825:                                                          // *                    Generate left bidiagonalizing vectors in U
3826:                                                          // *                    (Workspace: need 4*M, prefer 3*M+M*NB)
3827:                                                          // *
3828:                                                          this._dorgbr.Run("Q", M, M, M, ref U, offset_u, LDU
3829:                                                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3830:                                                          IWORK = IE + M;
3831:                                                          // *
3832:                                                          // *                    Perform bidiagonal QR iteration, computing left
3833:                                                          // *                    singular vectors of A in U and computing right
3834:                                                          // *                    singular vectors of A in VT
3835:                                                          // *                    (Workspace: need BDSPAC)
3836:                                                          // *
3837:                                                          this._dbdsqr.Run("U", M, N, M, 0, ref S, offset_s
3838:                                                                           , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
3839:                                                                           , 1, ref WORK, IWORK + o_work, ref INFO);
3840:                                                          // *
3841:                                                      }
3842:                                                      // *
3843:                                                  }
3844:                                              }
3845:                                          }
3846:                                          // *
3847:                                      }
3848:                                  }
3849:                              }
3850:                          }
3851:                      }
3852:                      // *
3853:                  }
3854:                  else
3855:                  {
3856:                      // *
3857:                      // *           N .LT. MNTHR
3858:                      // *
3859:                      // *           Path 10t(N greater than M, but not much larger)
3860:                      // *           Reduce to bidiagonal form without LQ decomposition
3861:                      // *
3862:                      IE = 1;
3863:                      ITAUQ = IE + M;
3864:                      ITAUP = ITAUQ + M;
3865:                      IWORK = ITAUP + M;
3866:                      // *
3867:                      // *           Bidiagonalize A
3868:                      // *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
3869:                      // *
3870:                      this._dgebrd.Run(M, N, ref A, offset_a, LDA, ref S, offset_s, ref WORK, IE + o_work
3871:                                       , ref WORK, ITAUQ + o_work, ref WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3872:                      if (WNTUAS)
3873:                      {
3874:                          // *
3875:                          // *              If left singular vectors desired in U, copy result to U
3876:                          // *              and generate left bidiagonalizing vectors in U
3877:                          // *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
3878:                          // *
3879:                          this._dlacpy.Run("L", M, M, A, offset_a, LDA, ref U, offset_u
3880:                                           , LDU);
3881:                          this._dorgbr.Run("Q", M, M, N, ref U, offset_u, LDU
3882:                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3883:                      }
3884:                      if (WNTVAS)
3885:                      {
3886:                          // *
3887:                          // *              If right singular vectors desired in VT, copy result to
3888:                          // *              VT and generate right bidiagonalizing vectors in VT
3889:                          // *              (Workspace: need 3*M+NRVT, prefer 3*M+NRVT*NB)
3890:                          // *
3891:                          this._dlacpy.Run("U", M, N, A, offset_a, LDA, ref VT, offset_vt
3892:                                           , LDVT);
3893:                          if (WNTVA) NRVT = N;
3894:                          if (WNTVS) NRVT = M;
3895:                          this._dorgbr.Run("P", NRVT, N, M, ref VT, offset_vt, LDVT
3896:                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3897:                      }
3898:                      if (WNTUO)
3899:                      {
3900:                          // *
3901:                          // *              If left singular vectors desired in A, generate left
3902:                          // *              bidiagonalizing vectors in A
3903:                          // *              (Workspace: need 4*M-1, prefer 3*M+(M-1)*NB)
3904:                          // *
3905:                          this._dorgbr.Run("Q", M, M, N, ref A, offset_a, LDA
3906:                                           , WORK, ITAUQ + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3907:                      }
3908:                      if (WNTVO)
3909:                      {
3910:                          // *
3911:                          // *              If right singular vectors desired in A, generate right
3912:                          // *              bidiagonalizing vectors in A
3913:                          // *              (Workspace: need 4*M, prefer 3*M+M*NB)
3914:                          // *
3915:                          this._dorgbr.Run("P", M, N, M, ref A, offset_a, LDA
3916:                                           , WORK, ITAUP + o_work, ref WORK, IWORK + o_work, LWORK - IWORK + 1, ref IERR);
3917:                      }
3918:                      IWORK = IE + M;
3919:                      if (WNTUAS || WNTUO) NRU = M;
3920:                      if (WNTUN) NRU = 0;
3921:                      if (WNTVAS || WNTVO) NCVT = N;
3922:                      if (WNTVN) NCVT = 0;
3923:                      if ((!WNTUO) && (!WNTVO))
3924:                      {
3925:                          // *
3926:                          // *              Perform bidiagonal QR iteration, if desired, computing
3927:                          // *              left singular vectors in U and computing right singular
3928:                          // *              vectors in VT
3929:                          // *              (Workspace: need BDSPAC)
3930:                          // *
3931:                          this._dbdsqr.Run("L", M, NCVT, NRU, 0, ref S, offset_s
3932:                                           , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref U, offset_u, LDU, ref DUM, offset_dum
3933:                                           , 1, ref WORK, IWORK + o_work, ref INFO);
3934:                      }
3935:                      else
3936:                      {
3937:                          if ((!WNTUO) && WNTVO)
3938:                          {
3939:                              // *
3940:                              // *              Perform bidiagonal QR iteration, if desired, computing
3941:                              // *              left singular vectors in U and computing right singular
3942:                              // *              vectors in A
3943:                              // *              (Workspace: need BDSPAC)
3944:                              // *
3945:                              this._dbdsqr.Run("L", M, NCVT, NRU, 0, ref S, offset_s
3946:                                               , ref WORK, IE + o_work, ref A, offset_a, LDA, ref U, offset_u, LDU, ref DUM, offset_dum
3947:                                               , 1, ref WORK, IWORK + o_work, ref INFO);
3948:                          }
3949:                          else
3950:                          {
3951:                              // *
3952:                              // *              Perform bidiagonal QR iteration, if desired, computing
3953:                              // *              left singular vectors in A and computing right singular
3954:                              // *              vectors in VT
3955:                              // *              (Workspace: need BDSPAC)
3956:                              // *
3957:                              this._dbdsqr.Run("L", M, NCVT, NRU, 0, ref S, offset_s
3958:                                               , ref WORK, IE + o_work, ref VT, offset_vt, LDVT, ref A, offset_a, LDA, ref DUM, offset_dum
3959:                                               , 1, ref WORK, IWORK + o_work, ref INFO);
3960:                          }
3961:                      }
3962:                      // *
3963:                  }
3964:                  // *
3965:              }
3966:              // *
3967:              // *     If DBDSQR failed to converge, copy unconverged superdiagonals
3968:              // *     to WORK( 2:MINMN )
3969:              // *
3970:              if (INFO != 0)
3971:              {
3972:                  if (IE > 2)
3973:                  {
3974:                      for (I = 1; I <= MINMN - 1; I++)
3975:                      {
3976:                          WORK[I + 1 + o_work] = WORK[I + IE - 1 + o_work];
3977:                      }
3978:                  }
3979:                  if (IE < 2)
3980:                  {
3981:                      for (I = MINMN - 1; I >= 1; I +=  - 1)
3982:                      {
3983:                          WORK[I + 1 + o_work] = WORK[I + IE - 1 + o_work];
3984:                      }
3985:                  }
3986:              }
3987:              // *
3988:              // *     Undo scaling if necessary
3989:              // *
3990:              if (ISCL == 1)
3991:              {
3992:                  if (ANRM > BIGNUM)
3993:                  {
3994:                      this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN
3995:                                       , 1, ref S, offset_s, MINMN, ref IERR);
3996:                  }
3997:                  if (INFO != 0 && ANRM > BIGNUM)
3998:                  {
3999:                      this._dlascl.Run("G", 0, 0, BIGNUM, ANRM, MINMN - 1
4000:                                       , 1, ref WORK, 2 + o_work, MINMN, ref IERR);
4001:                  }
4002:                  if (ANRM < SMLNUM)
4003:                  {
4004:                      this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN
4005:                                       , 1, ref S, offset_s, MINMN, ref IERR);
4006:                  }
4007:                  if (INFO != 0 && ANRM < SMLNUM)
4008:                  {
4009:                      this._dlascl.Run("G", 0, 0, SMLNUM, ANRM, MINMN - 1
4010:                                       , 1, ref WORK, 2 + o_work, MINMN, ref IERR);
4011:                  }
4012:              }
4013:              // *
4014:              // *     Return optimal workspace in WORK(1)
4015:              // *
4016:              WORK[1 + o_work] = MAXWRK;
4017:              // *
4018:              return;
4019:              // *
4020:              // *     End of DGESVD
4021:              // *
4022:   
4023:              #endregion
4024:   
4025:          }
4026:      }
4027:  }