`   1:  #region Translated by Jose Antonio De Santiago-Castillo.`
`   2:   `
`   3:  //Translated by Jose Antonio De Santiago-Castillo. `
`   4:  //E-mail:JAntonioDeSantiago@gmail.com`
`   5:  //Web: www.DotNumerics.com`
`   6:  //`
`   7:  //Fortran to C# Translation.`
`   8:  //Translated by:`
`   9:  //F2CSharp Version 0.71 (November 10, 2009)`
`  10:  //Code Optimizations: None`
`  11:  //`
`  12:  #endregion`
`  13:   `
`  14:  using System;`
`  15:  using DotNumerics.FortranLibrary;`
`  16:   `
`  17:  namespace DotNumerics.CSLapack`
`  18:  {`
`  19:      /// <summary>`
`  20:      /// -- LAPACK routine (version 3.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// November 2006`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DBDSDC computes the singular value decomposition (SVD) of a real`
`  27:      /// N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,`
`  28:      /// using a divide and conquer method, where S is a diagonal matrix`
`  29:      /// with non-negative diagonal elements (the singular values of B), and`
`  30:      /// U and VT are orthogonal matrices of left and right singular vectors,`
`  31:      /// respectively. DBDSDC can be used to compute all singular values,`
`  32:      /// and optionally, singular vectors or singular vectors in compact form.`
`  33:      /// `
`  34:      /// This code makes very mild assumptions about floating point`
`  35:      /// arithmetic. It will work on machines with a guard digit in`
`  36:      /// add/subtract, or on those binary machines without guard digits`
`  37:      /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.`
`  38:      /// It could conceivably fail on hexadecimal or decimal machines`
`  39:      /// without guard digits, but we know of none.  See DLASD3 for details.`
`  40:      /// `
`  41:      /// The code currently calls DLASDQ if singular values only are desired.`
`  42:      /// However, it can be slightly modified to compute singular values`
`  43:      /// using the divide and conquer method.`
`  44:      /// `
`  45:      ///</summary>`
`  46:      public class DBDSDC`
`  47:      {`
`  48:      `
`  49:   `
`  50:          #region Dependencies`
`  51:          `
`  52:          LSAME _lsame; ILAENV _ilaenv; DLAMCH _dlamch; DLANST _dlanst; DCOPY _dcopy; DLARTG _dlartg; DLASCL _dlascl; `
`  53:          DLASD0 _dlasd0;DLASDA _dlasda; DLASDQ _dlasdq; DLASET _dlaset; DLASR _dlasr; DSWAP _dswap; XERBLA _xerbla; `
`  54:   `
`  55:          #endregion`
`  56:   `
`  57:   `
`  58:          #region Fields`
`  59:          `
`  60:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double TWO = 2.0E+0; int DIFL = 0; int DIFR = 0; `
`  61:          int GIVCOL = 0;int GIVNUM = 0; int GIVPTR = 0; int I = 0; int IC = 0; int ICOMPQ = 0; int IERR = 0; int II = 0; `
`  62:          int IS = 0;int IU = 0; int IUPLO = 0; int IVT = 0; int J = 0; int K = 0; int KK = 0; int MLVL = 0; int NM1 = 0; `
`  63:          int NSIZE = 0;int PERM = 0; int POLES = 0; int QSTART = 0; int SMLSIZ = 0; int SMLSZP = 0; int SQRE = 0; int START = 0; `
`  64:          int WSTART = 0;int Z = 0; double CS = 0; double EPS = 0; double ORGNRM = 0; double P = 0; double R = 0; double SN = 0; `
`  65:   `
`  66:          #endregion`
`  67:   `
`  68:          public DBDSDC(LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANST dlanst, DCOPY dcopy, DLARTG dlartg, DLASCL dlascl, DLASD0 dlasd0, DLASDA dlasda, DLASDQ dlasdq`
`  69:                        , DLASET dlaset, DLASR dlasr, DSWAP dswap, XERBLA xerbla)`
`  70:          {`
`  71:      `
`  72:   `
`  73:              #region Set Dependencies`
`  74:              `
`  75:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy; `
`  76:              this._dlartg = dlartg;this._dlascl = dlascl; this._dlasd0 = dlasd0; this._dlasda = dlasda; this._dlasdq = dlasdq; `
`  77:              this._dlaset = dlaset;this._dlasr = dlasr; this._dswap = dswap; this._xerbla = xerbla; `
`  78:   `
`  79:              #endregion`
`  80:   `
`  81:          }`
`  82:      `
`  83:          public DBDSDC()`
`  84:          {`
`  85:      `
`  86:   `
`  87:              #region Dependencies (Initialization)`
`  88:              `
`  89:              LSAME lsame = new LSAME();`
`  90:              IEEECK ieeeck = new IEEECK();`
`  91:              IPARMQ iparmq = new IPARMQ();`
`  92:              DLAMC3 dlamc3 = new DLAMC3();`
`  93:              DLASSQ dlassq = new DLASSQ();`
`  94:              DCOPY dcopy = new DCOPY();`
`  95:              XERBLA xerbla = new XERBLA();`
`  96:              DLAMRG dlamrg = new DLAMRG();`
`  97:              DLAPY2 dlapy2 = new DLAPY2();`
`  98:              DROT drot = new DROT();`
`  99:              DNRM2 dnrm2 = new DNRM2();`
` 100:              DLASD5 dlasd5 = new DLASD5();`
` 101:              DLAS2 dlas2 = new DLAS2();`
` 102:              DLASQ5 dlasq5 = new DLASQ5();`
` 103:              DLAZQ4 dlazq4 = new DLAZQ4();`
` 104:              DSCAL dscal = new DSCAL();`
` 105:              DSWAP dswap = new DSWAP();`
` 106:              DLASDT dlasdt = new DLASDT();`
` 107:              DDOT ddot = new DDOT();`
` 108:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 109:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 110:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 111:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 112:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 113:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 114:              DLANST dlanst = new DLANST(lsame, dlassq);`
` 115:              DLARTG dlartg = new DLARTG(dlamch);`
` 116:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);`
` 117:              DLACPY dlacpy = new DLACPY(lsame);`
` 118:              DLASET dlaset = new DLASET(lsame);`
` 119:              DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);`
` 120:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
` 121:              DLAED6 dlaed6 = new DLAED6(dlamch);`
` 122:              DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);`
` 123:              DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);`
` 124:              DLASD1 dlasd1 = new DLASD1(dlamrg, dlascl, dlasd2, dlasd3, xerbla);`
` 125:              DLASQ6 dlasq6 = new DLASQ6(dlamch);`
` 126:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);`
` 127:              DLASRT dlasrt = new DLASRT(lsame, xerbla);`
` 128:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);`
` 129:              DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);`
` 130:              DLASR dlasr = new DLASR(lsame, xerbla);`
` 131:              DLASV2 dlasv2 = new DLASV2(dlamch);`
` 132:              DBDSQR dbdsqr = new DBDSQR(lsame, dlamch, dlartg, dlas2, dlasq1, dlasr, dlasv2, drot, dscal, dswap`
` 133:                                         , xerbla);`
` 134:              DLASDQ dlasdq = new DLASDQ(dbdsqr, dlartg, dlasr, dswap, xerbla, lsame);`
` 135:              DLASD0 dlasd0 = new DLASD0(dlasd1, dlasdq, dlasdt, xerbla);`
` 136:              DLASD7 dlasd7 = new DLASD7(dcopy, dlamrg, drot, xerbla, dlamch, dlapy2);`
` 137:              DLASD8 dlasd8 = new DLASD8(dcopy, dlascl, dlasd4, dlaset, xerbla, ddot, dlamc3, dnrm2);`
` 138:              DLASD6 dlasd6 = new DLASD6(dcopy, dlamrg, dlascl, dlasd7, dlasd8, xerbla);`
` 139:              DLASDA dlasda = new DLASDA(dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla);`
` 140:   `
` 141:              #endregion`
` 142:   `
` 143:   `
` 144:              #region Set Dependencies`
` 145:              `
` 146:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dcopy = dcopy; `
` 147:              this._dlartg = dlartg;this._dlascl = dlascl; this._dlasd0 = dlasd0; this._dlasda = dlasda; this._dlasdq = dlasdq; `
` 148:              this._dlaset = dlaset;this._dlasr = dlasr; this._dswap = dswap; this._xerbla = xerbla; `
` 149:   `
` 150:              #endregion`
` 151:   `
` 152:          }`
` 153:          /// <summary>`
` 154:          /// Purpose`
` 155:          /// =======`
` 156:          /// `
` 157:          /// DBDSDC computes the singular value decomposition (SVD) of a real`
` 158:          /// N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,`
` 159:          /// using a divide and conquer method, where S is a diagonal matrix`
` 160:          /// with non-negative diagonal elements (the singular values of B), and`
` 161:          /// U and VT are orthogonal matrices of left and right singular vectors,`
` 162:          /// respectively. DBDSDC can be used to compute all singular values,`
` 163:          /// and optionally, singular vectors or singular vectors in compact form.`
` 164:          /// `
` 165:          /// This code makes very mild assumptions about floating point`
` 166:          /// arithmetic. It will work on machines with a guard digit in`
` 167:          /// add/subtract, or on those binary machines without guard digits`
` 168:          /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.`
` 169:          /// It could conceivably fail on hexadecimal or decimal machines`
` 170:          /// without guard digits, but we know of none.  See DLASD3 for details.`
` 171:          /// `
` 172:          /// The code currently calls DLASDQ if singular values only are desired.`
` 173:          /// However, it can be slightly modified to compute singular values`
` 174:          /// using the divide and conquer method.`
` 175:          /// `
` 176:          ///</summary>`
` 177:          /// <param name="UPLO">`
` 178:          /// (input) CHARACTER*1`
` 179:          /// = 'U':  B is upper bidiagonal.`
` 180:          /// = 'L':  B is lower bidiagonal.`
` 181:          ///</param>`
` 182:          /// <param name="COMPQ">`
` 183:          /// (input) CHARACTER*1`
` 184:          /// Specifies whether singular vectors are to be computed`
` 185:          /// as follows:`
` 186:          /// = 'N':  Compute singular values only;`
` 187:          /// = 'P':  Compute singular values and compute singular`
` 188:          /// vectors in compact form;`
` 189:          /// = 'I':  Compute singular values and singular vectors.`
` 190:          ///</param>`
` 191:          /// <param name="N">`
` 192:          /// (input) INTEGER`
` 193:          /// The order of the matrix B.  N .GE. 0.`
` 194:          ///</param>`
` 195:          /// <param name="D">`
` 196:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 197:          /// On entry, the n diagonal elements of the bidiagonal matrix B.`
` 198:          /// On exit, if INFO=0, the singular values of B.`
` 199:          ///</param>`
` 200:          /// <param name="E">`
` 201:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 202:          /// On entry, the elements of E contain the offdiagonal`
` 203:          /// elements of the bidiagonal matrix whose SVD is desired.`
` 204:          /// On exit, E has been destroyed.`
` 205:          ///</param>`
` 206:          /// <param name="U">`
` 207:          /// (output) DOUBLE PRECISION array, dimension (LDU,N)`
` 208:          /// If  COMPQ = 'I', then:`
` 209:          /// On exit, if INFO = 0, U contains the left singular vectors`
` 210:          /// of the bidiagonal matrix.`
` 211:          /// For other values of COMPQ, U is not referenced.`
` 212:          ///</param>`
` 213:          /// <param name="LDU">`
` 214:          /// (input) INTEGER`
` 215:          /// The leading dimension of the array U.  LDU .GE. 1.`
` 216:          /// If singular vectors are desired, then LDU .GE. max( 1, N ).`
` 217:          ///</param>`
` 218:          /// <param name="VT">`
` 219:          /// (output) DOUBLE PRECISION array, dimension (LDVT,N)`
` 220:          /// If  COMPQ = 'I', then:`
` 221:          /// On exit, if INFO = 0, VT' contains the right singular`
` 222:          /// vectors of the bidiagonal matrix.`
` 223:          /// For other values of COMPQ, VT is not referenced.`
` 224:          ///</param>`
` 225:          /// <param name="LDVT">`
` 226:          /// (input) INTEGER`
` 227:          /// The leading dimension of the array VT.  LDVT .GE. 1.`
` 228:          /// If singular vectors are desired, then LDVT .GE. max( 1, N ).`
` 229:          ///</param>`
` 230:          /// <param name="Q">`
` 231:          /// (output) DOUBLE PRECISION array, dimension (LDQ)`
` 232:          /// If  COMPQ = 'P', then:`
` 233:          /// On exit, if INFO = 0, Q and IQ contain the left`
` 234:          /// and right singular vectors in a compact form,`
` 235:          /// requiring O(N log N) space instead of 2*N**2.`
` 236:          /// In particular, Q contains all the DOUBLE PRECISION data in`
` 237:          /// LDQ .GE. N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))`
` 238:          /// words of memory, where SMLSIZ is returned by ILAENV and`
` 239:          /// is equal to the maximum size of the subproblems at the`
` 240:          /// bottom of the computation tree (usually about 25).`
` 241:          /// For other values of COMPQ, Q is not referenced.`
` 242:          ///</param>`
` 243:          /// <param name="IQ">`
` 244:          /// (output) INTEGER array, dimension (LDIQ)`
` 245:          /// If  COMPQ = 'P', then:`
` 246:          /// On exit, if INFO = 0, Q and IQ contain the left`
` 247:          /// and right singular vectors in a compact form,`
` 248:          /// requiring O(N log N) space instead of 2*N**2.`
` 249:          /// In particular, IQ contains all INTEGER data in`
` 250:          /// LDIQ .GE. N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))`
` 251:          /// words of memory, where SMLSIZ is returned by ILAENV and`
` 252:          /// is equal to the maximum size of the subproblems at the`
` 253:          /// bottom of the computation tree (usually about 25).`
` 254:          /// For other values of COMPQ, IQ is not referenced.`
` 255:          ///</param>`
` 256:          /// <param name="WORK">`
` 257:          /// (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 258:          /// If COMPQ = 'N' then LWORK .GE. (4 * N).`
` 259:          /// If COMPQ = 'P' then LWORK .GE. (6 * N).`
` 260:          /// If COMPQ = 'I' then LWORK .GE. (3 * N**2 + 4 * N).`
` 261:          ///</param>`
` 262:          /// <param name="IWORK">`
` 263:          /// (workspace) INTEGER array, dimension (8*N)`
` 264:          ///</param>`
` 265:          /// <param name="INFO">`
` 266:          /// (output) INTEGER`
` 267:          /// = 0:  successful exit.`
` 268:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 269:          /// .GT. 0:  The algorithm failed to compute an singular value.`
` 270:          /// The update process of divide and conquer failed.`
` 271:          ///</param>`
` 272:          public void Run(string UPLO, string COMPQ, int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] U, int offset_u`
` 273:                           , int LDU, ref double[] VT, int offset_vt, int LDVT, ref double[] Q, int offset_q, ref int[] IQ, int offset_iq, ref double[] WORK, int offset_work`
` 274:                           , ref int[] IWORK, int offset_iwork, ref int INFO)`
` 275:          {`
` 276:   `
` 277:              #region Array Index Correction`
` 278:              `
` 279:               int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_u = -1 - LDU + offset_u;  int o_vt = -1 - LDVT + offset_vt; `
` 280:               int o_q = -1 + offset_q; int o_iq = -1 + offset_iq;  int o_work = -1 + offset_work;  int o_iwork = -1 + offset_iwork; `
` 281:   `
` 282:              #endregion`
` 283:   `
` 284:   `
` 285:              #region Strings`
` 286:              `
` 287:              UPLO = UPLO.Substring(0, 1);  COMPQ = COMPQ.Substring(0, 1);  `
` 288:   `
` 289:              #endregion`
` 290:   `
` 291:   `
` 292:              #region Prolog`
` 293:              `
` 294:              // *`
` 295:              // *  -- LAPACK routine (version 3.1) --`
` 296:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 297:              // *     November 2006`
` 298:              // *`
` 299:              // *     .. Scalar Arguments ..`
` 300:              // *     ..`
` 301:              // *     .. Array Arguments ..`
` 302:              // *     ..`
` 303:              // *`
` 304:              // *  Purpose`
` 305:              // *  =======`
` 306:              // *`
` 307:              // *  DBDSDC computes the singular value decomposition (SVD) of a real`
` 308:              // *  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,`
` 309:              // *  using a divide and conquer method, where S is a diagonal matrix`
` 310:              // *  with non-negative diagonal elements (the singular values of B), and`
` 311:              // *  U and VT are orthogonal matrices of left and right singular vectors,`
` 312:              // *  respectively. DBDSDC can be used to compute all singular values,`
` 313:              // *  and optionally, singular vectors or singular vectors in compact form.`
` 314:              // *`
` 315:              // *  This code makes very mild assumptions about floating point`
` 316:              // *  arithmetic. It will work on machines with a guard digit in`
` 317:              // *  add/subtract, or on those binary machines without guard digits`
` 318:              // *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.`
` 319:              // *  It could conceivably fail on hexadecimal or decimal machines`
` 320:              // *  without guard digits, but we know of none.  See DLASD3 for details.`
` 321:              // *`
` 322:              // *  The code currently calls DLASDQ if singular values only are desired.`
` 323:              // *  However, it can be slightly modified to compute singular values`
` 324:              // *  using the divide and conquer method.`
` 325:              // *`
` 326:              // *  Arguments`
` 327:              // *  =========`
` 328:              // *`
` 329:              // *  UPLO    (input) CHARACTER*1`
` 330:              // *          = 'U':  B is upper bidiagonal.`
` 331:              // *          = 'L':  B is lower bidiagonal.`
` 332:              // *`
` 333:              // *  COMPQ   (input) CHARACTER*1`
` 334:              // *          Specifies whether singular vectors are to be computed`
` 335:              // *          as follows:`
` 336:              // *          = 'N':  Compute singular values only;`
` 337:              // *          = 'P':  Compute singular values and compute singular`
` 338:              // *                  vectors in compact form;`
` 339:              // *          = 'I':  Compute singular values and singular vectors.`
` 340:              // *`
` 341:              // *  N       (input) INTEGER`
` 342:              // *          The order of the matrix B.  N >= 0.`
` 343:              // *`
` 344:              // *  D       (input/output) DOUBLE PRECISION array, dimension (N)`
` 345:              // *          On entry, the n diagonal elements of the bidiagonal matrix B.`
` 346:              // *          On exit, if INFO=0, the singular values of B.`
` 347:              // *`
` 348:              // *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 349:              // *          On entry, the elements of E contain the offdiagonal`
` 350:              // *          elements of the bidiagonal matrix whose SVD is desired.`
` 351:              // *          On exit, E has been destroyed.`
` 352:              // *`
` 353:              // *  U       (output) DOUBLE PRECISION array, dimension (LDU,N)`
` 354:              // *          If  COMPQ = 'I', then:`
` 355:              // *             On exit, if INFO = 0, U contains the left singular vectors`
` 356:              // *             of the bidiagonal matrix.`
` 357:              // *          For other values of COMPQ, U is not referenced.`
` 358:              // *`
` 359:              // *  LDU     (input) INTEGER`
` 360:              // *          The leading dimension of the array U.  LDU >= 1.`
` 361:              // *          If singular vectors are desired, then LDU >= max( 1, N ).`
` 362:              // *`
` 363:              // *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)`
` 364:              // *          If  COMPQ = 'I', then:`
` 365:              // *             On exit, if INFO = 0, VT' contains the right singular`
` 366:              // *             vectors of the bidiagonal matrix.`
` 367:              // *          For other values of COMPQ, VT is not referenced.`
` 368:              // *`
` 369:              // *  LDVT    (input) INTEGER`
` 370:              // *          The leading dimension of the array VT.  LDVT >= 1.`
` 371:              // *          If singular vectors are desired, then LDVT >= max( 1, N ).`
` 372:              // *`
` 373:              // *  Q       (output) DOUBLE PRECISION array, dimension (LDQ)`
` 374:              // *          If  COMPQ = 'P', then:`
` 375:              // *             On exit, if INFO = 0, Q and IQ contain the left`
` 376:              // *             and right singular vectors in a compact form,`
` 377:              // *             requiring O(N log N) space instead of 2*N**2.`
` 378:              // *             In particular, Q contains all the DOUBLE PRECISION data in`
` 379:              // *             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))`
` 380:              // *             words of memory, where SMLSIZ is returned by ILAENV and`
` 381:              // *             is equal to the maximum size of the subproblems at the`
` 382:              // *             bottom of the computation tree (usually about 25).`
` 383:              // *          For other values of COMPQ, Q is not referenced.`
` 384:              // *`
` 385:              // *  IQ      (output) INTEGER array, dimension (LDIQ)`
` 386:              // *          If  COMPQ = 'P', then:`
` 387:              // *             On exit, if INFO = 0, Q and IQ contain the left`
` 388:              // *             and right singular vectors in a compact form,`
` 389:              // *             requiring O(N log N) space instead of 2*N**2.`
` 390:              // *             In particular, IQ contains all INTEGER data in`
` 391:              // *             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))`
` 392:              // *             words of memory, where SMLSIZ is returned by ILAENV and`
` 393:              // *             is equal to the maximum size of the subproblems at the`
` 394:              // *             bottom of the computation tree (usually about 25).`
` 395:              // *          For other values of COMPQ, IQ is not referenced.`
` 396:              // *`
` 397:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))`
` 398:              // *          If COMPQ = 'N' then LWORK >= (4 * N).`
` 399:              // *          If COMPQ = 'P' then LWORK >= (6 * N).`
` 400:              // *          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).`
` 401:              // *`
` 402:              // *  IWORK   (workspace) INTEGER array, dimension (8*N)`
` 403:              // *`
` 404:              // *  INFO    (output) INTEGER`
` 405:              // *          = 0:  successful exit.`
` 406:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 407:              // *          > 0:  The algorithm failed to compute an singular value.`
` 408:              // *                The update process of divide and conquer failed.`
` 409:              // *`
` 410:              // *  Further Details`
` 411:              // *  ===============`
` 412:              // *`
` 413:              // *  Based on contributions by`
` 414:              // *     Ming Gu and Huan Ren, Computer Science Division, University of`
` 415:              // *     California at Berkeley, USA`
` 416:              // *`
` 417:              // *  =====================================================================`
` 418:              // *  Changed dimension statement in comment describing E from (N) to`
` 419:              // *  (N-1).  Sven, 17 Feb 05.`
` 420:              // *  =====================================================================`
` 421:              // *`
` 422:              // *     .. Parameters ..`
` 423:              // *     ..`
` 424:              // *     .. Local Scalars ..`
` 425:              // *     ..`
` 426:              // *     .. External Functions ..`
` 427:              // *     ..`
` 428:              // *     .. External Subroutines ..`
` 429:              // *     ..`
` 430:              // *     .. Intrinsic Functions ..`
` 431:              //      INTRINSIC          ABS, DBLE, INT, LOG, SIGN;`
` 432:              // *     ..`
` 433:              // *     .. Executable Statements ..`
` 434:              // *`
` 435:              // *     Test the input parameters.`
` 436:              // *`
` 437:   `
` 438:              #endregion`
` 439:   `
` 440:   `
` 441:              #region Body`
` 442:              `
` 443:              INFO = 0;`
` 444:              // *`
` 445:              IUPLO = 0;`
` 446:              if (this._lsame.Run(UPLO, "U")) IUPLO = 1;`
` 447:              if (this._lsame.Run(UPLO, "L")) IUPLO = 2;`
` 448:              if (this._lsame.Run(COMPQ, "N"))`
` 449:              {`
` 450:                  ICOMPQ = 0;`
` 451:              }`
` 452:              else`
` 453:              {`
` 454:                  if (this._lsame.Run(COMPQ, "P"))`
` 455:                  {`
` 456:                      ICOMPQ = 1;`
` 457:                  }`
` 458:                  else`
` 459:                  {`
` 460:                      if (this._lsame.Run(COMPQ, "I"))`
` 461:                      {`
` 462:                          ICOMPQ = 2;`
` 463:                      }`
` 464:                      else`
` 465:                      {`
` 466:                          ICOMPQ =  - 1;`
` 467:                      }`
` 468:                  }`
` 469:              }`
` 470:              if (IUPLO == 0)`
` 471:              {`
` 472:                  INFO =  - 1;`
` 473:              }`
` 474:              else`
` 475:              {`
` 476:                  if (ICOMPQ < 0)`
` 477:                  {`
` 478:                      INFO =  - 2;`
` 479:                  }`
` 480:                  else`
` 481:                  {`
` 482:                      if (N < 0)`
` 483:                      {`
` 484:                          INFO =  - 3;`
` 485:                      }`
` 486:                      else`
` 487:                      {`
` 488:                          if ((LDU < 1) || ((ICOMPQ == 2) && (LDU < N)))`
` 489:                          {`
` 490:                              INFO =  - 7;`
` 491:                          }`
` 492:                          else`
` 493:                          {`
` 494:                              if ((LDVT < 1) || ((ICOMPQ == 2) && (LDVT < N)))`
` 495:                              {`
` 496:                                  INFO =  - 9;`
` 497:                              }`
` 498:                          }`
` 499:                      }`
` 500:                  }`
` 501:              }`
` 502:              if (INFO != 0)`
` 503:              {`
` 504:                  this._xerbla.Run("DBDSDC",  - INFO);`
` 505:                  return;`
` 506:              }`
` 507:              // *`
` 508:              // *     Quick return if possible`
` 509:              // *`
` 510:              if (N == 0) return;`
` 511:              SMLSIZ = this._ilaenv.Run(9, "DBDSDC", " ", 0, 0, 0, 0);`
` 512:              if (N == 1)`
` 513:              {`
` 514:                  if (ICOMPQ == 1)`
` 515:                  {`
` 516:                      Q[1 + o_q] = FortranLib.Sign(ONE,D[1 + o_d]);`
` 517:                      Q[1 + SMLSIZ * N + o_q] = ONE;`
` 518:                  }`
` 519:                  else`
` 520:                  {`
` 521:                      if (ICOMPQ == 2)`
` 522:                      {`
` 523:                          U[1+1 * LDU + o_u] = FortranLib.Sign(ONE,D[1 + o_d]);`
` 524:                          VT[1+1 * LDVT + o_vt] = ONE;`
` 525:                      }`
` 526:                  }`
` 527:                  D[1 + o_d] = Math.Abs(D[1 + o_d]);`
` 528:                  return;`
` 529:              }`
` 530:              NM1 = N - 1;`
` 531:              // *`
` 532:              // *     If matrix lower bidiagonal, rotate to be upper bidiagonal`
` 533:              // *     by applying Givens rotations on the left`
` 534:              // *`
` 535:              WSTART = 1;`
` 536:              QSTART = 3;`
` 537:              if (ICOMPQ == 1)`
` 538:              {`
` 539:                  this._dcopy.Run(N, D, offset_d, 1, ref Q, 1 + o_q, 1);`
` 540:                  this._dcopy.Run(N - 1, E, offset_e, 1, ref Q, N + 1 + o_q, 1);`
` 541:              }`
` 542:              if (IUPLO == 2)`
` 543:              {`
` 544:                  QSTART = 5;`
` 545:                  WSTART = 2 * N - 1;`
` 546:                  for (I = 1; I <= N - 1; I++)`
` 547:                  {`
` 548:                      this._dlartg.Run(D[I + o_d], E[I + o_e], ref CS, ref SN, ref R);`
` 549:                      D[I + o_d] = R;`
` 550:                      E[I + o_e] = SN * D[I + 1 + o_d];`
` 551:                      D[I + 1 + o_d] = CS * D[I + 1 + o_d];`
` 552:                      if (ICOMPQ == 1)`
` 553:                      {`
` 554:                          Q[I + 2 * N + o_q] = CS;`
` 555:                          Q[I + 3 * N + o_q] = SN;`
` 556:                      }`
` 557:                      else`
` 558:                      {`
` 559:                          if (ICOMPQ == 2)`
` 560:                          {`
` 561:                              WORK[I + o_work] = CS;`
` 562:                              WORK[NM1 + I + o_work] =  - SN;`
` 563:                          }`
` 564:                      }`
` 565:                  }`
` 566:              }`
` 567:              // *`
` 568:              // *     If ICOMPQ = 0, use DLASDQ to compute the singular values.`
` 569:              // *`
` 570:              if (ICOMPQ == 0)`
` 571:              {`
` 572:                  this._dlasdq.Run("U", 0, N, 0, 0, 0`
` 573:                                   , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDVT, ref U, offset_u, LDU`
` 574:                                   , ref U, offset_u, LDU, ref WORK, WSTART + o_work, ref INFO);`
` 575:                  goto LABEL40;`
` 576:              }`
` 577:              // *`
` 578:              // *     If N is smaller than the minimum divide size SMLSIZ, then solve`
` 579:              // *     the problem with another solver.`
` 580:              // *`
` 581:              if (N <= SMLSIZ)`
` 582:              {`
` 583:                  if (ICOMPQ == 2)`
` 584:                  {`
` 585:                      this._dlaset.Run("A", N, N, ZERO, ONE, ref U, offset_u`
` 586:                                       , LDU);`
` 587:                      this._dlaset.Run("A", N, N, ZERO, ONE, ref VT, offset_vt`
` 588:                                       , LDVT);`
` 589:                      this._dlasdq.Run("U", 0, N, N, N, 0`
` 590:                                       , ref D, offset_d, ref E, offset_e, ref VT, offset_vt, LDVT, ref U, offset_u, LDU`
` 591:                                       , ref U, offset_u, LDU, ref WORK, WSTART + o_work, ref INFO);`
` 592:                  }`
` 593:                  else`
` 594:                  {`
` 595:                      if (ICOMPQ == 1)`
` 596:                      {`
` 597:                          IU = 1;`
` 598:                          IVT = IU + N;`
` 599:                          this._dlaset.Run("A", N, N, ZERO, ONE, ref Q, IU + (QSTART - 1) * N + o_q`
` 600:                                           , N);`
` 601:                          this._dlaset.Run("A", N, N, ZERO, ONE, ref Q, IVT + (QSTART - 1) * N + o_q`
` 602:                                           , N);`
` 603:                          this._dlasdq.Run("U", 0, N, N, N, 0`
` 604:                                           , ref D, offset_d, ref E, offset_e, ref Q, IVT + (QSTART - 1) * N + o_q, N, ref Q, IU + (QSTART - 1) * N + o_q, N`
` 605:                                           , ref Q, IU + (QSTART - 1) * N + o_q, N, ref WORK, WSTART + o_work, ref INFO);`
` 606:                      }`
` 607:                  }`
` 608:                  goto LABEL40;`
` 609:              }`
` 610:              // *`
` 611:              if (ICOMPQ == 2)`
` 612:              {`
` 613:                  this._dlaset.Run("A", N, N, ZERO, ONE, ref U, offset_u`
` 614:                                   , LDU);`
` 615:                  this._dlaset.Run("A", N, N, ZERO, ONE, ref VT, offset_vt`
` 616:                                   , LDVT);`
` 617:              }`
` 618:              // *`
` 619:              // *     Scale.`
` 620:              // *`
` 621:              ORGNRM = this._dlanst.Run("M", N, D, offset_d, E, offset_e);`
` 622:              if (ORGNRM == ZERO) return;`
` 623:              this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N`
` 624:                               , 1, ref D, offset_d, N, ref IERR);`
` 625:              this._dlascl.Run("G", 0, 0, ORGNRM, ONE, NM1`
` 626:                               , 1, ref E, offset_e, NM1, ref IERR);`
` 627:              // *`
` 628:              EPS = this._dlamch.Run("Epsilon");`
` 629:              // *`
` 630:              MLVL = Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(N) / Convert.ToDouble(SMLSIZ + 1)) / Math.Log(TWO))) + 1;`
` 631:              SMLSZP = SMLSIZ + 1;`
` 632:              // *`
` 633:              if (ICOMPQ == 1)`
` 634:              {`
` 635:                  IU = 1;`
` 636:                  IVT = 1 + SMLSIZ;`
` 637:                  DIFL = IVT + SMLSZP;`
` 638:                  DIFR = DIFL + MLVL;`
` 639:                  Z = DIFR + MLVL * 2;`
` 640:                  IC = Z + MLVL;`
` 641:                  IS = IC + 1;`
` 642:                  POLES = IS + 1;`
` 643:                  GIVNUM = POLES + 2 * MLVL;`
` 644:                  // *`
` 645:                  K = 1;`
` 646:                  GIVPTR = 2;`
` 647:                  PERM = 3;`
` 648:                  GIVCOL = PERM + MLVL;`
` 649:              }`
` 650:              // *`
` 651:              for (I = 1; I <= N; I++)`
` 652:              {`
` 653:                  if (Math.Abs(D[I + o_d]) < EPS)`
` 654:                  {`
` 655:                      D[I + o_d] = FortranLib.Sign(EPS,D[I + o_d]);`
` 656:                  }`
` 657:              }`
` 658:              // *`
` 659:              START = 1;`
` 660:              SQRE = 0;`
` 661:              // *`
` 662:              for (I = 1; I <= NM1; I++)`
` 663:              {`
` 664:                  if ((Math.Abs(E[I + o_e]) < EPS) || (I == NM1))`
` 665:                  {`
` 666:                      // *`
` 667:                      // *        Subproblem found. First determine its size and then`
` 668:                      // *        apply divide and conquer on it.`
` 669:                      // *`
` 670:                      if (I < NM1)`
` 671:                      {`
` 672:                          // *`
` 673:                          // *        A subproblem with E(I) small for I < NM1.`
` 674:                          // *`
` 675:                          NSIZE = I - START + 1;`
` 676:                      }`
` 677:                      else`
` 678:                      {`
` 679:                          if (Math.Abs(E[I + o_e]) >= EPS)`
` 680:                          {`
` 681:                              // *`
` 682:                              // *        A subproblem with E(NM1) not too small but I = NM1.`
` 683:                              // *`
` 684:                              NSIZE = N - START + 1;`
` 685:                          }`
` 686:                          else`
` 687:                          {`
` 688:                              // *`
` 689:                              // *        A subproblem with E(NM1) small. This implies an`
` 690:                              // *        1-by-1 subproblem at D(N). Solve this 1-by-1 problem`
` 691:                              // *        first.`
` 692:                              // *`
` 693:                              NSIZE = I - START + 1;`
` 694:                              if (ICOMPQ == 2)`
` 695:                              {`
` 696:                                  U[N+N * LDU + o_u] = FortranLib.Sign(ONE,D[N + o_d]);`
` 697:                                  VT[N+N * LDVT + o_vt] = ONE;`
` 698:                              }`
` 699:                              else`
` 700:                              {`
` 701:                                  if (ICOMPQ == 1)`
` 702:                                  {`
` 703:                                      Q[N + (QSTART - 1) * N + o_q] = FortranLib.Sign(ONE,D[N + o_d]);`
` 704:                                      Q[N + (SMLSIZ + QSTART - 1) * N + o_q] = ONE;`
` 705:                                  }`
` 706:                              }`
` 707:                              D[N + o_d] = Math.Abs(D[N + o_d]);`
` 708:                          }`
` 709:                      }`
` 710:                      if (ICOMPQ == 2)`
` 711:                      {`
` 712:                          this._dlasd0.Run(NSIZE, SQRE, ref D, START + o_d, ref E, START + o_e, ref U, START+START * LDU + o_u, LDU`
` 713:                                           , ref VT, START+START * LDVT + o_vt, LDVT, SMLSIZ, ref IWORK, offset_iwork, ref WORK, WSTART + o_work, ref INFO);`
` 714:                      }`
` 715:                      else`
` 716:                      {`
` 717:                          this._dlasda.Run(ICOMPQ, SMLSIZ, NSIZE, SQRE, ref D, START + o_d, ref E, START + o_e`
` 718:                                           , ref Q, START + (IU + QSTART - 2) * N + o_q, N, ref Q, START + (IVT + QSTART - 2) * N + o_q, ref IQ, START + K * N + o_iq, ref Q, START + (DIFL + QSTART - 2) * N + o_q, ref Q, START + (DIFR + QSTART - 2) * N + o_q`
` 719:                                           , ref Q, START + (Z + QSTART - 2) * N + o_q, ref Q, START + (POLES + QSTART - 2) * N + o_q, ref IQ, START + GIVPTR * N + o_iq, ref IQ, START + GIVCOL * N + o_iq, N, ref IQ, START + PERM * N + o_iq`
` 720:                                           , ref Q, START + (GIVNUM + QSTART - 2) * N + o_q, ref Q, START + (IC + QSTART - 2) * N + o_q, ref Q, START + (IS + QSTART - 2) * N + o_q, ref WORK, WSTART + o_work, ref IWORK, offset_iwork, ref INFO);`
` 721:                          if (INFO != 0)`
` 722:                          {`
` 723:                              return;`
` 724:                          }`
` 725:                      }`
` 726:                      START = I + 1;`
` 727:                  }`
` 728:              }`
` 729:              // *`
` 730:              // *     Unscale`
` 731:              // *`
` 732:              this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N`
` 733:                               , 1, ref D, offset_d, N, ref IERR);`
` 734:          LABEL40:;`
` 735:              // *`
` 736:              // *     Use Selection Sort to minimize swaps of singular vectors`
` 737:              // *`
` 738:              for (II = 2; II <= N; II++)`
` 739:              {`
` 740:                  I = II - 1;`
` 741:                  KK = I;`
` 742:                  P = D[I + o_d];`
` 743:                  for (J = II; J <= N; J++)`
` 744:                  {`
` 745:                      if (D[J + o_d] > P)`
` 746:                      {`
` 747:                          KK = J;`
` 748:                          P = D[J + o_d];`
` 749:                      }`
` 750:                  }`
` 751:                  if (KK != I)`
` 752:                  {`
` 753:                      D[KK + o_d] = D[I + o_d];`
` 754:                      D[I + o_d] = P;`
` 755:                      if (ICOMPQ == 1)`
` 756:                      {`
` 757:                          IQ[I + o_iq] = KK;`
` 758:                      }`
` 759:                      else`
` 760:                      {`
` 761:                          if (ICOMPQ == 2)`
` 762:                          {`
` 763:                              this._dswap.Run(N, ref U, 1+I * LDU + o_u, 1, ref U, 1+KK * LDU + o_u, 1);`
` 764:                              this._dswap.Run(N, ref VT, I+1 * LDVT + o_vt, LDVT, ref VT, KK+1 * LDVT + o_vt, LDVT);`
` 765:                          }`
` 766:                      }`
` 767:                  }`
` 768:                  else`
` 769:                  {`
` 770:                      if (ICOMPQ == 1)`
` 771:                      {`
` 772:                          IQ[I + o_iq] = I;`
` 773:                      }`
` 774:                  }`
` 775:              }`
` 776:              // *`
` 777:              // *     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO`
` 778:              // *`
` 779:              if (ICOMPQ == 1)`
` 780:              {`
` 781:                  if (IUPLO == 1)`
` 782:                  {`
` 783:                      IQ[N + o_iq] = 1;`
` 784:                  }`
` 785:                  else`
` 786:                  {`
` 787:                      IQ[N + o_iq] = 0;`
` 788:                  }`
` 789:              }`
` 790:              // *`
` 791:              // *     If B is lower bidiagonal, update U by those Givens rotations`
` 792:              // *     which rotated B to be upper bidiagonal`
` 793:              // *`
` 794:              if ((IUPLO == 2) && (ICOMPQ == 2))`
` 795:              {`
` 796:                  this._dlasr.Run("L", "V", "B", N, N, WORK, 1 + o_work`
` 797:                                  , WORK, N + o_work, ref U, offset_u, LDU);`
` 798:              }`
` 799:              // *`
` 800:              return;`
` 801:              // *`
` 802:              // *     End of DBDSDC`
` 803:              // *`
` 804:   `
` 805:              #endregion`
` 806:   `
` 807:          }`
` 808:      }`
` 809:  }`