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:      /// DLATRD reduces NB rows and columns of a real symmetric matrix A to
  27:      /// symmetric tridiagonal form by an orthogonal similarity
  28:      /// transformation Q' * A * Q, and returns the matrices V and W which are
  29:      /// needed to apply the transformation to the unreduced part of A.
  30:      /// 
  31:      /// If UPLO = 'U', DLATRD reduces the last NB rows and columns of a
  32:      /// matrix, of which the upper triangle is supplied;
  33:      /// if UPLO = 'L', DLATRD reduces the first NB rows and columns of a
  34:      /// matrix, of which the lower triangle is supplied.
  35:      /// 
  36:      /// This is an auxiliary routine called by DSYTRD.
  37:      /// 
  38:      ///</summary>
  39:      public class DLATRD
  40:      {
  41:      
  42:   
  43:          #region Dependencies
  44:          
  45:          DAXPY _daxpy; DGEMV _dgemv; DLARFG _dlarfg; DSCAL _dscal; DSYMV _dsymv; LSAME _lsame; DDOT _ddot; 
  46:   
  47:          #endregion
  48:   
  49:   
  50:          #region Fields
  51:          
  52:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double HALF = 0.5E+0; int I = 0; int IW = 0; 
  53:          double ALPHA = 0;
  54:   
  55:          #endregion
  56:   
  57:          public DLATRD(DAXPY daxpy, DGEMV dgemv, DLARFG dlarfg, DSCAL dscal, DSYMV dsymv, LSAME lsame, DDOT ddot)
  58:          {
  59:      
  60:   
  61:              #region Set Dependencies
  62:              
  63:              this._daxpy = daxpy; this._dgemv = dgemv; this._dlarfg = dlarfg; this._dscal = dscal; this._dsymv = dsymv; 
  64:              this._lsame = lsame;this._ddot = ddot; 
  65:   
  66:              #endregion
  67:   
  68:          }
  69:      
  70:          public DLATRD()
  71:          {
  72:      
  73:   
  74:              #region Dependencies (Initialization)
  75:              
  76:              DAXPY daxpy = new DAXPY();
  77:              LSAME lsame = new LSAME();
  78:              XERBLA xerbla = new XERBLA();
  79:              DLAMC3 dlamc3 = new DLAMC3();
  80:              DLAPY2 dlapy2 = new DLAPY2();
  81:              DNRM2 dnrm2 = new DNRM2();
  82:              DSCAL dscal = new DSCAL();
  83:              DDOT ddot = new DDOT();
  84:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  85:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  86:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  87:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  88:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  89:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  90:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  91:              DSYMV dsymv = new DSYMV(lsame, xerbla);
  92:   
  93:              #endregion
  94:   
  95:   
  96:              #region Set Dependencies
  97:              
  98:              this._daxpy = daxpy; this._dgemv = dgemv; this._dlarfg = dlarfg; this._dscal = dscal; this._dsymv = dsymv; 
  99:              this._lsame = lsame;this._ddot = ddot; 
 100:   
 101:              #endregion
 102:   
 103:          }
 104:          /// <summary>
 105:          /// Purpose
 106:          /// =======
 107:          /// 
 108:          /// DLATRD reduces NB rows and columns of a real symmetric matrix A to
 109:          /// symmetric tridiagonal form by an orthogonal similarity
 110:          /// transformation Q' * A * Q, and returns the matrices V and W which are
 111:          /// needed to apply the transformation to the unreduced part of A.
 112:          /// 
 113:          /// If UPLO = 'U', DLATRD reduces the last NB rows and columns of a
 114:          /// matrix, of which the upper triangle is supplied;
 115:          /// if UPLO = 'L', DLATRD reduces the first NB rows and columns of a
 116:          /// matrix, of which the lower triangle is supplied.
 117:          /// 
 118:          /// This is an auxiliary routine called by DSYTRD.
 119:          /// 
 120:          ///</summary>
 121:          /// <param name="UPLO">
 122:          /// (input) CHARACTER*1
 123:          /// Specifies whether the upper or lower triangular part of the
 124:          /// symmetric matrix A is stored:
 125:          /// = 'U': Upper triangular
 126:          /// = 'L': Lower triangular
 127:          ///</param>
 128:          /// <param name="N">
 129:          /// (input) INTEGER
 130:          /// The order of the matrix A.
 131:          ///</param>
 132:          /// <param name="NB">
 133:          /// (input) INTEGER
 134:          /// The number of rows and columns to be reduced.
 135:          ///</param>
 136:          /// <param name="A">
 137:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 138:          /// On entry, the symmetric matrix A.  If UPLO = 'U', the leading
 139:          /// n-by-n upper triangular part of A contains the upper
 140:          /// triangular part of the matrix A, and the strictly lower
 141:          /// triangular part of A is not referenced.  If UPLO = 'L', the
 142:          /// leading n-by-n lower triangular part of A contains the lower
 143:          /// triangular part of the matrix A, and the strictly upper
 144:          /// triangular part of A is not referenced.
 145:          /// On exit:
 146:          /// if UPLO = 'U', the last NB columns have been reduced to
 147:          /// tridiagonal form, with the diagonal elements overwriting
 148:          /// the diagonal elements of A; the elements above the diagonal
 149:          /// with the array TAU, represent the orthogonal matrix Q as a
 150:          /// product of elementary reflectors;
 151:          /// if UPLO = 'L', the first NB columns have been reduced to
 152:          /// tridiagonal form, with the diagonal elements overwriting
 153:          /// the diagonal elements of A; the elements below the diagonal
 154:          /// with the array TAU, represent the  orthogonal matrix Q as a
 155:          /// product of elementary reflectors.
 156:          /// See Further Details.
 157:          ///</param>
 158:          /// <param name="LDA">
 159:          /// (input) INTEGER
 160:          /// The leading dimension of the array A.  LDA .GE. (1,N).
 161:          ///</param>
 162:          /// <param name="E">
 163:          /// (output) DOUBLE PRECISION array, dimension (N-1)
 164:          /// If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
 165:          /// elements of the last NB columns of the reduced matrix;
 166:          /// if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
 167:          /// the first NB columns of the reduced matrix.
 168:          ///</param>
 169:          /// <param name="TAU">
 170:          /// (output) DOUBLE PRECISION array, dimension (N-1)
 171:          /// The scalar factors of the elementary reflectors, stored in
 172:          /// TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
 173:          /// See Further Details.
 174:          ///</param>
 175:          /// <param name="W">
 176:          /// (output) DOUBLE PRECISION array, dimension (LDW,NB)
 177:          /// The n-by-nb matrix W required to update the unreduced part
 178:          /// of A.
 179:          ///</param>
 180:          /// <param name="LDW">
 181:          /// (input) INTEGER
 182:          /// The leading dimension of the array W. LDW .GE. max(1,N).
 183:          ///</param>
 184:          public void Run(string UPLO, int N, int NB, ref double[] A, int offset_a, int LDA, ref double[] E, int offset_e
 185:                           , ref double[] TAU, int offset_tau, ref double[] W, int offset_w, int LDW)
 186:          {
 187:   
 188:              #region Array Index Correction
 189:              
 190:               int o_a = -1 - LDA + offset_a;  int o_e = -1 + offset_e;  int o_tau = -1 + offset_tau; 
 191:               int o_w = -1 - LDW + offset_w;
 192:   
 193:              #endregion
 194:   
 195:   
 196:              #region Strings
 197:              
 198:              UPLO = UPLO.Substring(0, 1);  
 199:   
 200:              #endregion
 201:   
 202:   
 203:              #region Prolog
 204:              
 205:              // *
 206:              // *  -- LAPACK auxiliary routine (version 3.1) --
 207:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 208:              // *     November 2006
 209:              // *
 210:              // *     .. Scalar Arguments ..
 211:              // *     ..
 212:              // *     .. Array Arguments ..
 213:              // *     ..
 214:              // *
 215:              // *  Purpose
 216:              // *  =======
 217:              // *
 218:              // *  DLATRD reduces NB rows and columns of a real symmetric matrix A to
 219:              // *  symmetric tridiagonal form by an orthogonal similarity
 220:              // *  transformation Q' * A * Q, and returns the matrices V and W which are
 221:              // *  needed to apply the transformation to the unreduced part of A.
 222:              // *
 223:              // *  If UPLO = 'U', DLATRD reduces the last NB rows and columns of a
 224:              // *  matrix, of which the upper triangle is supplied;
 225:              // *  if UPLO = 'L', DLATRD reduces the first NB rows and columns of a
 226:              // *  matrix, of which the lower triangle is supplied.
 227:              // *
 228:              // *  This is an auxiliary routine called by DSYTRD.
 229:              // *
 230:              // *  Arguments
 231:              // *  =========
 232:              // *
 233:              // *  UPLO    (input) CHARACTER*1
 234:              // *          Specifies whether the upper or lower triangular part of the
 235:              // *          symmetric matrix A is stored:
 236:              // *          = 'U': Upper triangular
 237:              // *          = 'L': Lower triangular
 238:              // *
 239:              // *  N       (input) INTEGER
 240:              // *          The order of the matrix A.
 241:              // *
 242:              // *  NB      (input) INTEGER
 243:              // *          The number of rows and columns to be reduced.
 244:              // *
 245:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 246:              // *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
 247:              // *          n-by-n upper triangular part of A contains the upper
 248:              // *          triangular part of the matrix A, and the strictly lower
 249:              // *          triangular part of A is not referenced.  If UPLO = 'L', the
 250:              // *          leading n-by-n lower triangular part of A contains the lower
 251:              // *          triangular part of the matrix A, and the strictly upper
 252:              // *          triangular part of A is not referenced.
 253:              // *          On exit:
 254:              // *          if UPLO = 'U', the last NB columns have been reduced to
 255:              // *            tridiagonal form, with the diagonal elements overwriting
 256:              // *            the diagonal elements of A; the elements above the diagonal
 257:              // *            with the array TAU, represent the orthogonal matrix Q as a
 258:              // *            product of elementary reflectors;
 259:              // *          if UPLO = 'L', the first NB columns have been reduced to
 260:              // *            tridiagonal form, with the diagonal elements overwriting
 261:              // *            the diagonal elements of A; the elements below the diagonal
 262:              // *            with the array TAU, represent the  orthogonal matrix Q as a
 263:              // *            product of elementary reflectors.
 264:              // *          See Further Details.
 265:              // *
 266:              // *  LDA     (input) INTEGER
 267:              // *          The leading dimension of the array A.  LDA >= (1,N).
 268:              // *
 269:              // *  E       (output) DOUBLE PRECISION array, dimension (N-1)
 270:              // *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
 271:              // *          elements of the last NB columns of the reduced matrix;
 272:              // *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
 273:              // *          the first NB columns of the reduced matrix.
 274:              // *
 275:              // *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
 276:              // *          The scalar factors of the elementary reflectors, stored in
 277:              // *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
 278:              // *          See Further Details.
 279:              // *
 280:              // *  W       (output) DOUBLE PRECISION array, dimension (LDW,NB)
 281:              // *          The n-by-nb matrix W required to update the unreduced part
 282:              // *          of A.
 283:              // *
 284:              // *  LDW     (input) INTEGER
 285:              // *          The leading dimension of the array W. LDW >= max(1,N).
 286:              // *
 287:              // *  Further Details
 288:              // *  ===============
 289:              // *
 290:              // *  If UPLO = 'U', the matrix Q is represented as a product of elementary
 291:              // *  reflectors
 292:              // *
 293:              // *     Q = H(n) H(n-1) . . . H(n-nb+1).
 294:              // *
 295:              // *  Each H(i) has the form
 296:              // *
 297:              // *     H(i) = I - tau * v * v'
 298:              // *
 299:              // *  where tau is a real scalar, and v is a real vector with
 300:              // *  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
 301:              // *  and tau in TAU(i-1).
 302:              // *
 303:              // *  If UPLO = 'L', the matrix Q is represented as a product of elementary
 304:              // *  reflectors
 305:              // *
 306:              // *     Q = H(1) H(2) . . . H(nb).
 307:              // *
 308:              // *  Each H(i) has the form
 309:              // *
 310:              // *     H(i) = I - tau * v * v'
 311:              // *
 312:              // *  where tau is a real scalar, and v is a real vector with
 313:              // *  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
 314:              // *  and tau in TAU(i).
 315:              // *
 316:              // *  The elements of the vectors v together form the n-by-nb matrix V
 317:              // *  which is needed, with W, to apply the transformation to the unreduced
 318:              // *  part of the matrix, using a symmetric rank-2k update of the form:
 319:              // *  A := A - V*W' - W*V'.
 320:              // *
 321:              // *  The contents of A on exit are illustrated by the following examples
 322:              // *  with n = 5 and nb = 2:
 323:              // *
 324:              // *  if UPLO = 'U':                       if UPLO = 'L':
 325:              // *
 326:              // *    (  a   a   a   v4  v5 )              (  d                  )
 327:              // *    (      a   a   v4  v5 )              (  1   d              )
 328:              // *    (          a   1   v5 )              (  v1  1   a          )
 329:              // *    (              d   1  )              (  v1  v2  a   a      )
 330:              // *    (                  d  )              (  v1  v2  a   a   a  )
 331:              // *
 332:              // *  where d denotes a diagonal element of the reduced matrix, a denotes
 333:              // *  an element of the original matrix that is unchanged, and vi denotes
 334:              // *  an element of the vector defining H(i).
 335:              // *
 336:              // *  =====================================================================
 337:              // *
 338:              // *     .. Parameters ..
 339:              // *     ..
 340:              // *     .. Local Scalars ..
 341:              // *     ..
 342:              // *     .. External Subroutines ..
 343:              // *     ..
 344:              // *     .. External Functions ..
 345:              // *     ..
 346:              // *     .. Intrinsic Functions ..
 347:              //      INTRINSIC          MIN;
 348:              // *     ..
 349:              // *     .. Executable Statements ..
 350:              // *
 351:              // *     Quick return if possible
 352:              // *
 353:   
 354:              #endregion
 355:   
 356:   
 357:              #region Body
 358:              
 359:              if (N <= 0) return;
 360:              // *
 361:              if (this._lsame.Run(UPLO, "U"))
 362:              {
 363:                  // *
 364:                  // *        Reduce last NB columns of upper triangle
 365:                  // *
 366:                  for (I = N; I >= N - NB + 1; I +=  - 1)
 367:                  {
 368:                      IW = I - N + NB;
 369:                      if (I < N)
 370:                      {
 371:                          // *
 372:                          // *              Update A(1:i,i)
 373:                          // *
 374:                          this._dgemv.Run("No transpose", I, N - I,  - ONE, A, 1+(I + 1) * LDA + o_a, LDA
 375:                                          , W, I+(IW + 1) * LDW + o_w, LDW, ONE, ref A, 1+I * LDA + o_a, 1);
 376:                          this._dgemv.Run("No transpose", I, N - I,  - ONE, W, 1+(IW + 1) * LDW + o_w, LDW
 377:                                          , A, I+(I + 1) * LDA + o_a, LDA, ONE, ref A, 1+I * LDA + o_a, 1);
 378:                      }
 379:                      if (I > 1)
 380:                      {
 381:                          // *
 382:                          // *              Generate elementary reflector H(i) to annihilate
 383:                          // *              A(1:i-2,i)
 384:                          // *
 385:                          this._dlarfg.Run(I - 1, ref A[I - 1+I * LDA + o_a], ref A, 1+I * LDA + o_a, 1, ref TAU[I - 1 + o_tau]);
 386:                          E[I - 1 + o_e] = A[I - 1+I * LDA + o_a];
 387:                          A[I - 1+I * LDA + o_a] = ONE;
 388:                          // *
 389:                          // *              Compute W(1:i-1,i)
 390:                          // *
 391:                          this._dsymv.Run("Upper", I - 1, ONE, A, offset_a, LDA, A, 1+I * LDA + o_a
 392:                                          , 1, ZERO, ref W, 1+IW * LDW + o_w, 1);
 393:                          if (I < N)
 394:                          {
 395:                              this._dgemv.Run("Transpose", I - 1, N - I, ONE, W, 1+(IW + 1) * LDW + o_w, LDW
 396:                                              , A, 1+I * LDA + o_a, 1, ZERO, ref W, I + 1+IW * LDW + o_w, 1);
 397:                              this._dgemv.Run("No transpose", I - 1, N - I,  - ONE, A, 1+(I + 1) * LDA + o_a, LDA
 398:                                              , W, I + 1+IW * LDW + o_w, 1, ONE, ref W, 1+IW * LDW + o_w, 1);
 399:                              this._dgemv.Run("Transpose", I - 1, N - I, ONE, A, 1+(I + 1) * LDA + o_a, LDA
 400:                                              , A, 1+I * LDA + o_a, 1, ZERO, ref W, I + 1+IW * LDW + o_w, 1);
 401:                              this._dgemv.Run("No transpose", I - 1, N - I,  - ONE, W, 1+(IW + 1) * LDW + o_w, LDW
 402:                                              , W, I + 1+IW * LDW + o_w, 1, ONE, ref W, 1+IW * LDW + o_w, 1);
 403:                          }
 404:                          this._dscal.Run(I - 1, TAU[I - 1 + o_tau], ref W, 1+IW * LDW + o_w, 1);
 405:                          ALPHA =  - HALF * TAU[I - 1 + o_tau] * this._ddot.Run(I - 1, W, 1+IW * LDW + o_w, 1, A, 1+I * LDA + o_a, 1);
 406:                          this._daxpy.Run(I - 1, ALPHA, A, 1+I * LDA + o_a, 1, ref W, 1+IW * LDW + o_w, 1);
 407:                      }
 408:                      // *
 409:                  }
 410:              }
 411:              else
 412:              {
 413:                  // *
 414:                  // *        Reduce first NB columns of lower triangle
 415:                  // *
 416:                  for (I = 1; I <= NB; I++)
 417:                  {
 418:                      // *
 419:                      // *           Update A(i:n,i)
 420:                      // *
 421:                      this._dgemv.Run("No transpose", N - I + 1, I - 1,  - ONE, A, I+1 * LDA + o_a, LDA
 422:                                      , W, I+1 * LDW + o_w, LDW, ONE, ref A, I+I * LDA + o_a, 1);
 423:                      this._dgemv.Run("No transpose", N - I + 1, I - 1,  - ONE, W, I+1 * LDW + o_w, LDW
 424:                                      , A, I+1 * LDA + o_a, LDA, ONE, ref A, I+I * LDA + o_a, 1);
 425:                      if (I < N)
 426:                      {
 427:                          // *
 428:                          // *              Generate elementary reflector H(i) to annihilate
 429:                          // *              A(i+2:n,i)
 430:                          // *
 431:                          this._dlarfg.Run(N - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, N)+I * LDA + o_a, 1, ref TAU[I + o_tau]);
 432:                          E[I + o_e] = A[I + 1+I * LDA + o_a];
 433:                          A[I + 1+I * LDA + o_a] = ONE;
 434:                          // *
 435:                          // *              Compute W(i+1:n,i)
 436:                          // *
 437:                          this._dsymv.Run("Lower", N - I, ONE, A, I + 1+(I + 1) * LDA + o_a, LDA, A, I + 1+I * LDA + o_a
 438:                                          , 1, ZERO, ref W, I + 1+I * LDW + o_w, 1);
 439:                          this._dgemv.Run("Transpose", N - I, I - 1, ONE, W, I + 1+1 * LDW + o_w, LDW
 440:                                          , A, I + 1+I * LDA + o_a, 1, ZERO, ref W, 1+I * LDW + o_w, 1);
 441:                          this._dgemv.Run("No transpose", N - I, I - 1,  - ONE, A, I + 1+1 * LDA + o_a, LDA
 442:                                          , W, 1+I * LDW + o_w, 1, ONE, ref W, I + 1+I * LDW + o_w, 1);
 443:                          this._dgemv.Run("Transpose", N - I, I - 1, ONE, A, I + 1+1 * LDA + o_a, LDA
 444:                                          , A, I + 1+I * LDA + o_a, 1, ZERO, ref W, 1+I * LDW + o_w, 1);
 445:                          this._dgemv.Run("No transpose", N - I, I - 1,  - ONE, W, I + 1+1 * LDW + o_w, LDW
 446:                                          , W, 1+I * LDW + o_w, 1, ONE, ref W, I + 1+I * LDW + o_w, 1);
 447:                          this._dscal.Run(N - I, TAU[I + o_tau], ref W, I + 1+I * LDW + o_w, 1);
 448:                          ALPHA =  - HALF * TAU[I + o_tau] * this._ddot.Run(N - I, W, I + 1+I * LDW + o_w, 1, A, I + 1+I * LDA + o_a, 1);
 449:                          this._daxpy.Run(N - I, ALPHA, A, I + 1+I * LDA + o_a, 1, ref W, I + 1+I * LDW + o_w, 1);
 450:                      }
 451:                      // *
 452:                  }
 453:              }
 454:              // *
 455:              return;
 456:              // *
 457:              // *     End of DLATRD
 458:              // *
 459:   
 460:              #endregion
 461:   
 462:          }
 463:      }
 464:  }