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:      /// DGEBRD reduces a general real M-by-N matrix A to upper or lower
  27:      /// bidiagonal form B by an orthogonal transformation: Q**T * 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 DGEBRD
  33:      {
  34:      
  35:   
  36:          #region Dependencies
  37:          
  38:          DGEBD2 _dgebd2; DGEMM _dgemm; DLABRD _dlabrd; XERBLA _xerbla; ILAENV _ilaenv; 
  39:   
  40:          #endregion
  41:   
  42:   
  43:          #region Fields
  44:          
  45:          const double ONE = 1.0E+0; bool LQUERY = false; int I = 0; int IINFO = 0; int J = 0; int LDWRKX = 0; int LDWRKY = 0; 
  46:          int LWKOPT = 0;int MINMN = 0; int NB = 0; int NBMIN = 0; int NX = 0; double WS = 0; 
  47:   
  48:          #endregion
  49:   
  50:          public DGEBRD(DGEBD2 dgebd2, DGEMM dgemm, DLABRD dlabrd, XERBLA xerbla, ILAENV ilaenv)
  51:          {
  52:      
  53:   
  54:              #region Set Dependencies
  55:              
  56:              this._dgebd2 = dgebd2; this._dgemm = dgemm; this._dlabrd = dlabrd; this._xerbla = xerbla; this._ilaenv = ilaenv; 
  57:   
  58:              #endregion
  59:   
  60:          }
  61:      
  62:          public DGEBRD()
  63:          {
  64:      
  65:   
  66:              #region Dependencies (Initialization)
  67:              
  68:              LSAME lsame = new LSAME();
  69:              XERBLA xerbla = new XERBLA();
  70:              DLAMC3 dlamc3 = new DLAMC3();
  71:              DLAPY2 dlapy2 = new DLAPY2();
  72:              DNRM2 dnrm2 = new DNRM2();
  73:              DSCAL dscal = new DSCAL();
  74:              IEEECK ieeeck = new IEEECK();
  75:              IPARMQ iparmq = new IPARMQ();
  76:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  77:              DGER dger = new DGER(xerbla);
  78:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  79:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  80:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  81:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  82:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  83:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  84:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  85:              DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);
  86:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  87:              DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);
  88:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  89:   
  90:              #endregion
  91:   
  92:   
  93:              #region Set Dependencies
  94:              
  95:              this._dgebd2 = dgebd2; this._dgemm = dgemm; this._dlabrd = dlabrd; this._xerbla = xerbla; this._ilaenv = ilaenv; 
  96:   
  97:              #endregion
  98:   
  99:          }
 100:          /// <summary>
 101:          /// Purpose
 102:          /// =======
 103:          /// 
 104:          /// DGEBRD reduces a general real M-by-N matrix A to upper or lower
 105:          /// bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
 106:          /// 
 107:          /// If m .GE. n, B is upper bidiagonal; if m .LT. n, B is lower bidiagonal.
 108:          /// 
 109:          ///</summary>
 110:          /// <param name="M">
 111:          /// (input) INTEGER
 112:          /// The number of rows in the matrix A.  M .GE. 0.
 113:          ///</param>
 114:          /// <param name="N">
 115:          /// (input) INTEGER
 116:          /// The number of columns in the matrix A.  N .GE. 0.
 117:          ///</param>
 118:          /// <param name="A">
 119:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 120:          /// On entry, the M-by-N general matrix to be reduced.
 121:          /// On exit,
 122:          /// if m .GE. n, the diagonal and the first superdiagonal are
 123:          /// overwritten with the upper bidiagonal matrix B; the
 124:          /// elements below the diagonal, with the array TAUQ, represent
 125:          /// the orthogonal matrix Q as a product of elementary
 126:          /// reflectors, and the elements above the first superdiagonal,
 127:          /// with the array TAUP, represent the orthogonal matrix P as
 128:          /// a product of elementary reflectors;
 129:          /// if m .LT. n, the diagonal and the first subdiagonal are
 130:          /// overwritten with the lower bidiagonal matrix B; the
 131:          /// elements below the first subdiagonal, with the array TAUQ,
 132:          /// represent the orthogonal matrix Q as a product of
 133:          /// elementary reflectors, and the elements above the diagonal,
 134:          /// with the array TAUP, represent the orthogonal matrix P as
 135:          /// a product of elementary reflectors.
 136:          /// See Further Details.
 137:          ///</param>
 138:          /// <param name="LDA">
 139:          /// (input) INTEGER
 140:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 141:          ///</param>
 142:          /// <param name="D">
 143:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 144:          /// The diagonal elements of the bidiagonal matrix B:
 145:          /// D(i) = A(i,i).
 146:          ///</param>
 147:          /// <param name="E">
 148:          /// (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
 149:          /// The off-diagonal elements of the bidiagonal matrix B:
 150:          /// if m .GE. n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
 151:          /// if m .LT. n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
 152:          ///</param>
 153:          /// <param name="TAUQ">
 154:          /// (output) DOUBLE PRECISION array dimension (min(M,N))
 155:          /// The scalar factors of the elementary reflectors which
 156:          /// represent the orthogonal matrix Q. See Further Details.
 157:          ///</param>
 158:          /// <param name="TAUP">
 159:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 160:          /// The scalar factors of the elementary reflectors which
 161:          /// represent the orthogonal matrix P. See Further Details.
 162:          ///</param>
 163:          /// <param name="WORK">
 164:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 165:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 166:          ///</param>
 167:          /// <param name="LWORK">
 168:          /// (input) INTEGER
 169:          /// The length of the array WORK.  LWORK .GE. max(1,M,N).
 170:          /// For optimum performance LWORK .GE. (M+N)*NB, where NB
 171:          /// is the optimal blocksize.
 172:          /// 
 173:          /// If LWORK = -1, then a workspace query is assumed; the routine
 174:          /// only calculates the optimal size of the WORK array, returns
 175:          /// this value as the first entry of the WORK array, and no error
 176:          /// message related to LWORK is issued by XERBLA.
 177:          ///</param>
 178:          /// <param name="INFO">
 179:          /// (output) INTEGER
 180:          /// = 0:  successful exit
 181:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 182:          ///</param>
 183:          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
 184:                           , ref double[] TAUQ, int offset_tauq, ref double[] TAUP, int offset_taup, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
 185:          {
 186:   
 187:              #region Array Index Correction
 188:              
 189:               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; 
 190:               int o_taup = -1 + offset_taup; int o_work = -1 + offset_work; 
 191:   
 192:              #endregion
 193:   
 194:   
 195:              #region Prolog
 196:              
 197:              // *
 198:              // *  -- LAPACK routine (version 3.1) --
 199:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 200:              // *     November 2006
 201:              // *
 202:              // *     .. Scalar Arguments ..
 203:              // *     ..
 204:              // *     .. Array Arguments ..
 205:              // *     ..
 206:              // *
 207:              // *  Purpose
 208:              // *  =======
 209:              // *
 210:              // *  DGEBRD reduces a general real M-by-N matrix A to upper or lower
 211:              // *  bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
 212:              // *
 213:              // *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
 214:              // *
 215:              // *  Arguments
 216:              // *  =========
 217:              // *
 218:              // *  M       (input) INTEGER
 219:              // *          The number of rows in the matrix A.  M >= 0.
 220:              // *
 221:              // *  N       (input) INTEGER
 222:              // *          The number of columns in the matrix A.  N >= 0.
 223:              // *
 224:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 225:              // *          On entry, the M-by-N general matrix to be reduced.
 226:              // *          On exit,
 227:              // *          if m >= n, the diagonal and the first superdiagonal are
 228:              // *            overwritten with the upper bidiagonal matrix B; the
 229:              // *            elements below the diagonal, with the array TAUQ, represent
 230:              // *            the orthogonal matrix Q as a product of elementary
 231:              // *            reflectors, and the elements above the first superdiagonal,
 232:              // *            with the array TAUP, represent the orthogonal matrix P as
 233:              // *            a product of elementary reflectors;
 234:              // *          if m < n, the diagonal and the first subdiagonal are
 235:              // *            overwritten with the lower bidiagonal matrix B; the
 236:              // *            elements below the first subdiagonal, with the array TAUQ,
 237:              // *            represent the orthogonal matrix Q as a product of
 238:              // *            elementary reflectors, and the elements above the diagonal,
 239:              // *            with the array TAUP, represent the orthogonal matrix P as
 240:              // *            a product of elementary reflectors.
 241:              // *          See Further Details.
 242:              // *
 243:              // *  LDA     (input) INTEGER
 244:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 245:              // *
 246:              // *  D       (output) DOUBLE PRECISION array, dimension (min(M,N))
 247:              // *          The diagonal elements of the bidiagonal matrix B:
 248:              // *          D(i) = A(i,i).
 249:              // *
 250:              // *  E       (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
 251:              // *          The off-diagonal elements of the bidiagonal matrix B:
 252:              // *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
 253:              // *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
 254:              // *
 255:              // *  TAUQ    (output) DOUBLE PRECISION array dimension (min(M,N))
 256:              // *          The scalar factors of the elementary reflectors which
 257:              // *          represent the orthogonal matrix Q. See Further Details.
 258:              // *
 259:              // *  TAUP    (output) DOUBLE PRECISION array, dimension (min(M,N))
 260:              // *          The scalar factors of the elementary reflectors which
 261:              // *          represent the orthogonal matrix P. See Further Details.
 262:              // *
 263:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 264:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 265:              // *
 266:              // *  LWORK   (input) INTEGER
 267:              // *          The length of the array WORK.  LWORK >= max(1,M,N).
 268:              // *          For optimum performance LWORK >= (M+N)*NB, where NB
 269:              // *          is the optimal blocksize.
 270:              // *
 271:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 272:              // *          only calculates the optimal size of the WORK array, returns
 273:              // *          this value as the first entry of the WORK array, and no error
 274:              // *          message related to LWORK is issued by XERBLA.
 275:              // *
 276:              // *  INFO    (output) INTEGER
 277:              // *          = 0:  successful exit
 278:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 279:              // *
 280:              // *  Further Details
 281:              // *  ===============
 282:              // *
 283:              // *  The matrices Q and P are represented as products of elementary
 284:              // *  reflectors:
 285:              // *
 286:              // *  If m >= n,
 287:              // *
 288:              // *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
 289:              // *
 290:              // *  Each H(i) and G(i) has the form:
 291:              // *
 292:              // *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
 293:              // *
 294:              // *  where tauq and taup are real scalars, and v and u are real vectors;
 295:              // *  v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
 296:              // *  u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
 297:              // *  tauq is stored in TAUQ(i) and taup in TAUP(i).
 298:              // *
 299:              // *  If m < n,
 300:              // *
 301:              // *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
 302:              // *
 303:              // *  Each H(i) and G(i) has the form:
 304:              // *
 305:              // *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
 306:              // *
 307:              // *  where tauq and taup are real scalars, and v and u are real vectors;
 308:              // *  v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
 309:              // *  u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
 310:              // *  tauq is stored in TAUQ(i) and taup in TAUP(i).
 311:              // *
 312:              // *  The contents of A on exit are illustrated by the following examples:
 313:              // *
 314:              // *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
 315:              // *
 316:              // *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
 317:              // *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
 318:              // *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
 319:              // *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
 320:              // *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
 321:              // *    (  v1  v2  v3  v4  v5 )
 322:              // *
 323:              // *  where d and e denote diagonal and off-diagonal elements of B, vi
 324:              // *  denotes an element of the vector defining H(i), and ui an element of
 325:              // *  the vector defining G(i).
 326:              // *
 327:              // *  =====================================================================
 328:              // *
 329:              // *     .. Parameters ..
 330:              // *     ..
 331:              // *     .. Local Scalars ..
 332:              // *     ..
 333:              // *     .. External Subroutines ..
 334:              // *     ..
 335:              // *     .. Intrinsic Functions ..
 336:              //      INTRINSIC          DBLE, MAX, MIN;
 337:              // *     ..
 338:              // *     .. External Functions ..
 339:              // *     ..
 340:              // *     .. Executable Statements ..
 341:              // *
 342:              // *     Test the input parameters
 343:              // *
 344:   
 345:              #endregion
 346:   
 347:   
 348:              #region Body
 349:              
 350:              INFO = 0;
 351:              NB = Math.Max(1, this._ilaenv.Run(1, "DGEBRD", " ", M, N,  - 1,  - 1));
 352:              LWKOPT = (M + N) * NB;
 353:              WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
 354:              LQUERY = (LWORK ==  - 1);
 355:              if (M < 0)
 356:              {
 357:                  INFO =  - 1;
 358:              }
 359:              else
 360:              {
 361:                  if (N < 0)
 362:                  {
 363:                      INFO =  - 2;
 364:                  }
 365:                  else
 366:                  {
 367:                      if (LDA < Math.Max(1, M))
 368:                      {
 369:                          INFO =  - 4;
 370:                      }
 371:                      else
 372:                      {
 373:                          if (LWORK < Math.Max(1, Math.Max(M, N)) && !LQUERY)
 374:                          {
 375:                              INFO =  - 10;
 376:                          }
 377:                      }
 378:                  }
 379:              }
 380:              if (INFO < 0)
 381:              {
 382:                  this._xerbla.Run("DGEBRD",  - INFO);
 383:                  return;
 384:              }
 385:              else
 386:              {
 387:                  if (LQUERY)
 388:                  {
 389:                      return;
 390:                  }
 391:              }
 392:              // *
 393:              // *     Quick return if possible
 394:              // *
 395:              MINMN = Math.Min(M, N);
 396:              if (MINMN == 0)
 397:              {
 398:                  WORK[1 + o_work] = 1;
 399:                  return;
 400:              }
 401:              // *
 402:              WS = Math.Max(M, N);
 403:              LDWRKX = M;
 404:              LDWRKY = N;
 405:              // *
 406:              if (NB > 1 && NB < MINMN)
 407:              {
 408:                  // *
 409:                  // *        Set the crossover point NX.
 410:                  // *
 411:                  NX = Math.Max(NB, this._ilaenv.Run(3, "DGEBRD", " ", M, N,  - 1,  - 1));
 412:                  // *
 413:                  // *        Determine when to switch from blocked to unblocked code.
 414:                  // *
 415:                  if (NX < MINMN)
 416:                  {
 417:                      WS = (M + N) * NB;
 418:                      if (LWORK < WS)
 419:                      {
 420:                          // *
 421:                          // *              Not enough work space for the optimal NB, consider using
 422:                          // *              a smaller block size.
 423:                          // *
 424:                          NBMIN = this._ilaenv.Run(2, "DGEBRD", " ", M, N,  - 1,  - 1);
 425:                          if (LWORK >= (M + N) * NBMIN)
 426:                          {
 427:                              NB = LWORK / (M + N);
 428:                          }
 429:                          else
 430:                          {
 431:                              NB = 1;
 432:                              NX = MINMN;
 433:                          }
 434:                      }
 435:                  }
 436:              }
 437:              else
 438:              {
 439:                  NX = MINMN;
 440:              }
 441:              // *
 442:              for (I = 1; (NB >= 0) ? (I <= MINMN - NX) : (I >= MINMN - NX); I += NB)
 443:              {
 444:                  // *
 445:                  // *        Reduce rows and columns i:i+nb-1 to bidiagonal form and return
 446:                  // *        the matrices X and Y which are needed to update the unreduced
 447:                  // *        part of the matrix
 448:                  // *
 449:                  this._dlabrd.Run(M - I + 1, N - I + 1, NB, ref A, I+I * LDA + o_a, LDA, ref D, I + o_d
 450:                                   , ref E, I + o_e, ref TAUQ, I + o_tauq, ref TAUP, I + o_taup, ref WORK, offset_work, LDWRKX, ref WORK, LDWRKX * NB + 1 + o_work
 451:                                   , LDWRKY);
 452:                  // *
 453:                  // *        Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
 454:                  // *        of the form  A := A - V*Y' - X*U'
 455:                  // *
 456:                  this._dgemm.Run("No transpose", "Transpose", M - I - NB + 1, N - I - NB + 1, NB,  - ONE
 457:                                  , A, I + NB+I * LDA + o_a, LDA, WORK, LDWRKX * NB + NB + 1 + o_work, LDWRKY, ONE, ref A, I + NB+(I + NB) * LDA + o_a
 458:                                  , LDA);
 459:                  this._dgemm.Run("No transpose", "No transpose", M - I - NB + 1, N - I - NB + 1, NB,  - ONE
 460:                                  , WORK, NB + 1 + o_work, LDWRKX, A, I+(I + NB) * LDA + o_a, LDA, ONE, ref A, I + NB+(I + NB) * LDA + o_a
 461:                                  , LDA);
 462:                  // *
 463:                  // *        Copy diagonal and off-diagonal elements of B back into A
 464:                  // *
 465:                  if (M >= N)
 466:                  {
 467:                      for (J = I; J <= I + NB - 1; J++)
 468:                      {
 469:                          A[J+J * LDA + o_a] = D[J + o_d];
 470:                          A[J+(J + 1) * LDA + o_a] = E[J + o_e];
 471:                      }
 472:                  }
 473:                  else
 474:                  {
 475:                      for (J = I; J <= I + NB - 1; J++)
 476:                      {
 477:                          A[J+J * LDA + o_a] = D[J + o_d];
 478:                          A[J + 1+J * LDA + o_a] = E[J + o_e];
 479:                      }
 480:                  }
 481:              }
 482:              // *
 483:              // *     Use unblocked code to reduce the remainder of the matrix
 484:              // *
 485:              this._dgebd2.Run(M - I + 1, N - I + 1, ref A, I+I * LDA + o_a, LDA, ref D, I + o_d, ref E, I + o_e
 486:                               , ref TAUQ, I + o_tauq, ref TAUP, I + o_taup, ref WORK, offset_work, ref IINFO);
 487:              WORK[1 + o_work] = WS;
 488:              return;
 489:              // *
 490:              // *     End of DGEBRD
 491:              // *
 492:   
 493:              #endregion
 494:   
 495:          }
 496:      }
 497:  }