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 driver routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DSYEV computes all eigenvalues and, optionally, eigenvectors of a
  27:      /// real symmetric matrix A.
  28:      /// 
  29:      ///</summary>
  30:      public class DSYEV
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          LSAME _lsame; ILAENV _ilaenv; DLAMCH _dlamch; DLANSY _dlansy; DLASCL _dlascl; DORGTR _dorgtr; DSCAL _dscal; 
  37:          DSTEQR _dsteqr;DSTERF _dsterf; DSYTRD _dsytrd; XERBLA _xerbla; 
  38:   
  39:          #endregion
  40:   
  41:   
  42:          #region Fields
  43:          
  44:          const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LOWER = false; bool LQUERY = false; bool WANTZ = false; 
  45:          int IINFO = 0;int IMAX = 0; int INDE = 0; int INDTAU = 0; int INDWRK = 0; int ISCALE = 0; int LLWORK = 0; int LWKOPT = 0; 
  46:          int NB = 0;double ANRM = 0; double BIGNUM = 0; double EPS = 0; double RMAX = 0; double RMIN = 0; double SAFMIN = 0; 
  47:          double SIGMA = 0;double SMLNUM = 0; 
  48:   
  49:          #endregion
  50:   
  51:          public DSYEV(LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANSY dlansy, DLASCL dlascl, DORGTR dorgtr, DSCAL dscal, DSTEQR dsteqr, DSTERF dsterf, DSYTRD dsytrd
  52:                       , XERBLA xerbla)
  53:          {
  54:      
  55:   
  56:              #region Set Dependencies
  57:              
  58:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlansy = dlansy; this._dlascl = dlascl; 
  59:              this._dorgtr = dorgtr;this._dscal = dscal; this._dsteqr = dsteqr; this._dsterf = dsterf; this._dsytrd = dsytrd; 
  60:              this._xerbla = xerbla;
  61:   
  62:              #endregion
  63:   
  64:          }
  65:      
  66:          public DSYEV()
  67:          {
  68:      
  69:   
  70:              #region Dependencies (Initialization)
  71:              
  72:              LSAME lsame = new LSAME();
  73:              IEEECK ieeeck = new IEEECK();
  74:              IPARMQ iparmq = new IPARMQ();
  75:              DLAMC3 dlamc3 = new DLAMC3();
  76:              DLASSQ dlassq = new DLASSQ();
  77:              XERBLA xerbla = new XERBLA();
  78:              DCOPY dcopy = new DCOPY();
  79:              DSCAL dscal = new DSCAL();
  80:              DLAPY2 dlapy2 = new DLAPY2();
  81:              DLAE2 dlae2 = new DLAE2();
  82:              DLAEV2 dlaev2 = new DLAEV2();
  83:              DSWAP dswap = new DSWAP();
  84:              DAXPY daxpy = new DAXPY();
  85:              DNRM2 dnrm2 = new DNRM2();
  86:              DDOT ddot = new DDOT();
  87:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  88:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  89:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  90:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  91:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  92:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  93:              DLANSY dlansy = new DLANSY(dlassq, lsame);
  94:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
  95:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  96:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  97:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
  98:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  99:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
 100:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
 101:              DGER dger = new DGER(xerbla);
 102:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 103:              DORG2L dorg2l = new DORG2L(dlarf, dscal, xerbla);
 104:              DORGQL dorgql = new DORGQL(dlarfb, dlarft, dorg2l, xerbla, ilaenv);
 105:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
 106:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
 107:              DORGTR dorgtr = new DORGTR(lsame, ilaenv, dorgql, dorgqr, xerbla);
 108:              DLANST dlanst = new DLANST(lsame, dlassq);
 109:              DLARTG dlartg = new DLARTG(dlamch);
 110:              DLASET dlaset = new DLASET(lsame);
 111:              DLASR dlasr = new DLASR(lsame, xerbla);
 112:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
 113:              DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
 114:                                         , dlasrt, dswap, xerbla);
 115:              DSTERF dsterf = new DSTERF(dlamch, dlanst, dlapy2, dlae2, dlascl, dlasrt, xerbla);
 116:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 117:              DSYMV dsymv = new DSYMV(lsame, xerbla);
 118:              DLATRD dlatrd = new DLATRD(daxpy, dgemv, dlarfg, dscal, dsymv, lsame, ddot);
 119:              DSYR2K dsyr2k = new DSYR2K(lsame, xerbla);
 120:              DSYR2 dsyr2 = new DSYR2(lsame, xerbla);
 121:              DSYTD2 dsytd2 = new DSYTD2(daxpy, dlarfg, dsymv, dsyr2, xerbla, lsame, ddot);
 122:              DSYTRD dsytrd = new DSYTRD(dlatrd, dsyr2k, dsytd2, xerbla, lsame, ilaenv);
 123:   
 124:              #endregion
 125:   
 126:   
 127:              #region Set Dependencies
 128:              
 129:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlansy = dlansy; this._dlascl = dlascl; 
 130:              this._dorgtr = dorgtr;this._dscal = dscal; this._dsteqr = dsteqr; this._dsterf = dsterf; this._dsytrd = dsytrd; 
 131:              this._xerbla = xerbla;
 132:   
 133:              #endregion
 134:   
 135:          }
 136:          /// <summary>
 137:          /// Purpose
 138:          /// =======
 139:          /// 
 140:          /// DSYEV computes all eigenvalues and, optionally, eigenvectors of a
 141:          /// real symmetric matrix A.
 142:          /// 
 143:          ///</summary>
 144:          /// <param name="JOBZ">
 145:          /// (input) CHARACTER*1
 146:          /// = 'N':  Compute eigenvalues only;
 147:          /// = 'V':  Compute eigenvalues and eigenvectors.
 148:          ///</param>
 149:          /// <param name="UPLO">
 150:          /// (input) CHARACTER*1
 151:          /// = 'U':  Upper triangle of A is stored;
 152:          /// = 'L':  Lower triangle of A is stored.
 153:          ///</param>
 154:          /// <param name="N">
 155:          /// (input) INTEGER
 156:          /// The order of the matrix A.  N .GE. 0.
 157:          ///</param>
 158:          /// <param name="A">
 159:          /// (input/output) DOUBLE PRECISION array, dimension (LDA, N)
 160:          /// On entry, the symmetric matrix A.  If UPLO = 'U', the
 161:          /// leading N-by-N upper triangular part of A contains the
 162:          /// upper triangular part of the matrix A.  If UPLO = 'L',
 163:          /// the leading N-by-N lower triangular part of A contains
 164:          /// the lower triangular part of the matrix A.
 165:          /// On exit, if JOBZ = 'V', then if INFO = 0, A contains the
 166:          /// orthonormal eigenvectors of the matrix A.
 167:          /// If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
 168:          /// or the upper triangle (if UPLO='U') of A, including the
 169:          /// diagonal, is destroyed.
 170:          ///</param>
 171:          /// <param name="LDA">
 172:          /// (input) INTEGER
 173:          /// The leading dimension of the array A.  LDA .GE. max(1,N).
 174:          ///</param>
 175:          /// <param name="W">
 176:          /// (output) DOUBLE PRECISION array, dimension (N)
 177:          /// If INFO = 0, the eigenvalues in ascending order.
 178:          ///</param>
 179:          /// <param name="WORK">
 180:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 181:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 182:          ///</param>
 183:          /// <param name="LWORK">
 184:          /// (input) INTEGER
 185:          /// The length of the array WORK.  LWORK .GE. max(1,3*N-1).
 186:          /// For optimal efficiency, LWORK .GE. (NB+2)*N,
 187:          /// where NB is the blocksize for DSYTRD returned by ILAENV.
 188:          /// 
 189:          /// If LWORK = -1, then a workspace query is assumed; the routine
 190:          /// only calculates the optimal size of the WORK array, returns
 191:          /// this value as the first entry of the WORK array, and no error
 192:          /// message related to LWORK is issued by XERBLA.
 193:          ///</param>
 194:          /// <param name="INFO">
 195:          /// (output) INTEGER
 196:          /// = 0:  successful exit
 197:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 198:          /// .GT. 0:  if INFO = i, the algorithm failed to converge; i
 199:          /// off-diagonal elements of an intermediate tridiagonal
 200:          /// form did not converge to zero.
 201:          ///</param>
 202:          public void Run(string JOBZ, string UPLO, int N, ref double[] A, int offset_a, int LDA, ref double[] W, int offset_w
 203:                           , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
 204:          {
 205:   
 206:              #region Array Index Correction
 207:              
 208:               int o_a = -1 - LDA + offset_a;  int o_w = -1 + offset_w;  int o_work = -1 + offset_work; 
 209:   
 210:              #endregion
 211:   
 212:   
 213:              #region Strings
 214:              
 215:              JOBZ = JOBZ.Substring(0, 1);  UPLO = UPLO.Substring(0, 1);  
 216:   
 217:              #endregion
 218:   
 219:   
 220:              #region Prolog
 221:              
 222:              // *
 223:              // *  -- LAPACK driver routine (version 3.1) --
 224:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 225:              // *     November 2006
 226:              // *
 227:              // *     .. Scalar Arguments ..
 228:              // *     ..
 229:              // *     .. Array Arguments ..
 230:              // *     ..
 231:              // *
 232:              // *  Purpose
 233:              // *  =======
 234:              // *
 235:              // *  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
 236:              // *  real symmetric matrix A.
 237:              // *
 238:              // *  Arguments
 239:              // *  =========
 240:              // *
 241:              // *  JOBZ    (input) CHARACTER*1
 242:              // *          = 'N':  Compute eigenvalues only;
 243:              // *          = 'V':  Compute eigenvalues and eigenvectors.
 244:              // *
 245:              // *  UPLO    (input) CHARACTER*1
 246:              // *          = 'U':  Upper triangle of A is stored;
 247:              // *          = 'L':  Lower triangle of A is stored.
 248:              // *
 249:              // *  N       (input) INTEGER
 250:              // *          The order of the matrix A.  N >= 0.
 251:              // *
 252:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
 253:              // *          On entry, the symmetric matrix A.  If UPLO = 'U', the
 254:              // *          leading N-by-N upper triangular part of A contains the
 255:              // *          upper triangular part of the matrix A.  If UPLO = 'L',
 256:              // *          the leading N-by-N lower triangular part of A contains
 257:              // *          the lower triangular part of the matrix A.
 258:              // *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
 259:              // *          orthonormal eigenvectors of the matrix A.
 260:              // *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
 261:              // *          or the upper triangle (if UPLO='U') of A, including the
 262:              // *          diagonal, is destroyed.
 263:              // *
 264:              // *  LDA     (input) INTEGER
 265:              // *          The leading dimension of the array A.  LDA >= max(1,N).
 266:              // *
 267:              // *  W       (output) DOUBLE PRECISION array, dimension (N)
 268:              // *          If INFO = 0, the eigenvalues in ascending order.
 269:              // *
 270:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 271:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 272:              // *
 273:              // *  LWORK   (input) INTEGER
 274:              // *          The length of the array WORK.  LWORK >= max(1,3*N-1).
 275:              // *          For optimal efficiency, LWORK >= (NB+2)*N,
 276:              // *          where NB is the blocksize for DSYTRD returned by ILAENV.
 277:              // *
 278:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 279:              // *          only calculates the optimal size of the WORK array, returns
 280:              // *          this value as the first entry of the WORK array, and no error
 281:              // *          message related to LWORK is issued by XERBLA.
 282:              // *
 283:              // *  INFO    (output) INTEGER
 284:              // *          = 0:  successful exit
 285:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 286:              // *          > 0:  if INFO = i, the algorithm failed to converge; i
 287:              // *                off-diagonal elements of an intermediate tridiagonal
 288:              // *                form did not converge to zero.
 289:              // *
 290:              // *  =====================================================================
 291:              // *
 292:              // *     .. Parameters ..
 293:              // *     ..
 294:              // *     .. Local Scalars ..
 295:              // *     ..
 296:              // *     .. External Functions ..
 297:              // *     ..
 298:              // *     .. External Subroutines ..
 299:              // *     ..
 300:              // *     .. Intrinsic Functions ..
 301:              //      INTRINSIC          MAX, SQRT;
 302:              // *     ..
 303:              // *     .. Executable Statements ..
 304:              // *
 305:              // *     Test the input parameters.
 306:              // *
 307:   
 308:              #endregion
 309:   
 310:   
 311:              #region Body
 312:              
 313:              WANTZ = this._lsame.Run(JOBZ, "V");
 314:              LOWER = this._lsame.Run(UPLO, "L");
 315:              LQUERY = (LWORK ==  - 1);
 316:              // *
 317:              INFO = 0;
 318:              if (!(WANTZ || this._lsame.Run(JOBZ, "N")))
 319:              {
 320:                  INFO =  - 1;
 321:              }
 322:              else
 323:              {
 324:                  if (!(LOWER || this._lsame.Run(UPLO, "U")))
 325:                  {
 326:                      INFO =  - 2;
 327:                  }
 328:                  else
 329:                  {
 330:                      if (N < 0)
 331:                      {
 332:                          INFO =  - 3;
 333:                      }
 334:                      else
 335:                      {
 336:                          if (LDA < Math.Max(1, N))
 337:                          {
 338:                              INFO =  - 5;
 339:                          }
 340:                      }
 341:                  }
 342:              }
 343:              // *
 344:              if (INFO == 0)
 345:              {
 346:                  NB = this._ilaenv.Run(1, "DSYTRD", UPLO, N,  - 1,  - 1,  - 1);
 347:                  LWKOPT = Math.Max(1, (NB + 2) * N);
 348:                  WORK[1 + o_work] = LWKOPT;
 349:                  // *
 350:                  if (LWORK < Math.Max(1, 3 * N - 1) && !LQUERY) INFO =  - 8;
 351:              }
 352:              // *
 353:              if (INFO != 0)
 354:              {
 355:                  this._xerbla.Run("DSYEV ",  - INFO);
 356:                  return;
 357:              }
 358:              else
 359:              {
 360:                  if (LQUERY)
 361:                  {
 362:                      return;
 363:                  }
 364:              }
 365:              // *
 366:              // *     Quick return if possible
 367:              // *
 368:              if (N == 0)
 369:              {
 370:                  return;
 371:              }
 372:              // *
 373:              if (N == 1)
 374:              {
 375:                  W[1 + o_w] = A[1+1 * LDA + o_a];
 376:                  WORK[1 + o_work] = 2;
 377:                  if (WANTZ) A[1+1 * LDA + o_a] = ONE;
 378:                  return;
 379:              }
 380:              // *
 381:              // *     Get machine constants.
 382:              // *
 383:              SAFMIN = this._dlamch.Run("Safe minimum");
 384:              EPS = this._dlamch.Run("Precision");
 385:              SMLNUM = SAFMIN / EPS;
 386:              BIGNUM = ONE / SMLNUM;
 387:              RMIN = Math.Sqrt(SMLNUM);
 388:              RMAX = Math.Sqrt(BIGNUM);
 389:              // *
 390:              // *     Scale matrix to allowable range, if necessary.
 391:              // *
 392:              ANRM = this._dlansy.Run("M", UPLO, N, A, offset_a, LDA, ref WORK, offset_work);
 393:              ISCALE = 0;
 394:              if (ANRM > ZERO && ANRM < RMIN)
 395:              {
 396:                  ISCALE = 1;
 397:                  SIGMA = RMIN / ANRM;
 398:              }
 399:              else
 400:              {
 401:                  if (ANRM > RMAX)
 402:                  {
 403:                      ISCALE = 1;
 404:                      SIGMA = RMAX / ANRM;
 405:                  }
 406:              }
 407:              if (ISCALE == 1)
 408:              {
 409:                  this._dlascl.Run(UPLO, 0, 0, ONE, SIGMA, N
 410:                                   , N, ref A, offset_a, LDA, ref INFO);
 411:              }
 412:              // *
 413:              // *     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
 414:              // *
 415:              INDE = 1;
 416:              INDTAU = INDE + N;
 417:              INDWRK = INDTAU + N;
 418:              LLWORK = LWORK - INDWRK + 1;
 419:              this._dsytrd.Run(UPLO, N, ref A, offset_a, LDA, ref W, offset_w, ref WORK, INDE + o_work
 420:                               , ref WORK, INDTAU + o_work, ref WORK, INDWRK + o_work, LLWORK, ref IINFO);
 421:              // *
 422:              // *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
 423:              // *     DORGTR to generate the orthogonal matrix, then call DSTEQR.
 424:              // *
 425:              if (!WANTZ)
 426:              {
 427:                  this._dsterf.Run(N, ref W, offset_w, ref WORK, INDE + o_work, ref INFO);
 428:              }
 429:              else
 430:              {
 431:                  this._dorgtr.Run(UPLO, N, ref A, offset_a, LDA, WORK, INDTAU + o_work, ref WORK, INDWRK + o_work
 432:                                   , LLWORK, ref IINFO);
 433:                  this._dsteqr.Run(JOBZ, N, ref W, offset_w, ref WORK, INDE + o_work, ref A, offset_a, LDA
 434:                                   , ref WORK, INDTAU + o_work, ref INFO);
 435:              }
 436:              // *
 437:              // *     If matrix was scaled, then rescale eigenvalues appropriately.
 438:              // *
 439:              if (ISCALE == 1)
 440:              {
 441:                  if (INFO == 0)
 442:                  {
 443:                      IMAX = N;
 444:                  }
 445:                  else
 446:                  {
 447:                      IMAX = INFO - 1;
 448:                  }
 449:                  this._dscal.Run(IMAX, ONE / SIGMA, ref W, offset_w, 1);
 450:              }
 451:              // *
 452:              // *     Set WORK(1) to optimal workspace size.
 453:              // *
 454:              WORK[1 + o_work] = LWKOPT;
 455:              // *
 456:              return;
 457:              // *
 458:              // *     End of DSYEV
 459:              // *
 460:   
 461:              #endregion
 462:   
 463:          }
 464:      }
 465:  }