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:      /// DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
  27:      /// where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
  28:      /// 
  29:      /// A related subroutine DLASD7 handles the case in which the singular
  30:      /// values (and the singular vectors in factored form) are desired.
  31:      /// 
  32:      /// DLASD1 computes the SVD as follows:
  33:      /// 
  34:      /// ( D1(in)  0    0     0 )
  35:      /// B = U(in) * (   Z1'   a   Z2'    b ) * VT(in)
  36:      /// (   0     0   D2(in) 0 )
  37:      /// 
  38:      /// = U(out) * ( D(out) 0) * VT(out)
  39:      /// 
  40:      /// where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
  41:      /// with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
  42:      /// elsewhere; and the entry b is empty if SQRE = 0.
  43:      /// 
  44:      /// The left singular vectors of the original matrix are stored in U, and
  45:      /// the transpose of the right singular vectors are stored in VT, and the
  46:      /// singular values are in D.  The algorithm consists of three stages:
  47:      /// 
  48:      /// The first stage consists of deflating the size of the problem
  49:      /// when there are multiple singular values or when there are zeros in
  50:      /// the Z vector.  For each such occurence the dimension of the
  51:      /// secular equation problem is reduced by one.  This stage is
  52:      /// performed by the routine DLASD2.
  53:      /// 
  54:      /// The second stage consists of calculating the updated
  55:      /// singular values. This is done by finding the square roots of the
  56:      /// roots of the secular equation via the routine DLASD4 (as called
  57:      /// by DLASD3). This routine also calculates the singular vectors of
  58:      /// the current problem.
  59:      /// 
  60:      /// The final stage consists of computing the updated singular vectors
  61:      /// directly using the updated singular values.  The singular vectors
  62:      /// for the current problem are multiplied with the singular vectors
  63:      /// from the overall problem.
  64:      /// 
  65:      ///</summary>
  66:      public class DLASD1
  67:      {
  68:      
  69:   
  70:          #region Dependencies
  71:          
  72:          DLAMRG _dlamrg; DLASCL _dlascl; DLASD2 _dlasd2; DLASD3 _dlasd3; XERBLA _xerbla; 
  73:   
  74:          #endregion
  75:   
  76:   
  77:          #region Fields
  78:          
  79:          const double ONE = 1.0E+0; const double ZERO = 0.0E+0; int COLTYP = 0; int I = 0; int IDX = 0; int IDXC = 0; int IDXP = 0; 
  80:          int IQ = 0;int ISIGMA = 0; int IU2 = 0; int IVT2 = 0; int IZ = 0; int K = 0; int LDQ = 0; int LDU2 = 0; int LDVT2 = 0; 
  81:          int M = 0;int N = 0; int N1 = 0; int N2 = 0; double ORGNRM = 0; 
  82:   
  83:          #endregion
  84:   
  85:          public DLASD1(DLAMRG dlamrg, DLASCL dlascl, DLASD2 dlasd2, DLASD3 dlasd3, XERBLA xerbla)
  86:          {
  87:      
  88:   
  89:              #region Set Dependencies
  90:              
  91:              this._dlamrg = dlamrg; this._dlascl = dlascl; this._dlasd2 = dlasd2; this._dlasd3 = dlasd3; this._xerbla = xerbla; 
  92:   
  93:              #endregion
  94:   
  95:          }
  96:      
  97:          public DLASD1()
  98:          {
  99:      
 100:   
 101:              #region Dependencies (Initialization)
 102:              
 103:              DLAMRG dlamrg = new DLAMRG();
 104:              LSAME lsame = new LSAME();
 105:              DLAMC3 dlamc3 = new DLAMC3();
 106:              XERBLA xerbla = new XERBLA();
 107:              DLAPY2 dlapy2 = new DLAPY2();
 108:              DCOPY dcopy = new DCOPY();
 109:              DROT drot = new DROT();
 110:              DNRM2 dnrm2 = new DNRM2();
 111:              DLASD5 dlasd5 = new DLASD5();
 112:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 113:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 114:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 115:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 116:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 117:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 118:              DLACPY dlacpy = new DLACPY(lsame);
 119:              DLASET dlaset = new DLASET(lsame);
 120:              DLASD2 dlasd2 = new DLASD2(dlamch, dlapy2, dcopy, dlacpy, dlamrg, dlaset, drot, xerbla);
 121:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 122:              DLAED6 dlaed6 = new DLAED6(dlamch);
 123:              DLASD4 dlasd4 = new DLASD4(dlaed6, dlasd5, dlamch);
 124:              DLASD3 dlasd3 = new DLASD3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlascl, dlasd4, xerbla);
 125:   
 126:              #endregion
 127:   
 128:   
 129:              #region Set Dependencies
 130:              
 131:              this._dlamrg = dlamrg; this._dlascl = dlascl; this._dlasd2 = dlasd2; this._dlasd3 = dlasd3; this._xerbla = xerbla; 
 132:   
 133:              #endregion
 134:   
 135:          }
 136:          /// <summary>
 137:          /// Purpose
 138:          /// =======
 139:          /// 
 140:          /// DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
 141:          /// where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
 142:          /// 
 143:          /// A related subroutine DLASD7 handles the case in which the singular
 144:          /// values (and the singular vectors in factored form) are desired.
 145:          /// 
 146:          /// DLASD1 computes the SVD as follows:
 147:          /// 
 148:          /// ( D1(in)  0    0     0 )
 149:          /// B = U(in) * (   Z1'   a   Z2'    b ) * VT(in)
 150:          /// (   0     0   D2(in) 0 )
 151:          /// 
 152:          /// = U(out) * ( D(out) 0) * VT(out)
 153:          /// 
 154:          /// where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
 155:          /// with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
 156:          /// elsewhere; and the entry b is empty if SQRE = 0.
 157:          /// 
 158:          /// The left singular vectors of the original matrix are stored in U, and
 159:          /// the transpose of the right singular vectors are stored in VT, and the
 160:          /// singular values are in D.  The algorithm consists of three stages:
 161:          /// 
 162:          /// The first stage consists of deflating the size of the problem
 163:          /// when there are multiple singular values or when there are zeros in
 164:          /// the Z vector.  For each such occurence the dimension of the
 165:          /// secular equation problem is reduced by one.  This stage is
 166:          /// performed by the routine DLASD2.
 167:          /// 
 168:          /// The second stage consists of calculating the updated
 169:          /// singular values. This is done by finding the square roots of the
 170:          /// roots of the secular equation via the routine DLASD4 (as called
 171:          /// by DLASD3). This routine also calculates the singular vectors of
 172:          /// the current problem.
 173:          /// 
 174:          /// The final stage consists of computing the updated singular vectors
 175:          /// directly using the updated singular values.  The singular vectors
 176:          /// for the current problem are multiplied with the singular vectors
 177:          /// from the overall problem.
 178:          /// 
 179:          ///</summary>
 180:          /// <param name="NL">
 181:          /// (input) INTEGER
 182:          /// The row dimension of the upper block.  NL .GE. 1.
 183:          ///</param>
 184:          /// <param name="NR">
 185:          /// (input) INTEGER
 186:          /// The row dimension of the lower block.  NR .GE. 1.
 187:          ///</param>
 188:          /// <param name="SQRE">
 189:          /// (input) INTEGER
 190:          /// = 0: the lower block is an NR-by-NR square matrix.
 191:          /// = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
 192:          /// 
 193:          /// The bidiagonal matrix has row dimension N = NL + NR + 1,
 194:          /// and column dimension M = N + SQRE.
 195:          ///</param>
 196:          /// <param name="D">
 197:          /// (input/output) DOUBLE PRECISION array,
 198:          /// dimension (N = NL+NR+1).
 199:          /// On entry D(1:NL,1:NL) contains the singular values of the
 200:          /// upper block; and D(NL+2:N) contains the singular values of
 201:          /// the lower block. On exit D(1:N) contains the singular values
 202:          /// of the modified matrix.
 203:          ///</param>
 204:          /// <param name="ALPHA">
 205:          /// (input/output) DOUBLE PRECISION
 206:          /// Contains the diagonal element associated with the added row.
 207:          ///</param>
 208:          /// <param name="BETA">
 209:          /// (input/output) DOUBLE PRECISION
 210:          /// Contains the off-diagonal element associated with the added
 211:          /// row.
 212:          ///</param>
 213:          /// <param name="U">
 214:          /// (input/output) DOUBLE PRECISION array, dimension(LDU,N)
 215:          /// On entry U(1:NL, 1:NL) contains the left singular vectors of
 216:          /// the upper block; U(NL+2:N, NL+2:N) contains the left singular
 217:          /// vectors of the lower block. On exit U contains the left
 218:          /// singular vectors of the bidiagonal matrix.
 219:          ///</param>
 220:          /// <param name="LDU">
 221:          /// (input) INTEGER
 222:          /// The leading dimension of the array U.  LDU .GE. max( 1, N ).
 223:          ///</param>
 224:          /// <param name="VT">
 225:          /// (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
 226:          /// where M = N + SQRE.
 227:          /// On entry VT(1:NL+1, 1:NL+1)' contains the right singular
 228:          /// vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
 229:          /// the right singular vectors of the lower block. On exit
 230:          /// VT' contains the right singular vectors of the
 231:          /// bidiagonal matrix.
 232:          ///</param>
 233:          /// <param name="LDVT">
 234:          /// (input) INTEGER
 235:          /// The leading dimension of the array VT.  LDVT .GE. max( 1, M ).
 236:          ///</param>
 237:          /// <param name="IDXQ">
 238:          /// (output) INTEGER array, dimension(N)
 239:          /// This contains the permutation which will reintegrate the
 240:          /// subproblem just solved back into sorted order, i.e.
 241:          /// D( IDXQ( I = 1, N ) ) will be in ascending order.
 242:          ///</param>
 243:          /// <param name="IWORK">
 244:          /// (workspace) INTEGER array, dimension( 4 * N )
 245:          ///</param>
 246:          /// <param name="WORK">
 247:          /// (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
 248:          ///</param>
 249:          /// <param name="INFO">
 250:          /// (output) INTEGER
 251:          /// = 0:  successful exit.
 252:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 253:          /// .GT. 0:  if INFO = 1, an singular value did not converge
 254:          ///</param>
 255:          public void Run(int NL, int NR, int SQRE, ref double[] D, int offset_d, ref double ALPHA, ref double BETA
 256:                           , ref double[] U, int offset_u, int LDU, ref double[] VT, int offset_vt, int LDVT, ref int[] IDXQ, int offset_idxq, ref int[] IWORK, int offset_iwork
 257:                           , ref double[] WORK, int offset_work, ref int INFO)
 258:          {
 259:   
 260:              #region Array Index Correction
 261:              
 262:               int o_d = -1 + offset_d;  int o_u = -1 - LDU + offset_u;  int o_vt = -1 - LDVT + offset_vt; 
 263:               int o_idxq = -1 + offset_idxq; int o_iwork = -1 + offset_iwork;  int o_work = -1 + offset_work; 
 264:   
 265:              #endregion
 266:   
 267:   
 268:              #region Prolog
 269:              
 270:              // *
 271:              // *  -- LAPACK auxiliary routine (version 3.1) --
 272:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 273:              // *     November 2006
 274:              // *
 275:              // *     .. Scalar Arguments ..
 276:              // *     ..
 277:              // *     .. Array Arguments ..
 278:              // *     ..
 279:              // *
 280:              // *  Purpose
 281:              // *  =======
 282:              // *
 283:              // *  DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
 284:              // *  where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
 285:              // *
 286:              // *  A related subroutine DLASD7 handles the case in which the singular
 287:              // *  values (and the singular vectors in factored form) are desired.
 288:              // *
 289:              // *  DLASD1 computes the SVD as follows:
 290:              // *
 291:              // *                ( D1(in)  0    0     0 )
 292:              // *    B = U(in) * (   Z1'   a   Z2'    b ) * VT(in)
 293:              // *                (   0     0   D2(in) 0 )
 294:              // *
 295:              // *      = U(out) * ( D(out) 0) * VT(out)
 296:              // *
 297:              // *  where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
 298:              // *  with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
 299:              // *  elsewhere; and the entry b is empty if SQRE = 0.
 300:              // *
 301:              // *  The left singular vectors of the original matrix are stored in U, and
 302:              // *  the transpose of the right singular vectors are stored in VT, and the
 303:              // *  singular values are in D.  The algorithm consists of three stages:
 304:              // *
 305:              // *     The first stage consists of deflating the size of the problem
 306:              // *     when there are multiple singular values or when there are zeros in
 307:              // *     the Z vector.  For each such occurence the dimension of the
 308:              // *     secular equation problem is reduced by one.  This stage is
 309:              // *     performed by the routine DLASD2.
 310:              // *
 311:              // *     The second stage consists of calculating the updated
 312:              // *     singular values. This is done by finding the square roots of the
 313:              // *     roots of the secular equation via the routine DLASD4 (as called
 314:              // *     by DLASD3). This routine also calculates the singular vectors of
 315:              // *     the current problem.
 316:              // *
 317:              // *     The final stage consists of computing the updated singular vectors
 318:              // *     directly using the updated singular values.  The singular vectors
 319:              // *     for the current problem are multiplied with the singular vectors
 320:              // *     from the overall problem.
 321:              // *
 322:              // *  Arguments
 323:              // *  =========
 324:              // *
 325:              // *  NL     (input) INTEGER
 326:              // *         The row dimension of the upper block.  NL >= 1.
 327:              // *
 328:              // *  NR     (input) INTEGER
 329:              // *         The row dimension of the lower block.  NR >= 1.
 330:              // *
 331:              // *  SQRE   (input) INTEGER
 332:              // *         = 0: the lower block is an NR-by-NR square matrix.
 333:              // *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
 334:              // *
 335:              // *         The bidiagonal matrix has row dimension N = NL + NR + 1,
 336:              // *         and column dimension M = N + SQRE.
 337:              // *
 338:              // *  D      (input/output) DOUBLE PRECISION array,
 339:              // *                        dimension (N = NL+NR+1).
 340:              // *         On entry D(1:NL,1:NL) contains the singular values of the
 341:              // *         upper block; and D(NL+2:N) contains the singular values of
 342:              // *         the lower block. On exit D(1:N) contains the singular values
 343:              // *         of the modified matrix.
 344:              // *
 345:              // *  ALPHA  (input/output) DOUBLE PRECISION
 346:              // *         Contains the diagonal element associated with the added row.
 347:              // *
 348:              // *  BETA   (input/output) DOUBLE PRECISION
 349:              // *         Contains the off-diagonal element associated with the added
 350:              // *         row.
 351:              // *
 352:              // *  U      (input/output) DOUBLE PRECISION array, dimension(LDU,N)
 353:              // *         On entry U(1:NL, 1:NL) contains the left singular vectors of
 354:              // *         the upper block; U(NL+2:N, NL+2:N) contains the left singular
 355:              // *         vectors of the lower block. On exit U contains the left
 356:              // *         singular vectors of the bidiagonal matrix.
 357:              // *
 358:              // *  LDU    (input) INTEGER
 359:              // *         The leading dimension of the array U.  LDU >= max( 1, N ).
 360:              // *
 361:              // *  VT     (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
 362:              // *         where M = N + SQRE.
 363:              // *         On entry VT(1:NL+1, 1:NL+1)' contains the right singular
 364:              // *         vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
 365:              // *         the right singular vectors of the lower block. On exit
 366:              // *         VT' contains the right singular vectors of the
 367:              // *         bidiagonal matrix.
 368:              // *
 369:              // *  LDVT   (input) INTEGER
 370:              // *         The leading dimension of the array VT.  LDVT >= max( 1, M ).
 371:              // *
 372:              // *  IDXQ  (output) INTEGER array, dimension(N)
 373:              // *         This contains the permutation which will reintegrate the
 374:              // *         subproblem just solved back into sorted order, i.e.
 375:              // *         D( IDXQ( I = 1, N ) ) will be in ascending order.
 376:              // *
 377:              // *  IWORK  (workspace) INTEGER array, dimension( 4 * N )
 378:              // *
 379:              // *  WORK   (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
 380:              // *
 381:              // *  INFO   (output) INTEGER
 382:              // *          = 0:  successful exit.
 383:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 384:              // *          > 0:  if INFO = 1, an singular value did not converge
 385:              // *
 386:              // *  Further Details
 387:              // *  ===============
 388:              // *
 389:              // *  Based on contributions by
 390:              // *     Ming Gu and Huan Ren, Computer Science Division, University of
 391:              // *     California at Berkeley, USA
 392:              // *
 393:              // *  =====================================================================
 394:              // *
 395:              // *     .. Parameters ..
 396:              // *
 397:              // *     ..
 398:              // *     .. Local Scalars ..
 399:              // *     ..
 400:              // *     .. External Subroutines ..
 401:              // *     ..
 402:              // *     .. Intrinsic Functions ..
 403:              //      INTRINSIC          ABS, MAX;
 404:              // *     ..
 405:              // *     .. Executable Statements ..
 406:              // *
 407:              // *     Test the input parameters.
 408:              // *
 409:   
 410:              #endregion
 411:   
 412:   
 413:              #region Body
 414:              
 415:              INFO = 0;
 416:              // *
 417:              if (NL < 1)
 418:              {
 419:                  INFO =  - 1;
 420:              }
 421:              else
 422:              {
 423:                  if (NR < 1)
 424:                  {
 425:                      INFO =  - 2;
 426:                  }
 427:                  else
 428:                  {
 429:                      if ((SQRE < 0) || (SQRE > 1))
 430:                      {
 431:                          INFO =  - 3;
 432:                      }
 433:                  }
 434:              }
 435:              if (INFO != 0)
 436:              {
 437:                  this._xerbla.Run("DLASD1",  - INFO);
 438:                  return;
 439:              }
 440:              // *
 441:              N = NL + NR + 1;
 442:              M = N + SQRE;
 443:              // *
 444:              // *     The following values are for bookkeeping purposes only.  They are
 445:              // *     integer pointers which indicate the portion of the workspace
 446:              // *     used by a particular array in DLASD2 and DLASD3.
 447:              // *
 448:              LDU2 = N;
 449:              LDVT2 = M;
 450:              // *
 451:              IZ = 1;
 452:              ISIGMA = IZ + M;
 453:              IU2 = ISIGMA + N;
 454:              IVT2 = IU2 + LDU2 * N;
 455:              IQ = IVT2 + LDVT2 * M;
 456:              // *
 457:              IDX = 1;
 458:              IDXC = IDX + N;
 459:              COLTYP = IDXC + N;
 460:              IDXP = COLTYP + N;
 461:              // *
 462:              // *     Scale.
 463:              // *
 464:              ORGNRM = Math.Max(Math.Abs(ALPHA), Math.Abs(BETA));
 465:              D[NL + 1 + o_d] = ZERO;
 466:              for (I = 1; I <= N; I++)
 467:              {
 468:                  if (Math.Abs(D[I + o_d]) > ORGNRM)
 469:                  {
 470:                      ORGNRM = Math.Abs(D[I + o_d]);
 471:                  }
 472:              }
 473:              this._dlascl.Run("G", 0, 0, ORGNRM, ONE, N
 474:                               , 1, ref D, offset_d, N, ref INFO);
 475:              ALPHA = ALPHA / ORGNRM;
 476:              BETA = BETA / ORGNRM;
 477:              // *
 478:              // *     Deflate singular values.
 479:              // *
 480:              this._dlasd2.Run(NL, NR, SQRE, ref K, ref D, offset_d, ref WORK, IZ + o_work
 481:                               , ALPHA, BETA, ref U, offset_u, LDU, ref VT, offset_vt, LDVT
 482:                               , ref WORK, ISIGMA + o_work, ref WORK, IU2 + o_work, LDU2, ref WORK, IVT2 + o_work, LDVT2, ref IWORK, IDXP + o_iwork
 483:                               , ref IWORK, IDX + o_iwork, ref IWORK, IDXC + o_iwork, ref IDXQ, offset_idxq, ref IWORK, COLTYP + o_iwork, ref INFO);
 484:              // *
 485:              // *     Solve Secular Equation and update singular vectors.
 486:              // *
 487:              LDQ = K;
 488:              this._dlasd3.Run(NL, NR, SQRE, K, ref D, offset_d, ref WORK, IQ + o_work
 489:                               , LDQ, ref WORK, ISIGMA + o_work, ref U, offset_u, LDU, WORK, IU2 + o_work, LDU2
 490:                               , ref VT, offset_vt, LDVT, ref WORK, IVT2 + o_work, LDVT2, IWORK, IDXC + o_iwork, IWORK, COLTYP + o_iwork
 491:                               , ref WORK, IZ + o_work, ref INFO);
 492:              if (INFO != 0)
 493:              {
 494:                  return;
 495:              }
 496:              // *
 497:              // *     Unscale.
 498:              // *
 499:              this._dlascl.Run("G", 0, 0, ONE, ORGNRM, N
 500:                               , 1, ref D, offset_d, N, ref INFO);
 501:              // *
 502:              // *     Prepare the IDXQ sorting permutation.
 503:              // *
 504:              N1 = K;
 505:              N2 = N - K;
 506:              this._dlamrg.Run(N1, N2, D, offset_d, 1,  - 1, ref IDXQ, offset_idxq);
 507:              // *
 508:              return;
 509:              // *
 510:              // *     End of DLASD1
 511:              // *
 512:   
 513:              #endregion
 514:   
 515:          }
 516:      }
 517:  }