Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   1:  #region Translated by Jose Antonio De Santiago-Castillo.
   2:   
   3:  //Translated by Jose Antonio De Santiago-Castillo. 
   4:  //E-mail:JAntonioDeSantiago@gmail.com
   5:  //Web: www.DotNumerics.com
   6:  //
   7:  //Fortran to C# Translation.
   8:  //Translated by:
   9:  //F2CSharp Version 0.71 (November 10, 2009)
  10:  //Code Optimizations: None
  11:  //
  12:  #endregion
  13:   
  14:  using System;
  15:  using DotNumerics.FortranLibrary;
  16:   
  17:  namespace DotNumerics.CSLapack
  18:  {
  19:      /// <summary>
  20:      /// -- LAPACK 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:  }