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