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:      /// DGETRI computes the inverse of a matrix using the LU factorization
  27:      /// computed by DGETRF.
  28:      /// 
  29:      /// This method inverts U and then computes inv(A) by solving the system
  30:      /// inv(A)*L = inv(U) for inv(A).
  31:      /// 
  32:      ///</summary>
  33:      public class DGETRI
  34:      {
  35:      
  36:   
  37:          #region Dependencies
  38:          
  39:          ILAENV _ilaenv; DGEMM _dgemm; DGEMV _dgemv; DSWAP _dswap; DTRSM _dtrsm; DTRTRI _dtrtri; XERBLA _xerbla; 
  40:   
  41:          #endregion
  42:   
  43:   
  44:          #region Fields
  45:          
  46:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; int I = 0; int IWS = 0; int J = 0; int JB = 0; 
  47:          int JJ = 0;int JP = 0; int LDWORK = 0; int LWKOPT = 0; int NB = 0; int NBMIN = 0; int NN = 0; 
  48:   
  49:          #endregion
  50:   
  51:          public DGETRI(ILAENV ilaenv, DGEMM dgemm, DGEMV dgemv, DSWAP dswap, DTRSM dtrsm, DTRTRI dtrtri, XERBLA xerbla)
  52:          {
  53:      
  54:   
  55:              #region Set Dependencies
  56:              
  57:              this._ilaenv = ilaenv; this._dgemm = dgemm; this._dgemv = dgemv; this._dswap = dswap; this._dtrsm = dtrsm; 
  58:              this._dtrtri = dtrtri;this._xerbla = xerbla; 
  59:   
  60:              #endregion
  61:   
  62:          }
  63:      
  64:          public DGETRI()
  65:          {
  66:      
  67:   
  68:              #region Dependencies (Initialization)
  69:              
  70:              IEEECK ieeeck = new IEEECK();
  71:              IPARMQ iparmq = new IPARMQ();
  72:              LSAME lsame = new LSAME();
  73:              XERBLA xerbla = new XERBLA();
  74:              DSWAP dswap = new DSWAP();
  75:              DSCAL dscal = new DSCAL();
  76:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  77:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  78:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  79:              DTRSM dtrsm = new DTRSM(lsame, xerbla);
  80:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  81:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  82:              DTRTI2 dtrti2 = new DTRTI2(lsame, dscal, dtrmv, xerbla);
  83:              DTRTRI dtrtri = new DTRTRI(lsame, ilaenv, dtrmm, dtrsm, dtrti2, xerbla);
  84:   
  85:              #endregion
  86:   
  87:   
  88:              #region Set Dependencies
  89:              
  90:              this._ilaenv = ilaenv; this._dgemm = dgemm; this._dgemv = dgemv; this._dswap = dswap; this._dtrsm = dtrsm; 
  91:              this._dtrtri = dtrtri;this._xerbla = xerbla; 
  92:   
  93:              #endregion
  94:   
  95:          }
  96:          /// <summary>
  97:          /// Purpose
  98:          /// =======
  99:          /// 
 100:          /// DGETRI computes the inverse of a matrix using the LU factorization
 101:          /// computed by DGETRF.
 102:          /// 
 103:          /// This method inverts U and then computes inv(A) by solving the system
 104:          /// inv(A)*L = inv(U) for inv(A).
 105:          /// 
 106:          ///</summary>
 107:          /// <param name="N">
 108:          /// (input) INTEGER
 109:          /// The order of 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 factors L and U from the factorization
 114:          /// A = P*L*U as computed by DGETRF.
 115:          /// On exit, if INFO = 0, the inverse of the original matrix A.
 116:          ///</param>
 117:          /// <param name="LDA">
 118:          /// (input) INTEGER
 119:          /// The leading dimension of the array A.  LDA .GE. max(1,N).
 120:          ///</param>
 121:          /// <param name="IPIV">
 122:          /// (input) INTEGER array, dimension (N)
 123:          /// The pivot indices from DGETRF; for 1.LE.i.LE.N, row i of the
 124:          /// matrix was interchanged with row IPIV(i).
 125:          ///</param>
 126:          /// <param name="WORK">
 127:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 128:          /// On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
 129:          ///</param>
 130:          /// <param name="LWORK">
 131:          /// (input) INTEGER
 132:          /// The dimension of the array WORK.  LWORK .GE. max(1,N).
 133:          /// For optimal performance LWORK .GE. N*NB, where NB is
 134:          /// the optimal blocksize returned by ILAENV.
 135:          /// 
 136:          /// If LWORK = -1, then a workspace query is assumed; the routine
 137:          /// only calculates the optimal size of the WORK array, returns
 138:          /// this value as the first entry of the WORK array, and no error
 139:          /// message related to LWORK is issued by XERBLA.
 140:          ///</param>
 141:          /// <param name="INFO">
 142:          /// (output) INTEGER
 143:          /// = 0:  successful exit
 144:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 145:          /// .GT. 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
 146:          /// singular and its inverse could not be computed.
 147:          ///</param>
 148:          public void Run(int N, ref double[] A, int offset_a, int LDA, int[] IPIV, int offset_ipiv, ref double[] WORK, int offset_work, int LWORK
 149:                           , ref int INFO)
 150:          {
 151:   
 152:              #region Array Index Correction
 153:              
 154:               int o_a = -1 - LDA + offset_a;  int o_ipiv = -1 + offset_ipiv;  int o_work = -1 + offset_work; 
 155:   
 156:              #endregion
 157:   
 158:   
 159:              #region Prolog
 160:              
 161:              // *
 162:              // *  -- LAPACK routine (version 3.1) --
 163:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 164:              // *     November 2006
 165:              // *
 166:              // *     .. Scalar Arguments ..
 167:              // *     ..
 168:              // *     .. Array Arguments ..
 169:              // *     ..
 170:              // *
 171:              // *  Purpose
 172:              // *  =======
 173:              // *
 174:              // *  DGETRI computes the inverse of a matrix using the LU factorization
 175:              // *  computed by DGETRF.
 176:              // *
 177:              // *  This method inverts U and then computes inv(A) by solving the system
 178:              // *  inv(A)*L = inv(U) for inv(A).
 179:              // *
 180:              // *  Arguments
 181:              // *  =========
 182:              // *
 183:              // *  N       (input) INTEGER
 184:              // *          The order of the matrix A.  N >= 0.
 185:              // *
 186:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 187:              // *          On entry, the factors L and U from the factorization
 188:              // *          A = P*L*U as computed by DGETRF.
 189:              // *          On exit, if INFO = 0, the inverse of the original matrix A.
 190:              // *
 191:              // *  LDA     (input) INTEGER
 192:              // *          The leading dimension of the array A.  LDA >= max(1,N).
 193:              // *
 194:              // *  IPIV    (input) INTEGER array, dimension (N)
 195:              // *          The pivot indices from DGETRF; for 1<=i<=N, row i of the
 196:              // *          matrix was interchanged with row IPIV(i).
 197:              // *
 198:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 199:              // *          On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
 200:              // *
 201:              // *  LWORK   (input) INTEGER
 202:              // *          The dimension of the array WORK.  LWORK >= max(1,N).
 203:              // *          For optimal performance LWORK >= N*NB, where NB is
 204:              // *          the optimal blocksize returned by ILAENV.
 205:              // *
 206:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 207:              // *          only calculates the optimal size of the WORK array, returns
 208:              // *          this value as the first entry of the WORK array, and no error
 209:              // *          message related to LWORK is issued by XERBLA.
 210:              // *
 211:              // *  INFO    (output) INTEGER
 212:              // *          = 0:  successful exit
 213:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 214:              // *          > 0:  if INFO = i, U(i,i) is exactly zero; the matrix is
 215:              // *                singular and its inverse could not be computed.
 216:              // *
 217:              // *  =====================================================================
 218:              // *
 219:              // *     .. Parameters ..
 220:              // *     ..
 221:              // *     .. Local Scalars ..
 222:              // *     ..
 223:              // *     .. External Functions ..
 224:              // *     ..
 225:              // *     .. External Subroutines ..
 226:              // *     ..
 227:              // *     .. Intrinsic Functions ..
 228:              //      INTRINSIC          MAX, MIN;
 229:              // *     ..
 230:              // *     .. Executable Statements ..
 231:              // *
 232:              // *     Test the input parameters.
 233:              // *
 234:   
 235:              #endregion
 236:   
 237:   
 238:              #region Body
 239:              
 240:              INFO = 0;
 241:              NB = this._ilaenv.Run(1, "DGETRI", " ", N,  - 1,  - 1,  - 1);
 242:              LWKOPT = N * NB;
 243:              WORK[1 + o_work] = LWKOPT;
 244:              LQUERY = (LWORK ==  - 1);
 245:              if (N < 0)
 246:              {
 247:                  INFO =  - 1;
 248:              }
 249:              else
 250:              {
 251:                  if (LDA < Math.Max(1, N))
 252:                  {
 253:                      INFO =  - 3;
 254:                  }
 255:                  else
 256:                  {
 257:                      if (LWORK < Math.Max(1, N) && !LQUERY)
 258:                      {
 259:                          INFO =  - 6;
 260:                      }
 261:                  }
 262:              }
 263:              if (INFO != 0)
 264:              {
 265:                  this._xerbla.Run("DGETRI",  - INFO);
 266:                  return;
 267:              }
 268:              else
 269:              {
 270:                  if (LQUERY)
 271:                  {
 272:                      return;
 273:                  }
 274:              }
 275:              // *
 276:              // *     Quick return if possible
 277:              // *
 278:              if (N == 0) return;
 279:              // *
 280:              // *     Form inv(U).  If INFO > 0 from DTRTRI, then U is singular,
 281:              // *     and the inverse is not computed.
 282:              // *
 283:              this._dtrtri.Run("Upper", "Non-unit", N, ref A, offset_a, LDA, ref INFO);
 284:              if (INFO > 0) return;
 285:              // *
 286:              NBMIN = 2;
 287:              LDWORK = N;
 288:              if (NB > 1 && NB < N)
 289:              {
 290:                  IWS = Math.Max(LDWORK * NB, 1);
 291:                  if (LWORK < IWS)
 292:                  {
 293:                      NB = LWORK / LDWORK;
 294:                      NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGETRI", " ", N,  - 1,  - 1,  - 1));
 295:                  }
 296:              }
 297:              else
 298:              {
 299:                  IWS = N;
 300:              }
 301:              // *
 302:              // *     Solve the equation inv(A)*L = inv(U) for inv(A).
 303:              // *
 304:              if (NB < NBMIN || NB >= N)
 305:              {
 306:                  // *
 307:                  // *        Use unblocked code.
 308:                  // *
 309:                  for (J = N; J >= 1; J +=  - 1)
 310:                  {
 311:                      // *
 312:                      // *           Copy current column of L to WORK and replace with zeros.
 313:                      // *
 314:                      for (I = J + 1; I <= N; I++)
 315:                      {
 316:                          WORK[I + o_work] = A[I+J * LDA + o_a];
 317:                          A[I+J * LDA + o_a] = ZERO;
 318:                      }
 319:                      // *
 320:                      // *           Compute current column of inv(A).
 321:                      // *
 322:                      if (J < N)
 323:                      {
 324:                          this._dgemv.Run("No transpose", N, N - J,  - ONE, A, 1+(J + 1) * LDA + o_a, LDA
 325:                                          , WORK, J + 1 + o_work, 1, ONE, ref A, 1+J * LDA + o_a, 1);
 326:                      }
 327:                  }
 328:              }
 329:              else
 330:              {
 331:                  // *
 332:                  // *        Use blocked code.
 333:                  // *
 334:                  NN = ((N - 1) / NB) * NB + 1;
 335:                  for (J = NN; ( - NB >= 0) ? (J <= 1) : (J >= 1); J +=  - NB)
 336:                  {
 337:                      JB = Math.Min(NB, N - J + 1);
 338:                      // *
 339:                      // *           Copy current block column of L to WORK and replace with
 340:                      // *           zeros.
 341:                      // *
 342:                      for (JJ = J; JJ <= J + JB - 1; JJ++)
 343:                      {
 344:                          for (I = JJ + 1; I <= N; I++)
 345:                          {
 346:                              WORK[I + (JJ - J) * LDWORK + o_work] = A[I+JJ * LDA + o_a];
 347:                              A[I+JJ * LDA + o_a] = ZERO;
 348:                          }
 349:                      }
 350:                      // *
 351:                      // *           Compute current block column of inv(A).
 352:                      // *
 353:                      if (J + JB <= N)
 354:                      {
 355:                          this._dgemm.Run("No transpose", "No transpose", N, JB, N - J - JB + 1,  - ONE
 356:                                          , A, 1+(J + JB) * LDA + o_a, LDA, WORK, J + JB + o_work, LDWORK, ONE, ref A, 1+J * LDA + o_a
 357:                                          , LDA);
 358:                      }
 359:                      this._dtrsm.Run("Right", "Lower", "No transpose", "Unit", N, JB
 360:                                      , ONE, WORK, J + o_work, LDWORK, ref A, 1+J * LDA + o_a, LDA);
 361:                  }
 362:              }
 363:              // *
 364:              // *     Apply column interchanges.
 365:              // *
 366:              for (J = N - 1; J >= 1; J +=  - 1)
 367:              {
 368:                  JP = IPIV[J + o_ipiv];
 369:                  if (JP != J) this._dswap.Run(N, ref A, 1+J * LDA + o_a, 1, ref A, 1+JP * LDA + o_a, 1);
 370:              }
 371:              // *
 372:              WORK[1 + o_work] = IWS;
 373:              return;
 374:              // *
 375:              // *     End of DGETRI
 376:              // *
 377:   
 378:              #endregion
 379:   
 380:          }
 381:      }
 382:  }