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 auxiliary routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLASD3 finds all the square roots of the roots of the secular
  27:      /// equation, as defined by the values in D and Z.  It makes the
  28:      /// appropriate calls to DLASD4 and then updates the singular
  29:      /// vectors by matrix multiplication.
  30:      /// 
  31:      /// This code makes very mild assumptions about floating point
  32:      /// arithmetic. It will work on machines with a guard digit in
  33:      /// add/subtract, or on those binary machines without guard digits
  34:      /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
  35:      /// It could conceivably fail on hexadecimal or decimal machines
  36:      /// without guard digits, but we know of none.
  37:      /// 
  38:      /// DLASD3 is called from DLASD1.
  39:      /// 
  40:      ///</summary>
  41:      public class DLASD3
  42:      {
  43:      
  44:   
  45:          #region Dependencies
  46:          
  47:          DLAMC3 _dlamc3; DNRM2 _dnrm2; DCOPY _dcopy; DGEMM _dgemm; DLACPY _dlacpy; DLASCL _dlascl; DLASD4 _dlasd4; XERBLA _xerbla; 
  48:   
  49:          #endregion
  50:   
  51:   
  52:          #region Fields
  53:          
  54:          const double ONE = 1.0E+0; const double ZERO = 0.0E+0; const double NEGONE =  - 1.0E+0; int CTEMP = 0; int I = 0; 
  55:          int J = 0;int JC = 0; int KTEMP = 0; int M = 0; int N = 0; int NLP1 = 0; int NLP2 = 0; int NRP1 = 0; double RHO = 0; 
  56:          double TEMP = 0;
  57:   
  58:          #endregion
  59:   
  60:          public DLASD3(DLAMC3 dlamc3, DNRM2 dnrm2, DCOPY dcopy, DGEMM dgemm, DLACPY dlacpy, DLASCL dlascl, DLASD4 dlasd4, XERBLA xerbla)
  61:          {
  62:      
  63:   
  64:              #region Set Dependencies
  65:              
  66:              this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy; 
  67:              this._dlascl = dlascl;this._dlasd4 = dlasd4; this._xerbla = xerbla; 
  68:   
  69:              #endregion
  70:   
  71:          }
  72:      
  73:          public DLASD3()
  74:          {
  75:      
  76:   
  77:              #region Dependencies (Initialization)
  78:              
  79:              DLAMC3 dlamc3 = new DLAMC3();
  80:              DNRM2 dnrm2 = new DNRM2();
  81:              DCOPY dcopy = new DCOPY();
  82:              LSAME lsame = new LSAME();
  83:              XERBLA xerbla = new XERBLA();
  84:              DLASD5 dlasd5 = new DLASD5();
  85:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  86:              DLACPY dlacpy = new DLACPY(lsame);
  87:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  88:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  89:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  90:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  91:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  92:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
  93:              DLAED6 dlaed6 = new DLAED6(dlamch);
  94:              DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
  95:   
  96:              #endregion
  97:   
  98:   
  99:              #region Set Dependencies
 100:              
 101:              this._dlamc3 = dlamc3; this._dnrm2 = dnrm2; this._dcopy = dcopy; this._dgemm = dgemm; this._dlacpy = dlacpy; 
 102:              this._dlascl = dlascl;this._dlasd4 = dlasd4; this._xerbla = xerbla; 
 103:   
 104:              #endregion
 105:   
 106:          }
 107:          /// <summary>
 108:          /// Purpose
 109:          /// =======
 110:          /// 
 111:          /// DLASD3 finds all the square roots of the roots of the secular
 112:          /// equation, as defined by the values in D and Z.  It makes the
 113:          /// appropriate calls to DLASD4 and then updates the singular
 114:          /// vectors by matrix multiplication.
 115:          /// 
 116:          /// This code makes very mild assumptions about floating point
 117:          /// arithmetic. It will work on machines with a guard digit in
 118:          /// add/subtract, or on those binary machines without guard digits
 119:          /// which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 120:          /// It could conceivably fail on hexadecimal or decimal machines
 121:          /// without guard digits, but we know of none.
 122:          /// 
 123:          /// DLASD3 is called from DLASD1.
 124:          /// 
 125:          ///</summary>
 126:          /// <param name="NL">
 127:          /// (input) INTEGER
 128:          /// The row dimension of the upper block.  NL .GE. 1.
 129:          ///</param>
 130:          /// <param name="NR">
 131:          /// (input) INTEGER
 132:          /// The row dimension of the lower block.  NR .GE. 1.
 133:          ///</param>
 134:          /// <param name="SQRE">
 135:          /// (input) INTEGER
 136:          /// = 0: the lower block is an NR-by-NR square matrix.
 137:          /// = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
 138:          /// 
 139:          /// The bidiagonal matrix has N = NL + NR + 1 rows and
 140:          /// M = N + SQRE .GE. N columns.
 141:          ///</param>
 142:          /// <param name="K">
 143:          /// (input) INTEGER
 144:          /// The size of the secular equation, 1 =.LT. K = .LT. N.
 145:          ///</param>
 146:          /// <param name="D">
 147:          /// (output) DOUBLE PRECISION array, dimension(K)
 148:          /// On exit the square roots of the roots of the secular equation,
 149:          /// in ascending order.
 150:          ///</param>
 151:          /// <param name="Q">
 152:          /// (workspace) DOUBLE PRECISION array,
 153:          /// dimension at least (LDQ,K).
 154:          ///</param>
 155:          /// <param name="LDQ">
 156:          /// (input) INTEGER
 157:          /// The leading dimension of the array Q.  LDQ .GE. K.
 158:          ///</param>
 159:          /// <param name="DSIGMA">
 160:          /// (input) DOUBLE PRECISION array, dimension(K)
 161:          /// The first K elements of this array contain the old roots
 162:          /// of the deflated updating problem.  These are the poles
 163:          /// of the secular equation.
 164:          ///</param>
 165:          /// <param name="U">
 166:          /// (output) DOUBLE PRECISION array, dimension (LDU, N)
 167:          /// The last N - K columns of this matrix contain the deflated
 168:          /// left singular vectors.
 169:          ///</param>
 170:          /// <param name="LDU">
 171:          /// (input) INTEGER
 172:          /// The leading dimension of the array U.  LDU .GE. N.
 173:          ///</param>
 174:          /// <param name="U2">
 175:          /// (input/output) DOUBLE PRECISION array, dimension (LDU2, N)
 176:          /// The first K columns of this matrix contain the non-deflated
 177:          /// left singular vectors for the split problem.
 178:          ///</param>
 179:          /// <param name="LDU2">
 180:          /// (input) INTEGER
 181:          /// The leading dimension of the array U2.  LDU2 .GE. N.
 182:          ///</param>
 183:          /// <param name="VT">
 184:          /// (output) DOUBLE PRECISION array, dimension (LDVT, M)
 185:          /// The last M - K columns of VT' contain the deflated
 186:          /// right singular vectors.
 187:          ///</param>
 188:          /// <param name="LDVT">
 189:          /// (input) INTEGER
 190:          /// The leading dimension of the array VT.  LDVT .GE. N.
 191:          ///</param>
 192:          /// <param name="VT2">
 193:          /// (input/output) DOUBLE PRECISION array, dimension (LDVT2, N)
 194:          /// The first K columns of VT2' contain the non-deflated
 195:          /// right singular vectors for the split problem.
 196:          ///</param>
 197:          /// <param name="LDVT2">
 198:          /// (input) INTEGER
 199:          /// The leading dimension of the array VT2.  LDVT2 .GE. N.
 200:          ///</param>
 201:          /// <param name="IDXC">
 202:          /// (input) INTEGER array, dimension ( N )
 203:          /// The permutation used to arrange the columns of U (and rows of
 204:          /// VT) into three groups:  the first group contains non-zero
 205:          /// entries only at and above (or before) NL +1; the second
 206:          /// contains non-zero entries only at and below (or after) NL+2;
 207:          /// and the third is dense. The first column of U and the row of
 208:          /// VT are treated separately, however.
 209:          /// 
 210:          /// The rows of the singular vectors found by DLASD4
 211:          /// must be likewise permuted before the matrix multiplies can
 212:          /// take place.
 213:          ///</param>
 214:          /// <param name="CTOT">
 215:          /// (input) INTEGER array, dimension ( 4 )
 216:          /// A count of the total number of the various types of columns
 217:          /// in U (or rows in VT), as described in IDXC. The fourth column
 218:          /// type is any column which has been deflated.
 219:          ///</param>
 220:          /// <param name="Z">
 221:          /// (input) DOUBLE PRECISION array, dimension (K)
 222:          /// The first K elements of this array contain the components
 223:          /// of the deflation-adjusted updating row vector.
 224:          ///</param>
 225:          /// <param name="INFO">
 226:          /// (output) INTEGER
 227:          /// = 0:  successful exit.
 228:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 229:          /// .GT. 0:  if INFO = 1, an singular value did not converge
 230:          ///</param>
 231:          public void Run(int NL, int NR, int SQRE, int K, ref double[] D, int offset_d, ref double[] Q, int offset_q
 232:                           , int LDQ, ref double[] DSIGMA, int offset_dsigma, ref double[] U, int offset_u, int LDU, double[] U2, int offset_u2, int LDU2
 233:                           , ref double[] VT, int offset_vt, int LDVT, ref double[] VT2, int offset_vt2, int LDVT2, int[] IDXC, int offset_idxc, int[] CTOT, int offset_ctot
 234:                           , ref double[] Z, int offset_z, ref int INFO)
 235:          {
 236:   
 237:              #region Array Index Correction
 238:              
 239:               int o_d = -1 + offset_d;  int o_q = -1 - LDQ + offset_q;  int o_dsigma = -1 + offset_dsigma; 
 240:               int o_u = -1 - LDU + offset_u; int o_u2 = -1 - LDU2 + offset_u2;  int o_vt = -1 - LDVT + offset_vt; 
 241:               int o_vt2 = -1 - LDVT2 + offset_vt2; int o_idxc = -1 + offset_idxc;  int o_ctot = -1 + offset_ctot; 
 242:               int o_z = -1 + offset_z;
 243:   
 244:              #endregion
 245:   
 246:   
 247:              #region Prolog
 248:              
 249:              // *
 250:              // *  -- LAPACK auxiliary routine (version 3.1) --
 251:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 252:              // *     November 2006
 253:              // *
 254:              // *     .. Scalar Arguments ..
 255:              // *     ..
 256:              // *     .. Array Arguments ..
 257:              // *     ..
 258:              // *
 259:              // *  Purpose
 260:              // *  =======
 261:              // *
 262:              // *  DLASD3 finds all the square roots of the roots of the secular
 263:              // *  equation, as defined by the values in D and Z.  It makes the
 264:              // *  appropriate calls to DLASD4 and then updates the singular
 265:              // *  vectors by matrix multiplication.
 266:              // *
 267:              // *  This code makes very mild assumptions about floating point
 268:              // *  arithmetic. It will work on machines with a guard digit in
 269:              // *  add/subtract, or on those binary machines without guard digits
 270:              // *  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
 271:              // *  It could conceivably fail on hexadecimal or decimal machines
 272:              // *  without guard digits, but we know of none.
 273:              // *
 274:              // *  DLASD3 is called from DLASD1.
 275:              // *
 276:              // *  Arguments
 277:              // *  =========
 278:              // *
 279:              // *  NL     (input) INTEGER
 280:              // *         The row dimension of the upper block.  NL >= 1.
 281:              // *
 282:              // *  NR     (input) INTEGER
 283:              // *         The row dimension of the lower block.  NR >= 1.
 284:              // *
 285:              // *  SQRE   (input) INTEGER
 286:              // *         = 0: the lower block is an NR-by-NR square matrix.
 287:              // *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
 288:              // *
 289:              // *         The bidiagonal matrix has N = NL + NR + 1 rows and
 290:              // *         M = N + SQRE >= N columns.
 291:              // *
 292:              // *  K      (input) INTEGER
 293:              // *         The size of the secular equation, 1 =< K = < N.
 294:              // *
 295:              // *  D      (output) DOUBLE PRECISION array, dimension(K)
 296:              // *         On exit the square roots of the roots of the secular equation,
 297:              // *         in ascending order.
 298:              // *
 299:              // *  Q      (workspace) DOUBLE PRECISION array,
 300:              // *                     dimension at least (LDQ,K).
 301:              // *
 302:              // *  LDQ    (input) INTEGER
 303:              // *         The leading dimension of the array Q.  LDQ >= K.
 304:              // *
 305:              // *  DSIGMA (input) DOUBLE PRECISION array, dimension(K)
 306:              // *         The first K elements of this array contain the old roots
 307:              // *         of the deflated updating problem.  These are the poles
 308:              // *         of the secular equation.
 309:              // *
 310:              // *  U      (output) DOUBLE PRECISION array, dimension (LDU, N)
 311:              // *         The last N - K columns of this matrix contain the deflated
 312:              // *         left singular vectors.
 313:              // *
 314:              // *  LDU    (input) INTEGER
 315:              // *         The leading dimension of the array U.  LDU >= N.
 316:              // *
 317:              // *  U2     (input/output) DOUBLE PRECISION array, dimension (LDU2, N)
 318:              // *         The first K columns of this matrix contain the non-deflated
 319:              // *         left singular vectors for the split problem.
 320:              // *
 321:              // *  LDU2   (input) INTEGER
 322:              // *         The leading dimension of the array U2.  LDU2 >= N.
 323:              // *
 324:              // *  VT     (output) DOUBLE PRECISION array, dimension (LDVT, M)
 325:              // *         The last M - K columns of VT' contain the deflated
 326:              // *         right singular vectors.
 327:              // *
 328:              // *  LDVT   (input) INTEGER
 329:              // *         The leading dimension of the array VT.  LDVT >= N.
 330:              // *
 331:              // *  VT2    (input/output) DOUBLE PRECISION array, dimension (LDVT2, N)
 332:              // *         The first K columns of VT2' contain the non-deflated
 333:              // *         right singular vectors for the split problem.
 334:              // *
 335:              // *  LDVT2  (input) INTEGER
 336:              // *         The leading dimension of the array VT2.  LDVT2 >= N.
 337:              // *
 338:              // *  IDXC   (input) INTEGER array, dimension ( N )
 339:              // *         The permutation used to arrange the columns of U (and rows of
 340:              // *         VT) into three groups:  the first group contains non-zero
 341:              // *         entries only at and above (or before) NL +1; the second
 342:              // *         contains non-zero entries only at and below (or after) NL+2;
 343:              // *         and the third is dense. The first column of U and the row of
 344:              // *         VT are treated separately, however.
 345:              // *
 346:              // *         The rows of the singular vectors found by DLASD4
 347:              // *         must be likewise permuted before the matrix multiplies can
 348:              // *         take place.
 349:              // *
 350:              // *  CTOT   (input) INTEGER array, dimension ( 4 )
 351:              // *         A count of the total number of the various types of columns
 352:              // *         in U (or rows in VT), as described in IDXC. The fourth column
 353:              // *         type is any column which has been deflated.
 354:              // *
 355:              // *  Z      (input) DOUBLE PRECISION array, dimension (K)
 356:              // *         The first K elements of this array contain the components
 357:              // *         of the deflation-adjusted updating row vector.
 358:              // *
 359:              // *  INFO   (output) INTEGER
 360:              // *         = 0:  successful exit.
 361:              // *         < 0:  if INFO = -i, the i-th argument had an illegal value.
 362:              // *         > 0:  if INFO = 1, an singular value did not converge
 363:              // *
 364:              // *  Further Details
 365:              // *  ===============
 366:              // *
 367:              // *  Based on contributions by
 368:              // *     Ming Gu and Huan Ren, Computer Science Division, University of
 369:              // *     California at Berkeley, USA
 370:              // *
 371:              // *  =====================================================================
 372:              // *
 373:              // *     .. Parameters ..
 374:              // *     ..
 375:              // *     .. Local Scalars ..
 376:              // *     ..
 377:              // *     .. External Functions ..
 378:              // *     ..
 379:              // *     .. External Subroutines ..
 380:              // *     ..
 381:              // *     .. Intrinsic Functions ..
 382:              //      INTRINSIC          ABS, SIGN, SQRT;
 383:              // *     ..
 384:              // *     .. Executable Statements ..
 385:              // *
 386:              // *     Test the input parameters.
 387:              // *
 388:   
 389:              #endregion
 390:   
 391:   
 392:              #region Body
 393:              
 394:              INFO = 0;
 395:              // *
 396:              if (NL < 1)
 397:              {
 398:                  INFO =  - 1;
 399:              }
 400:              else
 401:              {
 402:                  if (NR < 1)
 403:                  {
 404:                      INFO =  - 2;
 405:                  }
 406:                  else
 407:                  {
 408:                      if ((SQRE != 1) && (SQRE != 0))
 409:                      {
 410:                          INFO =  - 3;
 411:                      }
 412:                  }
 413:              }
 414:              // *
 415:              N = NL + NR + 1;
 416:              M = N + SQRE;
 417:              NLP1 = NL + 1;
 418:              NLP2 = NL + 2;
 419:              // *
 420:              if ((K < 1) || (K > N))
 421:              {
 422:                  INFO =  - 4;
 423:              }
 424:              else
 425:              {
 426:                  if (LDQ < K)
 427:                  {
 428:                      INFO =  - 7;
 429:                  }
 430:                  else
 431:                  {
 432:                      if (LDU < N)
 433:                      {
 434:                          INFO =  - 10;
 435:                      }
 436:                      else
 437:                      {
 438:                          if (LDU2 < N)
 439:                          {
 440:                              INFO =  - 12;
 441:                          }
 442:                          else
 443:                          {
 444:                              if (LDVT < M)
 445:                              {
 446:                                  INFO =  - 14;
 447:                              }
 448:                              else
 449:                              {
 450:                                  if (LDVT2 < M)
 451:                                  {
 452:                                      INFO =  - 16;
 453:                                  }
 454:                              }
 455:                          }
 456:                      }
 457:                  }
 458:              }
 459:              if (INFO != 0)
 460:              {
 461:                  this._xerbla.Run("DLASD3",  - INFO);
 462:                  return;
 463:              }
 464:              // *
 465:              // *     Quick return if possible
 466:              // *
 467:              if (K == 1)
 468:              {
 469:                  D[1 + o_d] = Math.Abs(Z[1 + o_z]);
 470:                  this._dcopy.Run(M, VT2, 1+1 * LDVT2 + o_vt2, LDVT2, ref VT, 1+1 * LDVT + o_vt, LDVT);
 471:                  if (Z[1 + o_z] > ZERO)
 472:                  {
 473:                      this._dcopy.Run(N, U2, 1+1 * LDU2 + o_u2, 1, ref U, 1+1 * LDU + o_u, 1);
 474:                  }
 475:                  else
 476:                  {
 477:                      for (I = 1; I <= N; I++)
 478:                      {
 479:                          U[I+1 * LDU + o_u] =  - U2[I+1 * LDU2 + o_u2];
 480:                      }
 481:                  }
 482:                  return;
 483:              }
 484:              // *
 485:              // *     Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
 486:              // *     be computed with high relative accuracy (barring over/underflow).
 487:              // *     This is a problem on machines without a guard digit in
 488:              // *     add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
 489:              // *     The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
 490:              // *     which on any of these machines zeros out the bottommost
 491:              // *     bit of DSIGMA(I) if it is 1; this makes the subsequent
 492:              // *     subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
 493:              // *     occurs. On binary machines with a guard digit (almost all
 494:              // *     machines) it does not change DSIGMA(I) at all. On hexadecimal
 495:              // *     and decimal machines with a guard digit, it slightly
 496:              // *     changes the bottommost bits of DSIGMA(I). It does not account
 497:              // *     for hexadecimal or decimal machines without guard digits
 498:              // *     (we know of none). We use a subroutine call to compute
 499:              // *     2*DSIGMA(I) to prevent optimizing compilers from eliminating
 500:              // *     this code.
 501:              // *
 502:              for (I = 1; I <= K; I++)
 503:              {
 504:                  DSIGMA[I + o_dsigma] = this._dlamc3.Run(DSIGMA[I + o_dsigma], DSIGMA[I + o_dsigma]) - DSIGMA[I + o_dsigma];
 505:              }
 506:              // *
 507:              // *     Keep a copy of Z.
 508:              // *
 509:              this._dcopy.Run(K, Z, offset_z, 1, ref Q, offset_q, 1);
 510:              // *
 511:              // *     Normalize Z.
 512:              // *
 513:              RHO = this._dnrm2.Run(K, Z, offset_z, 1);
 514:              this._dlascl.Run("G", 0, 0, RHO, ONE, K
 515:                               , 1, ref Z, offset_z, K, ref INFO);
 516:              RHO = RHO * RHO;
 517:              // *
 518:              // *     Find the new singular values.
 519:              // *
 520:              for (J = 1; J <= K; J++)
 521:              {
 522:                  this._dlasd4.Run(K, J, DSIGMA, offset_dsigma, Z, offset_z, ref U, 1+J * LDU + o_u, RHO
 523:                                   , ref D[J + o_d], ref VT, 1+J * LDVT + o_vt, ref INFO);
 524:                  // *
 525:                  // *        If the zero finder fails, the computation is terminated.
 526:                  // *
 527:                  if (INFO != 0)
 528:                  {
 529:                      return;
 530:                  }
 531:              }
 532:              // *
 533:              // *     Compute updated Z.
 534:              // *
 535:              for (I = 1; I <= K; I++)
 536:              {
 537:                  Z[I + o_z] = U[I+K * LDU + o_u] * VT[I+K * LDVT + o_vt];
 538:                  for (J = 1; J <= I - 1; J++)
 539:                  {
 540:                      Z[I + o_z] = Z[I + o_z] * (U[I+J * LDU + o_u] * VT[I+J * LDVT + o_vt] / (DSIGMA[I + o_dsigma] - DSIGMA[J + o_dsigma]) / (DSIGMA[I + o_dsigma] + DSIGMA[J + o_dsigma]));
 541:                  }
 542:                  for (J = I; J <= K - 1; J++)
 543:                  {
 544:                      Z[I + o_z] = Z[I + o_z] * (U[I+J * LDU + o_u] * VT[I+J * LDVT + o_vt] / (DSIGMA[I + o_dsigma] - DSIGMA[J + 1 + o_dsigma]) / (DSIGMA[I + o_dsigma] + DSIGMA[J + 1 + o_dsigma]));
 545:                  }
 546:                  Z[I + o_z] = FortranLib.Sign(Math.Sqrt(Math.Abs(Z[I + o_z])),Q[I+1 * LDQ + o_q]);
 547:              }
 548:              // *
 549:              // *     Compute left singular vectors of the modified diagonal matrix,
 550:              // *     and store related information for the right singular vectors.
 551:              // *
 552:              for (I = 1; I <= K; I++)
 553:              {
 554:                  VT[1+I * LDVT + o_vt] = Z[1 + o_z] / U[1+I * LDU + o_u] / VT[1+I * LDVT + o_vt];
 555:                  U[1+I * LDU + o_u] = NEGONE;
 556:                  for (J = 2; J <= K; J++)
 557:                  {
 558:                      VT[J+I * LDVT + o_vt] = Z[J + o_z] / U[J+I * LDU + o_u] / VT[J+I * LDVT + o_vt];
 559:                      U[J+I * LDU + o_u] = DSIGMA[J + o_dsigma] * VT[J+I * LDVT + o_vt];
 560:                  }
 561:                  TEMP = this._dnrm2.Run(K, U, 1+I * LDU + o_u, 1);
 562:                  Q[1+I * LDQ + o_q] = U[1+I * LDU + o_u] / TEMP;
 563:                  for (J = 2; J <= K; J++)
 564:                  {
 565:                      JC = IDXC[J + o_idxc];
 566:                      Q[J+I * LDQ + o_q] = U[JC+I * LDU + o_u] / TEMP;
 567:                  }
 568:              }
 569:              // *
 570:              // *     Update the left singular vector matrix.
 571:              // *
 572:              if (K == 2)
 573:              {
 574:                  this._dgemm.Run("N", "N", N, K, K, ONE
 575:                                  , U2, offset_u2, LDU2, Q, offset_q, LDQ, ZERO, ref U, offset_u
 576:                                  , LDU);
 577:                  goto LABEL100;
 578:              }
 579:              if (CTOT[1 + o_ctot] > 0)
 580:              {
 581:                  this._dgemm.Run("N", "N", NL, K, CTOT[1 + o_ctot], ONE
 582:                                  , U2, 1+2 * LDU2 + o_u2, LDU2, Q, 2+1 * LDQ + o_q, LDQ, ZERO, ref U, 1+1 * LDU + o_u
 583:                                  , LDU);
 584:                  if (CTOT[3 + o_ctot] > 0)
 585:                  {
 586:                      KTEMP = 2 + CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
 587:                      this._dgemm.Run("N", "N", NL, K, CTOT[3 + o_ctot], ONE
 588:                                      , U2, 1+KTEMP * LDU2 + o_u2, LDU2, Q, KTEMP+1 * LDQ + o_q, LDQ, ONE, ref U, 1+1 * LDU + o_u
 589:                                      , LDU);
 590:                  }
 591:              }
 592:              else
 593:              {
 594:                  if (CTOT[3 + o_ctot] > 0)
 595:                  {
 596:                      KTEMP = 2 + CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
 597:                      this._dgemm.Run("N", "N", NL, K, CTOT[3 + o_ctot], ONE
 598:                                      , U2, 1+KTEMP * LDU2 + o_u2, LDU2, Q, KTEMP+1 * LDQ + o_q, LDQ, ZERO, ref U, 1+1 * LDU + o_u
 599:                                      , LDU);
 600:                  }
 601:                  else
 602:                  {
 603:                      this._dlacpy.Run("F", NL, K, U2, offset_u2, LDU2, ref U, offset_u
 604:                                       , LDU);
 605:                  }
 606:              }
 607:              this._dcopy.Run(K, Q, 1+1 * LDQ + o_q, LDQ, ref U, NLP1+1 * LDU + o_u, LDU);
 608:              KTEMP = 2 + CTOT[1 + o_ctot];
 609:              CTEMP = CTOT[2 + o_ctot] + CTOT[3 + o_ctot];
 610:              this._dgemm.Run("N", "N", NR, K, CTEMP, ONE
 611:                              , U2, NLP2+KTEMP * LDU2 + o_u2, LDU2, Q, KTEMP+1 * LDQ + o_q, LDQ, ZERO, ref U, NLP2+1 * LDU + o_u
 612:                              , LDU);
 613:              // *
 614:              // *     Generate the right singular vectors.
 615:              // *
 616:          LABEL100:;
 617:              for (I = 1; I <= K; I++)
 618:              {
 619:                  TEMP = this._dnrm2.Run(K, VT, 1+I * LDVT + o_vt, 1);
 620:                  Q[I+1 * LDQ + o_q] = VT[1+I * LDVT + o_vt] / TEMP;
 621:                  for (J = 2; J <= K; J++)
 622:                  {
 623:                      JC = IDXC[J + o_idxc];
 624:                      Q[I+J * LDQ + o_q] = VT[JC+I * LDVT + o_vt] / TEMP;
 625:                  }
 626:              }
 627:              // *
 628:              // *     Update the right singular vector matrix.
 629:              // *
 630:              if (K == 2)
 631:              {
 632:                  this._dgemm.Run("N", "N", K, M, K, ONE
 633:                                  , Q, offset_q, LDQ, VT2, offset_vt2, LDVT2, ZERO, ref VT, offset_vt
 634:                                  , LDVT);
 635:                  return;
 636:              }
 637:              KTEMP = 1 + CTOT[1 + o_ctot];
 638:              this._dgemm.Run("N", "N", K, NLP1, KTEMP, ONE
 639:                              , Q, 1+1 * LDQ + o_q, LDQ, VT2, 1+1 * LDVT2 + o_vt2, LDVT2, ZERO, ref VT, 1+1 * LDVT + o_vt
 640:                              , LDVT);
 641:              KTEMP = 2 + CTOT[1 + o_ctot] + CTOT[2 + o_ctot];
 642:              if (KTEMP <= LDVT2)
 643:              {
 644:                  this._dgemm.Run("N", "N", K, NLP1, CTOT[3 + o_ctot], ONE
 645:                                  , Q, 1+KTEMP * LDQ + o_q, LDQ, VT2, KTEMP+1 * LDVT2 + o_vt2, LDVT2, ONE, ref VT, 1+1 * LDVT + o_vt
 646:                                  , LDVT);
 647:              }
 648:              // *
 649:              KTEMP = CTOT[1 + o_ctot] + 1;
 650:              NRP1 = NR + SQRE;
 651:              if (KTEMP > 1)
 652:              {
 653:                  for (I = 1; I <= K; I++)
 654:                  {
 655:                      Q[I+KTEMP * LDQ + o_q] = Q[I+1 * LDQ + o_q];
 656:                  }
 657:                  for (I = NLP2; I <= M; I++)
 658:                  {
 659:                      VT2[KTEMP+I * LDVT2 + o_vt2] = VT2[1+I * LDVT2 + o_vt2];
 660:                  }
 661:              }
 662:              CTEMP = 1 + CTOT[2 + o_ctot] + CTOT[3 + o_ctot];
 663:              this._dgemm.Run("N", "N", K, NRP1, CTEMP, ONE
 664:                              , Q, 1+KTEMP * LDQ + o_q, LDQ, VT2, KTEMP+NLP2 * LDVT2 + o_vt2, LDVT2, ZERO, ref VT, 1+NLP2 * LDVT + o_vt
 665:                              , LDVT);
 666:              // *
 667:              return;
 668:              // *
 669:              // *     End of DLASD3
 670:              // *
 671:   
 672:              #endregion
 673:   
 674:          }
 675:      }
 676:  }