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) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DGEBD2 reduces a real general m by n matrix A to upper or lower
  27:      /// bidiagonal form B by an orthogonal transformation: Q' * A * P = B.
  28:      /// 
  29:      /// If m .GE. n, B is upper bidiagonal; if m .LT. n, B is lower bidiagonal.
  30:      /// 
  31:      ///</summary>
  32:      public class DGEBD2
  33:      {
  34:      
  35:   
  36:          #region Dependencies
  37:          
  38:          DLARF _dlarf; DLARFG _dlarfg; XERBLA _xerbla; 
  39:   
  40:          #endregion
  41:   
  42:   
  43:          #region Fields
  44:          
  45:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; 
  46:   
  47:          #endregion
  48:   
  49:          public DGEBD2(DLARF dlarf, DLARFG dlarfg, XERBLA xerbla)
  50:          {
  51:      
  52:   
  53:              #region Set Dependencies
  54:              
  55:              this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla; 
  56:   
  57:              #endregion
  58:   
  59:          }
  60:      
  61:          public DGEBD2()
  62:          {
  63:      
  64:   
  65:              #region Dependencies (Initialization)
  66:              
  67:              LSAME lsame = new LSAME();
  68:              XERBLA xerbla = new XERBLA();
  69:              DLAMC3 dlamc3 = new DLAMC3();
  70:              DLAPY2 dlapy2 = new DLAPY2();
  71:              DNRM2 dnrm2 = new DNRM2();
  72:              DSCAL dscal = new DSCAL();
  73:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  74:              DGER dger = new DGER(xerbla);
  75:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  76:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  77:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  78:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  79:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  80:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  81:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  82:   
  83:              #endregion
  84:   
  85:   
  86:              #region Set Dependencies
  87:              
  88:              this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla; 
  89:   
  90:              #endregion
  91:   
  92:          }
  93:          /// <summary>
  94:          /// Purpose
  95:          /// =======
  96:          /// 
  97:          /// DGEBD2 reduces a real general m by n matrix A to upper or lower
  98:          /// bidiagonal form B by an orthogonal transformation: Q' * A * P = B.
  99:          /// 
 100:          /// If m .GE. n, B is upper bidiagonal; if m .LT. n, B is lower bidiagonal.
 101:          /// 
 102:          ///</summary>
 103:          /// <param name="M">
 104:          /// (input) INTEGER
 105:          /// The number of rows in the matrix A.  M .GE. 0.
 106:          ///</param>
 107:          /// <param name="N">
 108:          /// (input) INTEGER
 109:          /// The number of columns in the matrix A.  N .GE. 0.
 110:          ///</param>
 111:          /// <param name="A">
 112:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 113:          /// On entry, the m by n general matrix to be reduced.
 114:          /// On exit,
 115:          /// if m .GE. n, the diagonal and the first superdiagonal are
 116:          /// overwritten with the upper bidiagonal matrix B; the
 117:          /// elements below the diagonal, with the array TAUQ, represent
 118:          /// the orthogonal matrix Q as a product of elementary
 119:          /// reflectors, and the elements above the first superdiagonal,
 120:          /// with the array TAUP, represent the orthogonal matrix P as
 121:          /// a product of elementary reflectors;
 122:          /// if m .LT. n, the diagonal and the first subdiagonal are
 123:          /// overwritten with the lower bidiagonal matrix B; the
 124:          /// elements below the first subdiagonal, with the array TAUQ,
 125:          /// represent the orthogonal matrix Q as a product of
 126:          /// elementary reflectors, and the elements above the diagonal,
 127:          /// with the array TAUP, represent the orthogonal matrix P as
 128:          /// a product of elementary reflectors.
 129:          /// See Further Details.
 130:          ///</param>
 131:          /// <param name="LDA">
 132:          /// (input) INTEGER
 133:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 134:          ///</param>
 135:          /// <param name="D">
 136:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 137:          /// The diagonal elements of the bidiagonal matrix B:
 138:          /// D(i) = A(i,i).
 139:          ///</param>
 140:          /// <param name="E">
 141:          /// (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
 142:          /// The off-diagonal elements of the bidiagonal matrix B:
 143:          /// if m .GE. n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
 144:          /// if m .LT. n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
 145:          ///</param>
 146:          /// <param name="TAUQ">
 147:          /// (output) DOUBLE PRECISION array dimension (min(M,N))
 148:          /// The scalar factors of the elementary reflectors which
 149:          /// represent the orthogonal matrix Q. See Further Details.
 150:          ///</param>
 151:          /// <param name="TAUP">
 152:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 153:          /// The scalar factors of the elementary reflectors which
 154:          /// represent the orthogonal matrix P. See Further Details.
 155:          ///</param>
 156:          /// <param name="WORK">
 157:          /// (workspace) DOUBLE PRECISION array, dimension (max(M,N))
 158:          ///</param>
 159:          /// <param name="INFO">
 160:          /// (output) INTEGER
 161:          /// = 0: successful exit.
 162:          /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
 163:          ///</param>
 164:          public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref double[] D, int offset_d, ref double[] E, int offset_e
 165:                           , ref double[] TAUQ, int offset_tauq, ref double[] TAUP, int offset_taup, ref double[] WORK, int offset_work, ref int INFO)
 166:          {
 167:   
 168:              #region Array Index Correction
 169:              
 170:               int o_a = -1 - LDA + offset_a;  int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_tauq = -1 + offset_tauq; 
 171:               int o_taup = -1 + offset_taup; int o_work = -1 + offset_work; 
 172:   
 173:              #endregion
 174:   
 175:   
 176:              #region Prolog
 177:              
 178:              // *
 179:              // *  -- LAPACK routine (version 3.1) --
 180:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 181:              // *     November 2006
 182:              // *
 183:              // *     .. Scalar Arguments ..
 184:              // *     ..
 185:              // *     .. Array Arguments ..
 186:              // *     ..
 187:              // *
 188:              // *  Purpose
 189:              // *  =======
 190:              // *
 191:              // *  DGEBD2 reduces a real general m by n matrix A to upper or lower
 192:              // *  bidiagonal form B by an orthogonal transformation: Q' * A * P = B.
 193:              // *
 194:              // *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
 195:              // *
 196:              // *  Arguments
 197:              // *  =========
 198:              // *
 199:              // *  M       (input) INTEGER
 200:              // *          The number of rows in the matrix A.  M >= 0.
 201:              // *
 202:              // *  N       (input) INTEGER
 203:              // *          The number of columns in the matrix A.  N >= 0.
 204:              // *
 205:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 206:              // *          On entry, the m by n general matrix to be reduced.
 207:              // *          On exit,
 208:              // *          if m >= n, the diagonal and the first superdiagonal are
 209:              // *            overwritten with the upper bidiagonal matrix B; the
 210:              // *            elements below the diagonal, with the array TAUQ, represent
 211:              // *            the orthogonal matrix Q as a product of elementary
 212:              // *            reflectors, and the elements above the first superdiagonal,
 213:              // *            with the array TAUP, represent the orthogonal matrix P as
 214:              // *            a product of elementary reflectors;
 215:              // *          if m < n, the diagonal and the first subdiagonal are
 216:              // *            overwritten with the lower bidiagonal matrix B; the
 217:              // *            elements below the first subdiagonal, with the array TAUQ,
 218:              // *            represent the orthogonal matrix Q as a product of
 219:              // *            elementary reflectors, and the elements above the diagonal,
 220:              // *            with the array TAUP, represent the orthogonal matrix P as
 221:              // *            a product of elementary reflectors.
 222:              // *          See Further Details.
 223:              // *
 224:              // *  LDA     (input) INTEGER
 225:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 226:              // *
 227:              // *  D       (output) DOUBLE PRECISION array, dimension (min(M,N))
 228:              // *          The diagonal elements of the bidiagonal matrix B:
 229:              // *          D(i) = A(i,i).
 230:              // *
 231:              // *  E       (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
 232:              // *          The off-diagonal elements of the bidiagonal matrix B:
 233:              // *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
 234:              // *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
 235:              // *
 236:              // *  TAUQ    (output) DOUBLE PRECISION array dimension (min(M,N))
 237:              // *          The scalar factors of the elementary reflectors which
 238:              // *          represent the orthogonal matrix Q. See Further Details.
 239:              // *
 240:              // *  TAUP    (output) DOUBLE PRECISION array, dimension (min(M,N))
 241:              // *          The scalar factors of the elementary reflectors which
 242:              // *          represent the orthogonal matrix P. See Further Details.
 243:              // *
 244:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(M,N))
 245:              // *
 246:              // *  INFO    (output) INTEGER
 247:              // *          = 0: successful exit.
 248:              // *          < 0: if INFO = -i, the i-th argument had an illegal value.
 249:              // *
 250:              // *  Further Details
 251:              // *  ===============
 252:              // *
 253:              // *  The matrices Q and P are represented as products of elementary
 254:              // *  reflectors:
 255:              // *
 256:              // *  If m >= n,
 257:              // *
 258:              // *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
 259:              // *
 260:              // *  Each H(i) and G(i) has the form:
 261:              // *
 262:              // *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
 263:              // *
 264:              // *  where tauq and taup are real scalars, and v and u are real vectors;
 265:              // *  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
 266:              // *  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
 267:              // *  tauq is stored in TAUQ(i) and taup in TAUP(i).
 268:              // *
 269:              // *  If m < n,
 270:              // *
 271:              // *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
 272:              // *
 273:              // *  Each H(i) and G(i) has the form:
 274:              // *
 275:              // *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
 276:              // *
 277:              // *  where tauq and taup are real scalars, and v and u are real vectors;
 278:              // *  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
 279:              // *  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
 280:              // *  tauq is stored in TAUQ(i) and taup in TAUP(i).
 281:              // *
 282:              // *  The contents of A on exit are illustrated by the following examples:
 283:              // *
 284:              // *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
 285:              // *
 286:              // *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
 287:              // *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
 288:              // *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
 289:              // *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
 290:              // *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
 291:              // *    (  v1  v2  v3  v4  v5 )
 292:              // *
 293:              // *  where d and e denote diagonal and off-diagonal elements of B, vi
 294:              // *  denotes an element of the vector defining H(i), and ui an element of
 295:              // *  the vector defining G(i).
 296:              // *
 297:              // *  =====================================================================
 298:              // *
 299:              // *     .. Parameters ..
 300:              // *     ..
 301:              // *     .. Local Scalars ..
 302:              // *     ..
 303:              // *     .. External Subroutines ..
 304:              // *     ..
 305:              // *     .. Intrinsic Functions ..
 306:              //      INTRINSIC          MAX, MIN;
 307:              // *     ..
 308:              // *     .. Executable Statements ..
 309:              // *
 310:              // *     Test the input parameters
 311:              // *
 312:   
 313:              #endregion
 314:   
 315:   
 316:              #region Body
 317:              
 318:              INFO = 0;
 319:              if (M < 0)
 320:              {
 321:                  INFO =  - 1;
 322:              }
 323:              else
 324:              {
 325:                  if (N < 0)
 326:                  {
 327:                      INFO =  - 2;
 328:                  }
 329:                  else
 330:                  {
 331:                      if (LDA < Math.Max(1, M))
 332:                      {
 333:                          INFO =  - 4;
 334:                      }
 335:                  }
 336:              }
 337:              if (INFO < 0)
 338:              {
 339:                  this._xerbla.Run("DGEBD2",  - INFO);
 340:                  return;
 341:              }
 342:              // *
 343:              if (M >= N)
 344:              {
 345:                  // *
 346:                  // *        Reduce to upper bidiagonal form
 347:                  // *
 348:                  for (I = 1; I <= N; I++)
 349:                  {
 350:                      // *
 351:                      // *           Generate elementary reflector H(i) to annihilate A(i+1:m,i)
 352:                      // *
 353:                      this._dlarfg.Run(M - I + 1, ref A[I+I * LDA + o_a], ref A, Math.Min(I + 1, M)+I * LDA + o_a, 1, ref TAUQ[I + o_tauq]);
 354:                      D[I + o_d] = A[I+I * LDA + o_a];
 355:                      A[I+I * LDA + o_a] = ONE;
 356:                      // *
 357:                      // *           Apply H(i) to A(i:m,i+1:n) from the left
 358:                      // *
 359:                      if (I < N)
 360:                      {
 361:                          this._dlarf.Run("Left", M - I + 1, N - I, A, I+I * LDA + o_a, 1, TAUQ[I + o_tauq]
 362:                                          , ref A, I+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
 363:                      }
 364:                      A[I+I * LDA + o_a] = D[I + o_d];
 365:                      // *
 366:                      if (I < N)
 367:                      {
 368:                          // *
 369:                          // *              Generate elementary reflector G(i) to annihilate
 370:                          // *              A(i,i+2:n)
 371:                          // *
 372:                          this._dlarfg.Run(N - I, ref A[I+(I + 1) * LDA + o_a], ref A, I+Math.Min(I + 2, N) * LDA + o_a, LDA, ref TAUP[I + o_taup]);
 373:                          E[I + o_e] = A[I+(I + 1) * LDA + o_a];
 374:                          A[I+(I + 1) * LDA + o_a] = ONE;
 375:                          // *
 376:                          // *              Apply G(i) to A(i+1:m,i+1:n) from the right
 377:                          // *
 378:                          this._dlarf.Run("Right", M - I, N - I, A, I+(I + 1) * LDA + o_a, LDA, TAUP[I + o_taup]
 379:                                          , ref A, I + 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
 380:                          A[I+(I + 1) * LDA + o_a] = E[I + o_e];
 381:                      }
 382:                      else
 383:                      {
 384:                          TAUP[I + o_taup] = ZERO;
 385:                      }
 386:                  }
 387:              }
 388:              else
 389:              {
 390:                  // *
 391:                  // *        Reduce to lower bidiagonal form
 392:                  // *
 393:                  for (I = 1; I <= M; I++)
 394:                  {
 395:                      // *
 396:                      // *           Generate elementary reflector G(i) to annihilate A(i,i+1:n)
 397:                      // *
 398:                      this._dlarfg.Run(N - I + 1, ref A[I+I * LDA + o_a], ref A, I+Math.Min(I + 1, N) * LDA + o_a, LDA, ref TAUP[I + o_taup]);
 399:                      D[I + o_d] = A[I+I * LDA + o_a];
 400:                      A[I+I * LDA + o_a] = ONE;
 401:                      // *
 402:                      // *           Apply G(i) to A(i+1:m,i:n) from the right
 403:                      // *
 404:                      if (I < M)
 405:                      {
 406:                          this._dlarf.Run("Right", M - I, N - I + 1, A, I+I * LDA + o_a, LDA, TAUP[I + o_taup]
 407:                                          , ref A, I + 1+I * LDA + o_a, LDA, ref WORK, offset_work);
 408:                      }
 409:                      A[I+I * LDA + o_a] = D[I + o_d];
 410:                      // *
 411:                      if (I < M)
 412:                      {
 413:                          // *
 414:                          // *              Generate elementary reflector H(i) to annihilate
 415:                          // *              A(i+2:m,i)
 416:                          // *
 417:                          this._dlarfg.Run(M - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, M)+I * LDA + o_a, 1, ref TAUQ[I + o_tauq]);
 418:                          E[I + o_e] = A[I + 1+I * LDA + o_a];
 419:                          A[I + 1+I * LDA + o_a] = ONE;
 420:                          // *
 421:                          // *              Apply H(i) to A(i+1:m,i+1:n) from the left
 422:                          // *
 423:                          this._dlarf.Run("Left", M - I, N - I, A, I + 1+I * LDA + o_a, 1, TAUQ[I + o_tauq]
 424:                                          , ref A, I + 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
 425:                          A[I + 1+I * LDA + o_a] = E[I + o_e];
 426:                      }
 427:                      else
 428:                      {
 429:                          TAUQ[I + o_tauq] = ZERO;
 430:                      }
 431:                  }
 432:              }
 433:              return;
 434:              // *
 435:              // *     End of DGEBD2
 436:              // *
 437:   
 438:              #endregion
 439:   
 440:          }
 441:      }
 442:  }