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:      /// DGETF2 computes an LU factorization of a general m-by-n matrix A
  27:      /// using partial pivoting with row interchanges.
  28:      /// 
  29:      /// The factorization has the form
  30:      /// A = P * L * U
  31:      /// where P is a permutation matrix, L is lower triangular with unit
  32:      /// diagonal elements (lower trapezoidal if m .GT. n), and U is upper
  33:      /// triangular (upper trapezoidal if m .LT. n).
  34:      /// 
  35:      /// This is the right-looking Level 2 BLAS version of the algorithm.
  36:      /// 
  37:      ///</summary>
  38:      public class DGETF2
  39:      {
  40:      
  41:   
  42:          #region Dependencies
  43:          
  44:          DLAMCH _dlamch; IDAMAX _idamax; DGER _dger; DSCAL _dscal; DSWAP _dswap; XERBLA _xerbla; 
  45:   
  46:          #endregion
  47:   
  48:   
  49:          #region Fields
  50:          
  51:          const double ONE = 1.0E+0; const double ZERO = 0.0E+0; double SFMIN = 0; int I = 0; int J = 0; int JP = 0; 
  52:   
  53:          #endregion
  54:   
  55:          public DGETF2(DLAMCH dlamch, IDAMAX idamax, DGER dger, DSCAL dscal, DSWAP dswap, XERBLA xerbla)
  56:          {
  57:      
  58:   
  59:              #region Set Dependencies
  60:              
  61:              this._dlamch = dlamch; this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap; 
  62:              this._xerbla = xerbla;
  63:   
  64:              #endregion
  65:   
  66:          }
  67:      
  68:          public DGETF2()
  69:          {
  70:      
  71:   
  72:              #region Dependencies (Initialization)
  73:              
  74:              LSAME lsame = new LSAME();
  75:              DLAMC3 dlamc3 = new DLAMC3();
  76:              IDAMAX idamax = new IDAMAX();
  77:              XERBLA xerbla = new XERBLA();
  78:              DSCAL dscal = new DSCAL();
  79:              DSWAP dswap = new DSWAP();
  80:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  81:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  82:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  83:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  84:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  85:              DGER dger = new DGER(xerbla);
  86:   
  87:              #endregion
  88:   
  89:   
  90:              #region Set Dependencies
  91:              
  92:              this._dlamch = dlamch; this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap; 
  93:              this._xerbla = xerbla;
  94:   
  95:              #endregion
  96:   
  97:          }
  98:          /// <summary>
  99:          /// Purpose
 100:          /// =======
 101:          /// 
 102:          /// DGETF2 computes an LU factorization of a general m-by-n matrix A
 103:          /// using partial pivoting with row interchanges.
 104:          /// 
 105:          /// The factorization has the form
 106:          /// A = P * L * U
 107:          /// where P is a permutation matrix, L is lower triangular with unit
 108:          /// diagonal elements (lower trapezoidal if m .GT. n), and U is upper
 109:          /// triangular (upper trapezoidal if m .LT. n).
 110:          /// 
 111:          /// This is the right-looking Level 2 BLAS version of the algorithm.
 112:          /// 
 113:          ///</summary>
 114:          /// <param name="M">
 115:          /// (input) INTEGER
 116:          /// The number of rows of the matrix A.  M .GE. 0.
 117:          ///</param>
 118:          /// <param name="N">
 119:          /// (input) INTEGER
 120:          /// The number of columns of the matrix A.  N .GE. 0.
 121:          ///</param>
 122:          /// <param name="A">
 123:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 124:          /// On entry, the m by n matrix to be factored.
 125:          /// On exit, the factors L and U from the factorization
 126:          /// A = P*L*U; the unit diagonal elements of L are not stored.
 127:          ///</param>
 128:          /// <param name="LDA">
 129:          /// (input) INTEGER
 130:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 131:          ///</param>
 132:          /// <param name="IPIV">
 133:          /// (output) INTEGER array, dimension (min(M,N))
 134:          /// The pivot indices; for 1 .LE. i .LE. min(M,N), row i of the
 135:          /// matrix was interchanged with row IPIV(i).
 136:          ///</param>
 137:          /// <param name="INFO">
 138:          /// (output) INTEGER
 139:          /// = 0: successful exit
 140:          /// .LT. 0: if INFO = -k, the k-th argument had an illegal value
 141:          /// .GT. 0: if INFO = k, U(k,k) is exactly zero. The factorization
 142:          /// has been completed, but the factor U is exactly
 143:          /// singular, and division by zero will occur if it is used
 144:          /// to solve a system of equations.
 145:          ///</param>
 146:          public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref int[] IPIV, int offset_ipiv, ref int INFO)
 147:          {
 148:   
 149:              #region Array Index Correction
 150:              
 151:               int o_a = -1 - LDA + offset_a;  int o_ipiv = -1 + offset_ipiv; 
 152:   
 153:              #endregion
 154:   
 155:   
 156:              #region Prolog
 157:              
 158:              // *
 159:              // *  -- LAPACK routine (version 3.1) --
 160:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 161:              // *     November 2006
 162:              // *
 163:              // *     .. Scalar Arguments ..
 164:              // *     ..
 165:              // *     .. Array Arguments ..
 166:              // *     ..
 167:              // *
 168:              // *  Purpose
 169:              // *  =======
 170:              // *
 171:              // *  DGETF2 computes an LU factorization of a general m-by-n matrix A
 172:              // *  using partial pivoting with row interchanges.
 173:              // *
 174:              // *  The factorization has the form
 175:              // *     A = P * L * U
 176:              // *  where P is a permutation matrix, L is lower triangular with unit
 177:              // *  diagonal elements (lower trapezoidal if m > n), and U is upper
 178:              // *  triangular (upper trapezoidal if m < n).
 179:              // *
 180:              // *  This is the right-looking Level 2 BLAS version of the algorithm.
 181:              // *
 182:              // *  Arguments
 183:              // *  =========
 184:              // *
 185:              // *  M       (input) INTEGER
 186:              // *          The number of rows of the matrix A.  M >= 0.
 187:              // *
 188:              // *  N       (input) INTEGER
 189:              // *          The number of columns of the matrix A.  N >= 0.
 190:              // *
 191:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 192:              // *          On entry, the m by n matrix to be factored.
 193:              // *          On exit, the factors L and U from the factorization
 194:              // *          A = P*L*U; the unit diagonal elements of L are not stored.
 195:              // *
 196:              // *  LDA     (input) INTEGER
 197:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 198:              // *
 199:              // *  IPIV    (output) INTEGER array, dimension (min(M,N))
 200:              // *          The pivot indices; for 1 <= i <= min(M,N), row i of the
 201:              // *          matrix was interchanged with row IPIV(i).
 202:              // *
 203:              // *  INFO    (output) INTEGER
 204:              // *          = 0: successful exit
 205:              // *          < 0: if INFO = -k, the k-th argument had an illegal value
 206:              // *          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
 207:              // *               has been completed, but the factor U is exactly
 208:              // *               singular, and division by zero will occur if it is used
 209:              // *               to solve a system of equations.
 210:              // *
 211:              // *  =====================================================================
 212:              // *
 213:              // *     .. Parameters ..
 214:              // *     ..
 215:              // *     .. Local Scalars ..
 216:              // *     ..
 217:              // *     .. External Functions ..
 218:              // *     ..
 219:              // *     .. External Subroutines ..
 220:              // *     ..
 221:              // *     .. Intrinsic Functions ..
 222:              //      INTRINSIC          MAX, MIN;
 223:              // *     ..
 224:              // *     .. Executable Statements ..
 225:              // *
 226:              // *     Test the input parameters.
 227:              // *
 228:   
 229:              #endregion
 230:   
 231:   
 232:              #region Body
 233:              
 234:              INFO = 0;
 235:              if (M < 0)
 236:              {
 237:                  INFO =  - 1;
 238:              }
 239:              else
 240:              {
 241:                  if (N < 0)
 242:                  {
 243:                      INFO =  - 2;
 244:                  }
 245:                  else
 246:                  {
 247:                      if (LDA < Math.Max(1, M))
 248:                      {
 249:                          INFO =  - 4;
 250:                      }
 251:                  }
 252:              }
 253:              if (INFO != 0)
 254:              {
 255:                  this._xerbla.Run("DGETF2",  - INFO);
 256:                  return;
 257:              }
 258:              // *
 259:              // *     Quick return if possible
 260:              // *
 261:              if (M == 0 || N == 0) return;
 262:              // *
 263:              // *     Compute machine safe minimum 
 264:              // * 
 265:              SFMIN = this._dlamch.Run("S");
 266:              // *
 267:              for (J = 1; J <= Math.Min(M, N); J++)
 268:              {
 269:                  // *
 270:                  // *        Find pivot and test for singularity.
 271:                  // *
 272:                  JP = J - 1 + this._idamax.Run(M - J + 1, A, J+J * LDA + o_a, 1);
 273:                  IPIV[J + o_ipiv] = JP;
 274:                  if (A[JP+J * LDA + o_a] != ZERO)
 275:                  {
 276:                      // *
 277:                      // *           Apply the interchange to columns 1:N.
 278:                      // *
 279:                      if (JP != J) this._dswap.Run(N, ref A, J+1 * LDA + o_a, LDA, ref A, JP+1 * LDA + o_a, LDA);
 280:                      // *
 281:                      // *           Compute elements J+1:M of J-th column.
 282:                      // *
 283:                      if (J < M)
 284:                      {
 285:                          if (Math.Abs(A[J+J * LDA + o_a]) >= SFMIN)
 286:                          {
 287:                              this._dscal.Run(M - J, ONE / A[J+J * LDA + o_a], ref A, J + 1+J * LDA + o_a, 1);
 288:                          }
 289:                          else
 290:                          {
 291:                              for (I = 1; I <= M - J; I++)
 292:                              {
 293:                                  A[J + I+J * LDA + o_a] = A[J + I+J * LDA + o_a] / A[J+J * LDA + o_a];
 294:                              }
 295:                          }
 296:                      }
 297:                      // *
 298:                  }
 299:                  else
 300:                  {
 301:                      if (INFO == 0)
 302:                      {
 303:                          // *
 304:                          INFO = J;
 305:                      }
 306:                  }
 307:                  // *
 308:                  if (J < Math.Min(M, N))
 309:                  {
 310:                      // *
 311:                      // *           Update trailing submatrix.
 312:                      // *
 313:                      this._dger.Run(M - J, N - J,  - ONE, A, J + 1+J * LDA + o_a, 1, A, J+(J + 1) * LDA + o_a
 314:                                     , LDA, ref A, J + 1+(J + 1) * LDA + o_a, LDA);
 315:                  }
 316:              }
 317:              return;
 318:              // *
 319:              // *     End of DGETF2
 320:              // *
 321:   
 322:              #endregion
 323:   
 324:          }
 325:      }
 326:  }