`   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.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// January 2007`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DBDSQR computes the singular values and, optionally, the right and/or`
`  27:      /// left singular vectors from the singular value decomposition (SVD) of`
`  28:      /// a real N-by-N (upper or lower) bidiagonal matrix B using the implicit`
`  29:      /// zero-shift QR algorithm.  The SVD of B has the form`
`  30:      /// `
`  31:      /// B = Q * S * P**T`
`  32:      /// `
`  33:      /// where S is the diagonal matrix of singular values, Q is an orthogonal`
`  34:      /// matrix of left singular vectors, and P is an orthogonal matrix of`
`  35:      /// right singular vectors.  If left singular vectors are requested, this`
`  36:      /// subroutine actually returns U*Q instead of Q, and, if right singular`
`  37:      /// vectors are requested, this subroutine returns P**T*VT instead of`
`  38:      /// P**T, for given real input matrices U and VT.  When U and VT are the`
`  39:      /// orthogonal matrices that reduce a general matrix A to bidiagonal`
`  40:      /// form:  A = U*B*VT, as computed by DGEBRD, then`
`  41:      /// `
`  42:      /// A = (U*Q) * S * (P**T*VT)`
`  43:      /// `
`  44:      /// is the SVD of A.  Optionally, the subroutine may also compute Q**T*C`
`  45:      /// for a given real input matrix C.`
`  46:      /// `
`  47:      /// See "Computing  Small Singular Values of Bidiagonal Matrices With`
`  48:      /// Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,`
`  49:      /// LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,`
`  50:      /// no. 5, pp. 873-912, Sept 1990) and`
`  51:      /// "Accurate singular values and differential qd algorithms," by`
`  52:      /// B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics`
`  53:      /// Department, University of California at Berkeley, July 1992`
`  54:      /// for a detailed description of the algorithm.`
`  55:      /// `
`  56:      ///</summary>`
`  57:      public class DBDSQR`
`  58:      {`
`  59:      `
`  60:   `
`  61:          #region Dependencies`
`  62:          `
`  63:          LSAME _lsame; DLAMCH _dlamch; DLARTG _dlartg; DLAS2 _dlas2; DLASQ1 _dlasq1; DLASR _dlasr; DLASV2 _dlasv2; DROT _drot; `
`  64:          DSCAL _dscal;DSWAP _dswap; XERBLA _xerbla; `
`  65:   `
`  66:          #endregion`
`  67:   `
`  68:   `
`  69:          #region Fields`
`  70:          `
`  71:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double NEGONE =  - 1.0E0; const double HNDRTH = 0.01E0; `
`  72:          const double TEN = 10.0E0;const double HNDRD = 100.0E0; const double MEIGTH =  - 0.125E0; const int MAXITR = 6; `
`  73:          bool LOWER = false;bool ROTATE = false; int I = 0; int IDIR = 0; int ISUB = 0; int ITER = 0; int J = 0; int LL = 0; `
`  74:          int LLL = 0;int M = 0; int MAXIT = 0; int NM1 = 0; int NM12 = 0; int NM13 = 0; int OLDLL = 0; int OLDM = 0; `
`  75:          double ABSE = 0;double ABSS = 0; double COSL = 0; double COSR = 0; double CS = 0; double EPS = 0; double F = 0; `
`  76:          double G = 0;double H = 0; double MU = 0; double OLDCS = 0; double OLDSN = 0; double R = 0; double SHIFT = 0; `
`  77:          double SIGMN = 0;double SIGMX = 0; double SINL = 0; double SINR = 0; double SLL = 0; double SMAX = 0; double SMIN = 0; `
`  78:          double SMINL = 0;double SMINOA = 0; double SN = 0; double THRESH = 0; double TOL = 0; double TOLMUL = 0; double UNFL = 0; `
`  79:   `
`  80:          #endregion`
`  81:   `
`  82:          public DBDSQR(LSAME lsame, DLAMCH dlamch, DLARTG dlartg, DLAS2 dlas2, DLASQ1 dlasq1, DLASR dlasr, DLASV2 dlasv2, DROT drot, DSCAL dscal, DSWAP dswap`
`  83:                        , XERBLA xerbla)`
`  84:          {`
`  85:      `
`  86:   `
`  87:              #region Set Dependencies`
`  88:              `
`  89:              this._lsame = lsame; this._dlamch = dlamch; this._dlartg = dlartg; this._dlas2 = dlas2; this._dlasq1 = dlasq1; `
`  90:              this._dlasr = dlasr;this._dlasv2 = dlasv2; this._drot = drot; this._dscal = dscal; this._dswap = dswap; `
`  91:              this._xerbla = xerbla;`
`  92:   `
`  93:              #endregion`
`  94:   `
`  95:          }`
`  96:      `
`  97:          public DBDSQR()`
`  98:          {`
`  99:      `
` 100:   `
` 101:              #region Dependencies (Initialization)`
` 102:              `
` 103:              LSAME lsame = new LSAME();`
` 104:              DLAMC3 dlamc3 = new DLAMC3();`
` 105:              DLAS2 dlas2 = new DLAS2();`
` 106:              DCOPY dcopy = new DCOPY();`
` 107:              XERBLA xerbla = new XERBLA();`
` 108:              DLASQ5 dlasq5 = new DLASQ5();`
` 109:              DLAZQ4 dlazq4 = new DLAZQ4();`
` 110:              IEEECK ieeeck = new IEEECK();`
` 111:              IPARMQ iparmq = new IPARMQ();`
` 112:              DROT drot = new DROT();`
` 113:              DSCAL dscal = new DSCAL();`
` 114:              DSWAP dswap = new DSWAP();`
` 115:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 116:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 117:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 118:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 119:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 120:              DLARTG dlartg = new DLARTG(dlamch);`
` 121:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);`
` 122:              DLASQ6 dlasq6 = new DLASQ6(dlamch);`
` 123:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);`
` 124:              DLASRT dlasrt = new DLASRT(lsame, xerbla);`
` 125:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 126:              DLASQ2 dlasq2 = new DLASQ2(dlazq3, dlasrt, xerbla, dlamch, ilaenv);`
` 127:              DLASQ1 dlasq1 = new DLASQ1(dcopy, dlas2, dlascl, dlasq2, dlasrt, xerbla, dlamch);`
` 128:              DLASR dlasr = new DLASR(lsame, xerbla);`
` 129:              DLASV2 dlasv2 = new DLASV2(dlamch);`
` 130:   `
` 131:              #endregion`
` 132:   `
` 133:   `
` 134:              #region Set Dependencies`
` 135:              `
` 136:              this._lsame = lsame; this._dlamch = dlamch; this._dlartg = dlartg; this._dlas2 = dlas2; this._dlasq1 = dlasq1; `
` 137:              this._dlasr = dlasr;this._dlasv2 = dlasv2; this._drot = drot; this._dscal = dscal; this._dswap = dswap; `
` 138:              this._xerbla = xerbla;`
` 139:   `
` 140:              #endregion`
` 141:   `
` 142:          }`
` 143:          /// <summary>`
` 144:          /// Purpose`
` 145:          /// =======`
` 146:          /// `
` 147:          /// DBDSQR computes the singular values and, optionally, the right and/or`
` 148:          /// left singular vectors from the singular value decomposition (SVD) of`
` 149:          /// a real N-by-N (upper or lower) bidiagonal matrix B using the implicit`
` 150:          /// zero-shift QR algorithm.  The SVD of B has the form`
` 151:          /// `
` 152:          /// B = Q * S * P**T`
` 153:          /// `
` 154:          /// where S is the diagonal matrix of singular values, Q is an orthogonal`
` 155:          /// matrix of left singular vectors, and P is an orthogonal matrix of`
` 156:          /// right singular vectors.  If left singular vectors are requested, this`
` 157:          /// subroutine actually returns U*Q instead of Q, and, if right singular`
` 158:          /// vectors are requested, this subroutine returns P**T*VT instead of`
` 159:          /// P**T, for given real input matrices U and VT.  When U and VT are the`
` 160:          /// orthogonal matrices that reduce a general matrix A to bidiagonal`
` 161:          /// form:  A = U*B*VT, as computed by DGEBRD, then`
` 162:          /// `
` 163:          /// A = (U*Q) * S * (P**T*VT)`
` 164:          /// `
` 165:          /// is the SVD of A.  Optionally, the subroutine may also compute Q**T*C`
` 166:          /// for a given real input matrix C.`
` 167:          /// `
` 168:          /// See "Computing  Small Singular Values of Bidiagonal Matrices With`
` 169:          /// Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,`
` 170:          /// LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,`
` 171:          /// no. 5, pp. 873-912, Sept 1990) and`
` 172:          /// "Accurate singular values and differential qd algorithms," by`
` 173:          /// B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics`
` 174:          /// Department, University of California at Berkeley, July 1992`
` 175:          /// for a detailed description of the algorithm.`
` 176:          /// `
` 177:          ///</summary>`
` 178:          /// <param name="UPLO">`
` 179:          /// (input) CHARACTER*1`
` 180:          /// = 'U':  B is upper bidiagonal;`
` 181:          /// = 'L':  B is lower bidiagonal.`
` 182:          ///</param>`
` 183:          /// <param name="N">`
` 184:          /// (input) INTEGER`
` 185:          /// The order of the matrix B.  N .GE. 0.`
` 186:          ///</param>`
` 187:          /// <param name="NCVT">`
` 188:          /// (input) INTEGER`
` 189:          /// The number of columns of the matrix VT. NCVT .GE. 0.`
` 190:          ///</param>`
` 191:          /// <param name="NRU">`
` 192:          /// (input) INTEGER`
` 193:          /// The number of rows of the matrix U. NRU .GE. 0.`
` 194:          ///</param>`
` 195:          /// <param name="NCC">`
` 196:          /// (input) INTEGER`
` 197:          /// The number of columns of the matrix C. NCC .GE. 0.`
` 198:          ///</param>`
` 199:          /// <param name="D">`
` 200:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 201:          /// On entry, the n diagonal elements of the bidiagonal matrix B.`
` 202:          /// On exit, if INFO=0, the singular values of B in decreasing`
` 203:          /// order.`
` 204:          ///</param>`
` 205:          /// <param name="E">`
` 206:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 207:          /// On entry, the N-1 offdiagonal elements of the bidiagonal`
` 208:          /// matrix B. `
` 209:          /// On exit, if INFO = 0, E is destroyed; if INFO .GT. 0, D and E`
` 210:          /// will contain the diagonal and superdiagonal elements of a`
` 211:          /// bidiagonal matrix orthogonally equivalent to the one given`
` 212:          /// as input.`
` 213:          ///</param>`
` 214:          /// <param name="VT">`
` 215:          /// (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)`
` 216:          /// On entry, an N-by-NCVT matrix VT.`
` 217:          /// On exit, VT is overwritten by P**T * VT.`
` 218:          /// Not referenced if NCVT = 0.`
` 219:          ///</param>`
` 220:          /// <param name="LDVT">`
` 221:          /// (input) INTEGER`
` 222:          /// The leading dimension of the array VT.`
` 223:          /// LDVT .GE. max(1,N) if NCVT .GT. 0; LDVT .GE. 1 if NCVT = 0.`
` 224:          ///</param>`
` 225:          /// <param name="U">`
` 226:          /// (input/output) DOUBLE PRECISION array, dimension (LDU, N)`
` 227:          /// On entry, an NRU-by-N matrix U.`
` 228:          /// On exit, U is overwritten by U * Q.`
` 229:          /// Not referenced if NRU = 0.`
` 230:          ///</param>`
` 231:          /// <param name="LDU">`
` 232:          /// (input) INTEGER`
` 233:          /// The leading dimension of the array U.  LDU .GE. max(1,NRU).`
` 234:          ///</param>`
` 235:          /// <param name="C">`
` 236:          /// (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)`
` 237:          /// On entry, an N-by-NCC matrix C.`
` 238:          /// On exit, C is overwritten by Q**T * C.`
` 239:          /// Not referenced if NCC = 0.`
` 240:          ///</param>`
` 241:          /// <param name="LDC">`
` 242:          /// (input) INTEGER`
` 243:          /// The leading dimension of the array C.`
` 244:          /// LDC .GE. max(1,N) if NCC .GT. 0; LDC .GE.1 if NCC = 0.`
` 245:          ///</param>`
` 246:          /// <param name="WORK">`
` 247:          /// (workspace) DOUBLE PRECISION array, dimension (2*N)`
` 248:          /// if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise`
` 249:          ///</param>`
` 250:          /// <param name="INFO">`
` 251:          /// (output) INTEGER`
` 252:          /// = 0:  successful exit`
` 253:          /// .LT. 0:  If INFO = -i, the i-th argument had an illegal value`
` 254:          /// .GT. 0:  the algorithm did not converge; D and E contain the`
` 255:          /// elements of a bidiagonal matrix which is orthogonally`
` 256:          /// similar to the input matrix B;  if INFO = i, i`
` 257:          /// elements of E have not converged to zero.`
` 258:          ///</param>`
` 259:          public void Run(string UPLO, int N, int NCVT, int NRU, int NCC, ref double[] D, int offset_d`
` 260:                           , ref double[] E, int offset_e, ref double[] VT, int offset_vt, int LDVT, ref double[] U, int offset_u, int LDU, ref double[] C, int offset_c`
` 261:                           , int LDC, ref double[] WORK, int offset_work, ref int INFO)`
` 262:          {`
` 263:   `
` 264:              #region Array Index Correction`
` 265:              `
` 266:               int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_vt = -1 - LDVT + offset_vt;  int o_u = -1 - LDU + offset_u; `
` 267:               int o_c = -1 - LDC + offset_c; int o_work = -1 + offset_work; `
` 268:   `
` 269:              #endregion`
` 270:   `
` 271:   `
` 272:              #region Strings`
` 273:              `
` 274:              UPLO = UPLO.Substring(0, 1);  `
` 275:   `
` 276:              #endregion`
` 277:   `
` 278:   `
` 279:              #region Prolog`
` 280:              `
` 281:              // *`
` 282:              // *  -- LAPACK routine (version 3.1.1) --`
` 283:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 284:              // *     January 2007`
` 285:              // *`
` 286:              // *     .. Scalar Arguments ..`
` 287:              // *     ..`
` 288:              // *     .. Array Arguments ..`
` 289:              // *     ..`
` 290:              // *`
` 291:              // *  Purpose`
` 292:              // *  =======`
` 293:              // *`
` 294:              // *  DBDSQR computes the singular values and, optionally, the right and/or`
` 295:              // *  left singular vectors from the singular value decomposition (SVD) of`
` 296:              // *  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit`
` 297:              // *  zero-shift QR algorithm.  The SVD of B has the form`
` 298:              // * `
` 299:              // *     B = Q * S * P**T`
` 300:              // * `
` 301:              // *  where S is the diagonal matrix of singular values, Q is an orthogonal`
` 302:              // *  matrix of left singular vectors, and P is an orthogonal matrix of`
` 303:              // *  right singular vectors.  If left singular vectors are requested, this`
` 304:              // *  subroutine actually returns U*Q instead of Q, and, if right singular`
` 305:              // *  vectors are requested, this subroutine returns P**T*VT instead of`
` 306:              // *  P**T, for given real input matrices U and VT.  When U and VT are the`
` 307:              // *  orthogonal matrices that reduce a general matrix A to bidiagonal`
` 308:              // *  form:  A = U*B*VT, as computed by DGEBRD, then`
` 309:              // *`
` 310:              // *     A = (U*Q) * S * (P**T*VT)`
` 311:              // *`
` 312:              // *  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C`
` 313:              // *  for a given real input matrix C.`
` 314:              // *`
` 315:              // *  See "Computing  Small Singular Values of Bidiagonal Matrices With`
` 316:              // *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,`
` 317:              // *  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,`
` 318:              // *  no. 5, pp. 873-912, Sept 1990) and`
` 319:              // *  "Accurate singular values and differential qd algorithms," by`
` 320:              // *  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics`
` 321:              // *  Department, University of California at Berkeley, July 1992`
` 322:              // *  for a detailed description of the algorithm.`
` 323:              // *`
` 324:              // *  Arguments`
` 325:              // *  =========`
` 326:              // *`
` 327:              // *  UPLO    (input) CHARACTER*1`
` 328:              // *          = 'U':  B is upper bidiagonal;`
` 329:              // *          = 'L':  B is lower bidiagonal.`
` 330:              // *`
` 331:              // *  N       (input) INTEGER`
` 332:              // *          The order of the matrix B.  N >= 0.`
` 333:              // *`
` 334:              // *  NCVT    (input) INTEGER`
` 335:              // *          The number of columns of the matrix VT. NCVT >= 0.`
` 336:              // *`
` 337:              // *  NRU     (input) INTEGER`
` 338:              // *          The number of rows of the matrix U. NRU >= 0.`
` 339:              // *`
` 340:              // *  NCC     (input) INTEGER`
` 341:              // *          The number of columns of the matrix C. NCC >= 0.`
` 342:              // *`
` 343:              // *  D       (input/output) DOUBLE PRECISION array, dimension (N)`
` 344:              // *          On entry, the n diagonal elements of the bidiagonal matrix B.`
` 345:              // *          On exit, if INFO=0, the singular values of B in decreasing`
` 346:              // *          order.`
` 347:              // *`
` 348:              // *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 349:              // *          On entry, the N-1 offdiagonal elements of the bidiagonal`
` 350:              // *          matrix B. `
` 351:              // *          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E`
` 352:              // *          will contain the diagonal and superdiagonal elements of a`
` 353:              // *          bidiagonal matrix orthogonally equivalent to the one given`
` 354:              // *          as input.`
` 355:              // *`
` 356:              // *  VT      (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)`
` 357:              // *          On entry, an N-by-NCVT matrix VT.`
` 358:              // *          On exit, VT is overwritten by P**T * VT.`
` 359:              // *          Not referenced if NCVT = 0.`
` 360:              // *`
` 361:              // *  LDVT    (input) INTEGER`
` 362:              // *          The leading dimension of the array VT.`
` 363:              // *          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.`
` 364:              // *`
` 365:              // *  U       (input/output) DOUBLE PRECISION array, dimension (LDU, N)`
` 366:              // *          On entry, an NRU-by-N matrix U.`
` 367:              // *          On exit, U is overwritten by U * Q.`
` 368:              // *          Not referenced if NRU = 0.`
` 369:              // *`
` 370:              // *  LDU     (input) INTEGER`
` 371:              // *          The leading dimension of the array U.  LDU >= max(1,NRU).`
` 372:              // *`
` 373:              // *  C       (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)`
` 374:              // *          On entry, an N-by-NCC matrix C.`
` 375:              // *          On exit, C is overwritten by Q**T * C.`
` 376:              // *          Not referenced if NCC = 0.`
` 377:              // *`
` 378:              // *  LDC     (input) INTEGER`
` 379:              // *          The leading dimension of the array C.`
` 380:              // *          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.`
` 381:              // *`
` 382:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)`
` 383:              // *          if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise`
` 384:              // *`
` 385:              // *  INFO    (output) INTEGER`
` 386:              // *          = 0:  successful exit`
` 387:              // *          < 0:  If INFO = -i, the i-th argument had an illegal value`
` 388:              // *          > 0:  the algorithm did not converge; D and E contain the`
` 389:              // *                elements of a bidiagonal matrix which is orthogonally`
` 390:              // *                similar to the input matrix B;  if INFO = i, i`
` 391:              // *                elements of E have not converged to zero.`
` 392:              // *`
` 393:              // *  Internal Parameters`
` 394:              // *  ===================`
` 395:              // *`
` 396:              // *  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))`
` 397:              // *          TOLMUL controls the convergence criterion of the QR loop.`
` 398:              // *          If it is positive, TOLMUL*EPS is the desired relative`
` 399:              // *             precision in the computed singular values.`
` 400:              // *          If it is negative, abs(TOLMUL*EPS*sigma_max) is the`
` 401:              // *             desired absolute accuracy in the computed singular`
` 402:              // *             values (corresponds to relative accuracy`
` 403:              // *             abs(TOLMUL*EPS) in the largest singular value.`
` 404:              // *          abs(TOLMUL) should be between 1 and 1/EPS, and preferably`
` 405:              // *             between 10 (for fast convergence) and .1/EPS`
` 406:              // *             (for there to be some accuracy in the results).`
` 407:              // *          Default is to lose at either one eighth or 2 of the`
` 408:              // *             available decimal digits in each computed singular value`
` 409:              // *             (whichever is smaller).`
` 410:              // *`
` 411:              // *  MAXITR  INTEGER, default = 6`
` 412:              // *          MAXITR controls the maximum number of passes of the`
` 413:              // *          algorithm through its inner loop. The algorithms stops`
` 414:              // *          (and so fails to converge) if the number of passes`
` 415:              // *          through the inner loop exceeds MAXITR*N**2.`
` 416:              // *`
` 417:              // *  =====================================================================`
` 418:              // *`
` 419:              // *     .. Parameters ..`
` 420:              // *     ..`
` 421:              // *     .. Local Scalars ..`
` 422:              // *     ..`
` 423:              // *     .. External Functions ..`
` 424:              // *     ..`
` 425:              // *     .. External Subroutines ..`
` 426:              // *     ..`
` 427:              // *     .. Intrinsic Functions ..`
` 428:              //      INTRINSIC          ABS, DBLE, MAX, MIN, SIGN, SQRT;`
` 429:              // *     ..`
` 430:              // *     .. Executable Statements ..`
` 431:              // *`
` 432:              // *     Test the input parameters.`
` 433:              // *`
` 434:   `
` 435:              #endregion`
` 436:   `
` 437:   `
` 438:              #region Body`
` 439:              `
` 440:              INFO = 0;`
` 441:              LOWER = this._lsame.Run(UPLO, "L");`
` 442:              if (!this._lsame.Run(UPLO, "U") && !LOWER)`
` 443:              {`
` 444:                  INFO =  - 1;`
` 445:              }`
` 446:              else`
` 447:              {`
` 448:                  if (N < 0)`
` 449:                  {`
` 450:                      INFO =  - 2;`
` 451:                  }`
` 452:                  else`
` 453:                  {`
` 454:                      if (NCVT < 0)`
` 455:                      {`
` 456:                          INFO =  - 3;`
` 457:                      }`
` 458:                      else`
` 459:                      {`
` 460:                          if (NRU < 0)`
` 461:                          {`
` 462:                              INFO =  - 4;`
` 463:                          }`
` 464:                          else`
` 465:                          {`
` 466:                              if (NCC < 0)`
` 467:                              {`
` 468:                                  INFO =  - 5;`
` 469:                              }`
` 470:                              else`
` 471:                              {`
` 472:                                  if ((NCVT == 0 && LDVT < 1) || (NCVT > 0 && LDVT < Math.Max(1, N)))`
` 473:                                  {`
` 474:                                      INFO =  - 9;`
` 475:                                  }`
` 476:                                  else`
` 477:                                  {`
` 478:                                      if (LDU < Math.Max(1, NRU))`
` 479:                                      {`
` 480:                                          INFO =  - 11;`
` 481:                                      }`
` 482:                                      else`
` 483:                                      {`
` 484:                                          if ((NCC == 0 && LDC < 1) || (NCC > 0 && LDC < Math.Max(1, N)))`
` 485:                                          {`
` 486:                                              INFO =  - 13;`
` 487:                                          }`
` 488:                                      }`
` 489:                                  }`
` 490:                              }`
` 491:                          }`
` 492:                      }`
` 493:                  }`
` 494:              }`
` 495:              if (INFO != 0)`
` 496:              {`
` 497:                  this._xerbla.Run("DBDSQR",  - INFO);`
` 498:                  return;`
` 499:              }`
` 500:              if (N == 0) return;`
` 501:              if (N == 1) goto LABEL160;`
` 502:              // *`
` 503:              // *     ROTATE is true if any singular vectors desired, false otherwise`
` 504:              // *`
` 505:              ROTATE = (NCVT > 0) || (NRU > 0) || (NCC > 0);`
` 506:              // *`
` 507:              // *     If no singular vectors desired, use qd algorithm`
` 508:              // *`
` 509:              if (!ROTATE)`
` 510:              {`
` 511:                  this._dlasq1.Run(N, ref D, offset_d, E, offset_e, ref WORK, offset_work, ref INFO);`
` 512:                  return;`
` 513:              }`
` 514:              // *`
` 515:              NM1 = N - 1;`
` 516:              NM12 = NM1 + NM1;`
` 517:              NM13 = NM12 + NM1;`
` 518:              IDIR = 0;`
` 519:              // *`
` 520:              // *     Get machine constants`
` 521:              // *`
` 522:              EPS = this._dlamch.Run("Epsilon");`
` 523:              UNFL = this._dlamch.Run("Safe minimum");`
` 524:              // *`
` 525:              // *     If matrix lower bidiagonal, rotate to be upper bidiagonal`
` 526:              // *     by applying Givens rotations on the left`
` 527:              // *`
` 528:              if (LOWER)`
` 529:              {`
` 530:                  for (I = 1; I <= N - 1; I++)`
` 531:                  {`
` 532:                      this._dlartg.Run(D[I + o_d], E[I + o_e], ref CS, ref SN, ref R);`
` 533:                      D[I + o_d] = R;`
` 534:                      E[I + o_e] = SN * D[I + 1 + o_d];`
` 535:                      D[I + 1 + o_d] = CS * D[I + 1 + o_d];`
` 536:                      WORK[I + o_work] = CS;`
` 537:                      WORK[NM1 + I + o_work] = SN;`
` 538:                  }`
` 539:                  // *`
` 540:                  // *        Update singular vectors if desired`
` 541:                  // *`
` 542:                  if (NRU > 0)`
` 543:                  {`
` 544:                      this._dlasr.Run("R", "V", "F", NRU, N, WORK, 1 + o_work`
` 545:                                      , WORK, N + o_work, ref U, offset_u, LDU);`
` 546:                  }`
` 547:                  if (NCC > 0)`
` 548:                  {`
` 549:                      this._dlasr.Run("L", "V", "F", N, NCC, WORK, 1 + o_work`
` 550:                                      , WORK, N + o_work, ref C, offset_c, LDC);`
` 551:                  }`
` 552:              }`
` 553:              // *`
` 554:              // *     Compute singular values to relative accuracy TOL`
` 555:              // *     (By setting TOL to be negative, algorithm will compute`
` 556:              // *     singular values to absolute accuracy ABS(TOL)*norm(input matrix))`
` 557:              // *`
` 558:              TOLMUL = Math.Max(TEN, Math.Min(HNDRD, Math.Pow(EPS,MEIGTH)));`
` 559:              TOL = TOLMUL * EPS;`
` 560:              // *`
` 561:              // *     Compute approximate maximum, minimum singular values`
` 562:              // *`
` 563:              SMAX = ZERO;`
` 564:              for (I = 1; I <= N; I++)`
` 565:              {`
` 566:                  SMAX = Math.Max(SMAX, Math.Abs(D[I + o_d]));`
` 567:              }`
` 568:              for (I = 1; I <= N - 1; I++)`
` 569:              {`
` 570:                  SMAX = Math.Max(SMAX, Math.Abs(E[I + o_e]));`
` 571:              }`
` 572:              SMINL = ZERO;`
` 573:              if (TOL >= ZERO)`
` 574:              {`
` 575:                  // *`
` 576:                  // *        Relative accuracy desired`
` 577:                  // *`
` 578:                  SMINOA = Math.Abs(D[1 + o_d]);`
` 579:                  if (SMINOA == ZERO) goto LABEL50;`
` 580:                  MU = SMINOA;`
` 581:                  for (I = 2; I <= N; I++)`
` 582:                  {`
` 583:                      MU = Math.Abs(D[I + o_d]) * (MU / (MU + Math.Abs(E[I - 1 + o_e])));`
` 584:                      SMINOA = Math.Min(SMINOA, MU);`
` 585:                      if (SMINOA == ZERO) goto LABEL50;`
` 586:                  }`
` 587:              LABEL50:;`
` 588:                  SMINOA = SMINOA / Math.Sqrt(Convert.ToDouble(N));`
` 589:                  THRESH = Math.Max(TOL * SMINOA, MAXITR * N * N * UNFL);`
` 590:              }`
` 591:              else`
` 592:              {`
` 593:                  // *`
` 594:                  // *        Absolute accuracy desired`
` 595:                  // *`
` 596:                  THRESH = Math.Max(Math.Abs(TOL) * SMAX, MAXITR * N * N * UNFL);`
` 597:              }`
` 598:              // *`
` 599:              // *     Prepare for main iteration loop for the singular values`
` 600:              // *     (MAXIT is the maximum number of passes through the inner`
` 601:              // *     loop permitted before nonconvergence signalled.)`
` 602:              // *`
` 603:              MAXIT = MAXITR * N * N;`
` 604:              ITER = 0;`
` 605:              OLDLL =  - 1;`
` 606:              OLDM =  - 1;`
` 607:              // *`
` 608:              // *     M points to last element of unconverged part of matrix`
` 609:              // *`
` 610:              M = N;`
` 611:              // *`
` 612:              // *     Begin main iteration loop`
` 613:              // *`
` 614:          LABEL60:;`
` 615:              // *`
` 616:              // *     Check for convergence or exceeding iteration count`
` 617:              // *`
` 618:              if (M <= 1) goto LABEL160;`
` 619:              if (ITER > MAXIT) goto LABEL200;`
` 620:              // *`
` 621:              // *     Find diagonal block of matrix to work on`
` 622:              // *`
` 623:              if (TOL < ZERO && Math.Abs(D[M + o_d]) <= THRESH) D[M + o_d] = ZERO;`
` 624:              SMAX = Math.Abs(D[M + o_d]);`
` 625:              SMIN = SMAX;`
` 626:              for (LLL = 1; LLL <= M - 1; LLL++)`
` 627:              {`
` 628:                  LL = M - LLL;`
` 629:                  ABSS = Math.Abs(D[LL + o_d]);`
` 630:                  ABSE = Math.Abs(E[LL + o_e]);`
` 631:                  if (TOL < ZERO && ABSS <= THRESH) D[LL + o_d] = ZERO;`
` 632:                  if (ABSE <= THRESH) goto LABEL80;`
` 633:                  SMIN = Math.Min(SMIN, ABSS);`
` 634:                  SMAX = Math.Max(SMAX, Math.Max(ABSS, ABSE));`
` 635:              }`
` 636:              LL = 0;`
` 637:              goto LABEL90;`
` 638:          LABEL80:;`
` 639:              E[LL + o_e] = ZERO;`
` 640:              // *`
` 641:              // *     Matrix splits since E(LL) = 0`
` 642:              // *`
` 643:              if (LL == M - 1)`
` 644:              {`
` 645:                  // *`
` 646:                  // *        Convergence of bottom singular value, return to top of loop`
` 647:                  // *`
` 648:                  M = M - 1;`
` 649:                  goto LABEL60;`
` 650:              }`
` 651:          LABEL90:;`
` 652:              LL = LL + 1;`
` 653:              // *`
` 654:              // *     E(LL) through E(M-1) are nonzero, E(LL-1) is zero`
` 655:              // *`
` 656:              if (LL == M - 1)`
` 657:              {`
` 658:                  // *`
` 659:                  // *        2 by 2 block, handle separately`
` 660:                  // *`
` 661:                  this._dlasv2.Run(D[M - 1 + o_d], E[M - 1 + o_e], D[M + o_d], ref SIGMN, ref SIGMX, ref SINR`
` 662:                                   , ref COSR, ref SINL, ref COSL);`
` 663:                  D[M - 1 + o_d] = SIGMX;`
` 664:                  E[M - 1 + o_e] = ZERO;`
` 665:                  D[M + o_d] = SIGMN;`
` 666:                  // *`
` 667:                  // *        Compute singular vectors, if desired`
` 668:                  // *`
` 669:                  if (NCVT > 0)`
` 670:                  {`
` 671:                      this._drot.Run(NCVT, ref VT, M - 1+1 * LDVT + o_vt, LDVT, ref VT, M+1 * LDVT + o_vt, LDVT, COSR`
` 672:                                     , SINR);`
` 673:                  }`
` 674:                  if (NRU > 0)`
` 675:                  {`
` 676:                      this._drot.Run(NRU, ref U, 1+(M - 1) * LDU + o_u, 1, ref U, 1+M * LDU + o_u, 1, COSL`
` 677:                                     , SINL);`
` 678:                  }`
` 679:                  if (NCC > 0)`
` 680:                  {`
` 681:                      this._drot.Run(NCC, ref C, M - 1+1 * LDC + o_c, LDC, ref C, M+1 * LDC + o_c, LDC, COSL`
` 682:                                     , SINL);`
` 683:                  }`
` 684:                  M = M - 2;`
` 685:                  goto LABEL60;`
` 686:              }`
` 687:              // *`
` 688:              // *     If working on new submatrix, choose shift direction`
` 689:              // *     (from larger end diagonal element towards smaller)`
` 690:              // *`
` 691:              if (LL > OLDM || M < OLDLL)`
` 692:              {`
` 693:                  if (Math.Abs(D[LL + o_d]) >= Math.Abs(D[M + o_d]))`
` 694:                  {`
` 695:                      // *`
` 696:                      // *           Chase bulge from top (big end) to bottom (small end)`
` 697:                      // *`
` 698:                      IDIR = 1;`
` 699:                  }`
` 700:                  else`
` 701:                  {`
` 702:                      // *`
` 703:                      // *           Chase bulge from bottom (big end) to top (small end)`
` 704:                      // *`
` 705:                      IDIR = 2;`
` 706:                  }`
` 707:              }`
` 708:              // *`
` 709:              // *     Apply convergence tests`
` 710:              // *`
` 711:              if (IDIR == 1)`
` 712:              {`
` 713:                  // *`
` 714:                  // *        Run convergence test in forward direction`
` 715:                  // *        First apply standard test to bottom of matrix`
` 716:                  // *`
` 717:                  if (Math.Abs(E[M - 1 + o_e]) <= Math.Abs(TOL) * Math.Abs(D[M + o_d]) || (TOL < ZERO && Math.Abs(E[M - 1 + o_e]) <= THRESH))`
` 718:                  {`
` 719:                      E[M - 1 + o_e] = ZERO;`
` 720:                      goto LABEL60;`
` 721:                  }`
` 722:                  // *`
` 723:                  if (TOL >= ZERO)`
` 724:                  {`
` 725:                      // *`
` 726:                      // *           If relative accuracy desired,`
` 727:                      // *           apply convergence criterion forward`
` 728:                      // *`
` 729:                      MU = Math.Abs(D[LL + o_d]);`
` 730:                      SMINL = MU;`
` 731:                      for (LLL = LL; LLL <= M - 1; LLL++)`
` 732:                      {`
` 733:                          if (Math.Abs(E[LLL + o_e]) <= TOL * MU)`
` 734:                          {`
` 735:                              E[LLL + o_e] = ZERO;`
` 736:                              goto LABEL60;`
` 737:                          }`
` 738:                          MU = Math.Abs(D[LLL + 1 + o_d]) * (MU / (MU + Math.Abs(E[LLL + o_e])));`
` 739:                          SMINL = Math.Min(SMINL, MU);`
` 740:                      }`
` 741:                  }`
` 742:                  // *`
` 743:              }`
` 744:              else`
` 745:              {`
` 746:                  // *`
` 747:                  // *        Run convergence test in backward direction`
` 748:                  // *        First apply standard test to top of matrix`
` 749:                  // *`
` 750:                  if (Math.Abs(E[LL + o_e]) <= Math.Abs(TOL) * Math.Abs(D[LL + o_d]) || (TOL < ZERO && Math.Abs(E[LL + o_e]) <= THRESH))`
` 751:                  {`
` 752:                      E[LL + o_e] = ZERO;`
` 753:                      goto LABEL60;`
` 754:                  }`
` 755:                  // *`
` 756:                  if (TOL >= ZERO)`
` 757:                  {`
` 758:                      // *`
` 759:                      // *           If relative accuracy desired,`
` 760:                      // *           apply convergence criterion backward`
` 761:                      // *`
` 762:                      MU = Math.Abs(D[M + o_d]);`
` 763:                      SMINL = MU;`
` 764:                      for (LLL = M - 1; LLL >= LL; LLL +=  - 1)`
` 765:                      {`
` 766:                          if (Math.Abs(E[LLL + o_e]) <= TOL * MU)`
` 767:                          {`
` 768:                              E[LLL + o_e] = ZERO;`
` 769:                              goto LABEL60;`
` 770:                          }`
` 771:                          MU = Math.Abs(D[LLL + o_d]) * (MU / (MU + Math.Abs(E[LLL + o_e])));`
` 772:                          SMINL = Math.Min(SMINL, MU);`
` 773:                      }`
` 774:                  }`
` 775:              }`
` 776:              OLDLL = LL;`
` 777:              OLDM = M;`
` 778:              // *`
` 779:              // *     Compute shift.  First, test if shifting would ruin relative`
` 780:              // *     accuracy, and if so set the shift to zero.`
` 781:              // *`
` 782:              if (TOL >= ZERO && N * TOL * (SMINL / SMAX) <= Math.Max(EPS, HNDRTH * TOL))`
` 783:              {`
` 784:                  // *`
` 785:                  // *        Use a zero shift to avoid loss of relative accuracy`
` 786:                  // *`
` 787:                  SHIFT = ZERO;`
` 788:              }`
` 789:              else`
` 790:              {`
` 791:                  // *`
` 792:                  // *        Compute the shift from 2-by-2 block at end of matrix`
` 793:                  // *`
` 794:                  if (IDIR == 1)`
` 795:                  {`
` 796:                      SLL = Math.Abs(D[LL + o_d]);`
` 797:                      this._dlas2.Run(D[M - 1 + o_d], E[M - 1 + o_e], D[M + o_d], ref SHIFT, ref R);`
` 798:                  }`
` 799:                  else`
` 800:                  {`
` 801:                      SLL = Math.Abs(D[M + o_d]);`
` 802:                      this._dlas2.Run(D[LL + o_d], E[LL + o_e], D[LL + 1 + o_d], ref SHIFT, ref R);`
` 803:                  }`
` 804:                  // *`
` 805:                  // *        Test if shift negligible, and if so set to zero`
` 806:                  // *`
` 807:                  if (SLL > ZERO)`
` 808:                  {`
` 809:                      if (Math.Pow(SHIFT / SLL,2) < EPS) SHIFT = ZERO;`
` 810:                  }`
` 811:              }`
` 812:              // *`
` 813:              // *     Increment iteration count`
` 814:              // *`
` 815:              ITER = ITER + M - LL;`
` 816:              // *`
` 817:              // *     If SHIFT = 0, do simplified QR iteration`
` 818:              // *`
` 819:              if (SHIFT == ZERO)`
` 820:              {`
` 821:                  if (IDIR == 1)`
` 822:                  {`
` 823:                      // *`
` 824:                      // *           Chase bulge from top to bottom`
` 825:                      // *           Save cosines and sines for later singular vector updates`
` 826:                      // *`
` 827:                      CS = ONE;`
` 828:                      OLDCS = ONE;`
` 829:                      for (I = LL; I <= M - 1; I++)`
` 830:                      {`
` 831:                          this._dlartg.Run(D[I + o_d] * CS, E[I + o_e], ref CS, ref SN, ref R);`
` 832:                          if (I > LL) E[I - 1 + o_e] = OLDSN * R;`
` 833:                          this._dlartg.Run(OLDCS * R, D[I + 1 + o_d] * SN, ref OLDCS, ref OLDSN, ref D[I + o_d]);`
` 834:                          WORK[I - LL + 1 + o_work] = CS;`
` 835:                          WORK[I - LL + 1 + NM1 + o_work] = SN;`
` 836:                          WORK[I - LL + 1 + NM12 + o_work] = OLDCS;`
` 837:                          WORK[I - LL + 1 + NM13 + o_work] = OLDSN;`
` 838:                      }`
` 839:                      H = D[M + o_d] * CS;`
` 840:                      D[M + o_d] = H * OLDCS;`
` 841:                      E[M - 1 + o_e] = H * OLDSN;`
` 842:                      // *`
` 843:                      // *           Update singular vectors`
` 844:                      // *`
` 845:                      if (NCVT > 0)`
` 846:                      {`
` 847:                          this._dlasr.Run("L", "V", "F", M - LL + 1, NCVT, WORK, 1 + o_work`
` 848:                                          , WORK, N + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);`
` 849:                      }`
` 850:                      if (NRU > 0)`
` 851:                      {`
` 852:                          this._dlasr.Run("R", "V", "F", NRU, M - LL + 1, WORK, NM12 + 1 + o_work`
` 853:                                          , WORK, NM13 + 1 + o_work, ref U, 1+LL * LDU + o_u, LDU);`
` 854:                      }`
` 855:                      if (NCC > 0)`
` 856:                      {`
` 857:                          this._dlasr.Run("L", "V", "F", M - LL + 1, NCC, WORK, NM12 + 1 + o_work`
` 858:                                          , WORK, NM13 + 1 + o_work, ref C, LL+1 * LDC + o_c, LDC);`
` 859:                      }`
` 860:                      // *`
` 861:                      // *           Test convergence`
` 862:                      // *`
` 863:                      if (Math.Abs(E[M - 1 + o_e]) <= THRESH) E[M - 1 + o_e] = ZERO;`
` 864:                      // *`
` 865:                  }`
` 866:                  else`
` 867:                  {`
` 868:                      // *`
` 869:                      // *           Chase bulge from bottom to top`
` 870:                      // *           Save cosines and sines for later singular vector updates`
` 871:                      // *`
` 872:                      CS = ONE;`
` 873:                      OLDCS = ONE;`
` 874:                      for (I = M; I >= LL + 1; I +=  - 1)`
` 875:                      {`
` 876:                          this._dlartg.Run(D[I + o_d] * CS, E[I - 1 + o_e], ref CS, ref SN, ref R);`
` 877:                          if (I < M) E[I + o_e] = OLDSN * R;`
` 878:                          this._dlartg.Run(OLDCS * R, D[I - 1 + o_d] * SN, ref OLDCS, ref OLDSN, ref D[I + o_d]);`
` 879:                          WORK[I - LL + o_work] = CS;`
` 880:                          WORK[I - LL + NM1 + o_work] =  - SN;`
` 881:                          WORK[I - LL + NM12 + o_work] = OLDCS;`
` 882:                          WORK[I - LL + NM13 + o_work] =  - OLDSN;`
` 883:                      }`
` 884:                      H = D[LL + o_d] * CS;`
` 885:                      D[LL + o_d] = H * OLDCS;`
` 886:                      E[LL + o_e] = H * OLDSN;`
` 887:                      // *`
` 888:                      // *           Update singular vectors`
` 889:                      // *`
` 890:                      if (NCVT > 0)`
` 891:                      {`
` 892:                          this._dlasr.Run("L", "V", "B", M - LL + 1, NCVT, WORK, NM12 + 1 + o_work`
` 893:                                          , WORK, NM13 + 1 + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);`
` 894:                      }`
` 895:                      if (NRU > 0)`
` 896:                      {`
` 897:                          this._dlasr.Run("R", "V", "B", NRU, M - LL + 1, WORK, 1 + o_work`
` 898:                                          , WORK, N + o_work, ref U, 1+LL * LDU + o_u, LDU);`
` 899:                      }`
` 900:                      if (NCC > 0)`
` 901:                      {`
` 902:                          this._dlasr.Run("L", "V", "B", M - LL + 1, NCC, WORK, 1 + o_work`
` 903:                                          , WORK, N + o_work, ref C, LL+1 * LDC + o_c, LDC);`
` 904:                      }`
` 905:                      // *`
` 906:                      // *           Test convergence`
` 907:                      // *`
` 908:                      if (Math.Abs(E[LL + o_e]) <= THRESH) E[LL + o_e] = ZERO;`
` 909:                  }`
` 910:              }`
` 911:              else`
` 912:              {`
` 913:                  // *`
` 914:                  // *        Use nonzero shift`
` 915:                  // *`
` 916:                  if (IDIR == 1)`
` 917:                  {`
` 918:                      // *`
` 919:                      // *           Chase bulge from top to bottom`
` 920:                      // *           Save cosines and sines for later singular vector updates`
` 921:                      // *`
` 922:                      F = (Math.Abs(D[LL + o_d]) - SHIFT) * (FortranLib.Sign(ONE,D[LL + o_d]) + SHIFT / D[LL + o_d]);`
` 923:                      G = E[LL + o_e];`
` 924:                      for (I = LL; I <= M - 1; I++)`
` 925:                      {`
` 926:                          this._dlartg.Run(F, G, ref COSR, ref SINR, ref R);`
` 927:                          if (I > LL) E[I - 1 + o_e] = R;`
` 928:                          F = COSR * D[I + o_d] + SINR * E[I + o_e];`
` 929:                          E[I + o_e] = COSR * E[I + o_e] - SINR * D[I + o_d];`
` 930:                          G = SINR * D[I + 1 + o_d];`
` 931:                          D[I + 1 + o_d] = COSR * D[I + 1 + o_d];`
` 932:                          this._dlartg.Run(F, G, ref COSL, ref SINL, ref R);`
` 933:                          D[I + o_d] = R;`
` 934:                          F = COSL * E[I + o_e] + SINL * D[I + 1 + o_d];`
` 935:                          D[I + 1 + o_d] = COSL * D[I + 1 + o_d] - SINL * E[I + o_e];`
` 936:                          if (I < M - 1)`
` 937:                          {`
` 938:                              G = SINL * E[I + 1 + o_e];`
` 939:                              E[I + 1 + o_e] = COSL * E[I + 1 + o_e];`
` 940:                          }`
` 941:                          WORK[I - LL + 1 + o_work] = COSR;`
` 942:                          WORK[I - LL + 1 + NM1 + o_work] = SINR;`
` 943:                          WORK[I - LL + 1 + NM12 + o_work] = COSL;`
` 944:                          WORK[I - LL + 1 + NM13 + o_work] = SINL;`
` 945:                      }`
` 946:                      E[M - 1 + o_e] = F;`
` 947:                      // *`
` 948:                      // *           Update singular vectors`
` 949:                      // *`
` 950:                      if (NCVT > 0)`
` 951:                      {`
` 952:                          this._dlasr.Run("L", "V", "F", M - LL + 1, NCVT, WORK, 1 + o_work`
` 953:                                          , WORK, N + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);`
` 954:                      }`
` 955:                      if (NRU > 0)`
` 956:                      {`
` 957:                          this._dlasr.Run("R", "V", "F", NRU, M - LL + 1, WORK, NM12 + 1 + o_work`
` 958:                                          , WORK, NM13 + 1 + o_work, ref U, 1+LL * LDU + o_u, LDU);`
` 959:                      }`
` 960:                      if (NCC > 0)`
` 961:                      {`
` 962:                          this._dlasr.Run("L", "V", "F", M - LL + 1, NCC, WORK, NM12 + 1 + o_work`
` 963:                                          , WORK, NM13 + 1 + o_work, ref C, LL+1 * LDC + o_c, LDC);`
` 964:                      }`
` 965:                      // *`
` 966:                      // *           Test convergence`
` 967:                      // *`
` 968:                      if (Math.Abs(E[M - 1 + o_e]) <= THRESH) E[M - 1 + o_e] = ZERO;`
` 969:                      // *`
` 970:                  }`
` 971:                  else`
` 972:                  {`
` 973:                      // *`
` 974:                      // *           Chase bulge from bottom to top`
` 975:                      // *           Save cosines and sines for later singular vector updates`
` 976:                      // *`
` 977:                      F = (Math.Abs(D[M + o_d]) - SHIFT) * (FortranLib.Sign(ONE,D[M + o_d]) + SHIFT / D[M + o_d]);`
` 978:                      G = E[M - 1 + o_e];`
` 979:                      for (I = M; I >= LL + 1; I +=  - 1)`
` 980:                      {`
` 981:                          this._dlartg.Run(F, G, ref COSR, ref SINR, ref R);`
` 982:                          if (I < M) E[I + o_e] = R;`
` 983:                          F = COSR * D[I + o_d] + SINR * E[I - 1 + o_e];`
` 984:                          E[I - 1 + o_e] = COSR * E[I - 1 + o_e] - SINR * D[I + o_d];`
` 985:                          G = SINR * D[I - 1 + o_d];`
` 986:                          D[I - 1 + o_d] = COSR * D[I - 1 + o_d];`
` 987:                          this._dlartg.Run(F, G, ref COSL, ref SINL, ref R);`
` 988:                          D[I + o_d] = R;`
` 989:                          F = COSL * E[I - 1 + o_e] + SINL * D[I - 1 + o_d];`
` 990:                          D[I - 1 + o_d] = COSL * D[I - 1 + o_d] - SINL * E[I - 1 + o_e];`
` 991:                          if (I > LL + 1)`
` 992:                          {`
` 993:                              G = SINL * E[I - 2 + o_e];`
` 994:                              E[I - 2 + o_e] = COSL * E[I - 2 + o_e];`
` 995:                          }`
` 996:                          WORK[I - LL + o_work] = COSR;`
` 997:                          WORK[I - LL + NM1 + o_work] =  - SINR;`
` 998:                          WORK[I - LL + NM12 + o_work] = COSL;`
` 999:                          WORK[I - LL + NM13 + o_work] =  - SINL;`
`1000:                      }`
`1001:                      E[LL + o_e] = F;`
`1002:                      // *`
`1003:                      // *           Test convergence`
`1004:                      // *`
`1005:                      if (Math.Abs(E[LL + o_e]) <= THRESH) E[LL + o_e] = ZERO;`
`1006:                      // *`
`1007:                      // *           Update singular vectors if desired`
`1008:                      // *`
`1009:                      if (NCVT > 0)`
`1010:                      {`
`1011:                          this._dlasr.Run("L", "V", "B", M - LL + 1, NCVT, WORK, NM12 + 1 + o_work`
`1012:                                          , WORK, NM13 + 1 + o_work, ref VT, LL+1 * LDVT + o_vt, LDVT);`
`1013:                      }`
`1014:                      if (NRU > 0)`
`1015:                      {`
`1016:                          this._dlasr.Run("R", "V", "B", NRU, M - LL + 1, WORK, 1 + o_work`
`1017:                                          , WORK, N + o_work, ref U, 1+LL * LDU + o_u, LDU);`
`1018:                      }`
`1019:                      if (NCC > 0)`
`1020:                      {`
`1021:                          this._dlasr.Run("L", "V", "B", M - LL + 1, NCC, WORK, 1 + o_work`
`1022:                                          , WORK, N + o_work, ref C, LL+1 * LDC + o_c, LDC);`
`1023:                      }`
`1024:                  }`
`1025:              }`
`1026:              // *`
`1027:              // *     QR iteration finished, go back and check convergence`
`1028:              // *`
`1029:              goto LABEL60;`
`1030:              // *`
`1031:              // *     All singular values converged, so make them positive`
`1032:              // *`
`1033:          LABEL160:;`
`1034:              for (I = 1; I <= N; I++)`
`1035:              {`
`1036:                  if (D[I + o_d] < ZERO)`
`1037:                  {`
`1038:                      D[I + o_d] =  - D[I + o_d];`
`1039:                      // *`
`1040:                      // *           Change sign of singular vectors, if desired`
`1041:                      // *`
`1042:                      if (NCVT > 0) this._dscal.Run(NCVT, NEGONE, ref VT, I+1 * LDVT + o_vt, LDVT);`
`1043:                  }`
`1044:              }`
`1045:              // *`
`1046:              // *     Sort the singular values into decreasing order (insertion sort on`
`1047:              // *     singular values, but only one transposition per singular vector)`
`1048:              // *`
`1049:              for (I = 1; I <= N - 1; I++)`
`1050:              {`
`1051:                  // *`
`1052:                  // *        Scan for smallest D(I)`
`1053:                  // *`
`1054:                  ISUB = 1;`
`1055:                  SMIN = D[1 + o_d];`
`1056:                  for (J = 2; J <= N + 1 - I; J++)`
`1057:                  {`
`1058:                      if (D[J + o_d] <= SMIN)`
`1059:                      {`
`1060:                          ISUB = J;`
`1061:                          SMIN = D[J + o_d];`
`1062:                      }`
`1063:                  }`
`1064:                  if (ISUB != N + 1 - I)`
`1065:                  {`
`1066:                      // *`
`1067:                      // *           Swap singular values and vectors`
`1068:                      // *`
`1069:                      D[ISUB + o_d] = D[N + 1 - I + o_d];`
`1070:                      D[N + 1 - I + o_d] = SMIN;`
`1071:                      if (NCVT > 0) this._dswap.Run(NCVT, ref VT, ISUB+1 * LDVT + o_vt, LDVT, ref VT, N + 1 - I+1 * LDVT + o_vt, LDVT);`
`1072:                      if (NRU > 0) this._dswap.Run(NRU, ref U, 1+ISUB * LDU + o_u, 1, ref U, 1+(N + 1 - I) * LDU + o_u, 1);`
`1073:                      if (NCC > 0) this._dswap.Run(NCC, ref C, ISUB+1 * LDC + o_c, LDC, ref C, N + 1 - I+1 * LDC + o_c, LDC);`
`1074:                  }`
`1075:              }`
`1076:              goto LABEL220;`
`1077:              // *`
`1078:              // *     Maximum number of iterations exceeded, failure to converge`
`1079:              // *`
`1080:          LABEL200:;`
`1081:              INFO = 0;`
`1082:              for (I = 1; I <= N - 1; I++)`
`1083:              {`
`1084:                  if (E[I + o_e] != ZERO) INFO = INFO + 1;`
`1085:              }`
`1086:          LABEL220:;`
`1087:              return;`
`1088:              // *`
`1089:              // *     End of DBDSQR`
`1090:              // *`
`1091:   `
`1092:              #endregion`
`1093:   `
`1094:          }`
`1095:      }`
`1096:  }`