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:      /// DGBTF2 computes an LU factorization of a real m-by-n band matrix A
  27:      /// using partial pivoting with row interchanges.
  28:      /// 
  29:      /// This is the unblocked version of the algorithm, calling Level 2 BLAS.
  30:      /// 
  31:      ///</summary>
  32:      public class DGBTF2
  33:      {
  34:      
  35:   
  36:          #region Dependencies
  37:          
  38:          IDAMAX _idamax; DGER _dger; DSCAL _dscal; DSWAP _dswap; XERBLA _xerbla; 
  39:   
  40:          #endregion
  41:   
  42:   
  43:          #region Fields
  44:          
  45:          const double ONE = 1.0E+0; const double ZERO = 0.0E+0; int I = 0; int J = 0; int JP = 0; int JU = 0; int KM = 0; 
  46:          int KV = 0;
  47:   
  48:          #endregion
  49:   
  50:          public DGBTF2(IDAMAX idamax, DGER dger, DSCAL dscal, DSWAP dswap, XERBLA xerbla)
  51:          {
  52:      
  53:   
  54:              #region Set Dependencies
  55:              
  56:              this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap; this._xerbla = xerbla; 
  57:   
  58:              #endregion
  59:   
  60:          }
  61:      
  62:          public DGBTF2()
  63:          {
  64:      
  65:   
  66:              #region Dependencies (Initialization)
  67:              
  68:              IDAMAX idamax = new IDAMAX();
  69:              XERBLA xerbla = new XERBLA();
  70:              DSCAL dscal = new DSCAL();
  71:              DSWAP dswap = new DSWAP();
  72:              DGER dger = new DGER(xerbla);
  73:   
  74:              #endregion
  75:   
  76:   
  77:              #region Set Dependencies
  78:              
  79:              this._idamax = idamax; this._dger = dger; this._dscal = dscal; this._dswap = dswap; this._xerbla = xerbla; 
  80:   
  81:              #endregion
  82:   
  83:          }
  84:          /// <summary>
  85:          /// Purpose
  86:          /// =======
  87:          /// 
  88:          /// DGBTF2 computes an LU factorization of a real m-by-n band matrix A
  89:          /// using partial pivoting with row interchanges.
  90:          /// 
  91:          /// This is the unblocked version of the algorithm, calling Level 2 BLAS.
  92:          /// 
  93:          ///</summary>
  94:          /// <param name="M">
  95:          /// (input) INTEGER
  96:          /// The number of rows of the matrix A.  M .GE. 0.
  97:          ///</param>
  98:          /// <param name="N">
  99:          /// (input) INTEGER
 100:          /// The number of columns of the matrix A.  N .GE. 0.
 101:          ///</param>
 102:          /// <param name="KL">
 103:          /// (input) INTEGER
 104:          /// The number of subdiagonals within the band of A.  KL .GE. 0.
 105:          ///</param>
 106:          /// <param name="KU">
 107:          /// (input) INTEGER
 108:          /// The number of superdiagonals within the band of A.  KU .GE. 0.
 109:          ///</param>
 110:          /// <param name="AB">
 111:          /// (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 112:          /// On entry, the matrix A in band storage, in rows KL+1 to
 113:          /// 2*KL+KU+1; rows 1 to KL of the array need not be set.
 114:          /// The j-th column of A is stored in the j-th column of the
 115:          /// array AB as follows:
 116:          /// AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku).LE.i.LE.min(m,j+kl)
 117:          /// 
 118:          /// On exit, details of the factorization: U is stored as an
 119:          /// upper triangular band matrix with KL+KU superdiagonals in
 120:          /// rows 1 to KL+KU+1, and the multipliers used during the
 121:          /// factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
 122:          /// See below for further details.
 123:          ///</param>
 124:          /// <param name="LDAB">
 125:          /// (input) INTEGER
 126:          /// The leading dimension of the array AB.  LDAB .GE. 2*KL+KU+1.
 127:          ///</param>
 128:          /// <param name="IPIV">
 129:          /// (output) INTEGER array, dimension (min(M,N))
 130:          /// The pivot indices; for 1 .LE. i .LE. min(M,N), row i of the
 131:          /// matrix was interchanged with row IPIV(i).
 132:          ///</param>
 133:          /// <param name="INFO">
 134:          /// (output) INTEGER
 135:          /// = 0: successful exit
 136:          /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
 137:          /// .GT. 0: if INFO = +i, U(i,i) is exactly zero. The factorization
 138:          /// has been completed, but the factor U is exactly
 139:          /// singular, and division by zero will occur if it is used
 140:          /// to solve a system of equations.
 141:          ///</param>
 142:          public void Run(int M, int N, int KL, int KU, ref double[] AB, int offset_ab, int LDAB
 143:                           , ref int[] IPIV, int offset_ipiv, ref int INFO)
 144:          {
 145:   
 146:              #region Array Index Correction
 147:              
 148:               int o_ab = -1 - LDAB + offset_ab;  int o_ipiv = -1 + offset_ipiv; 
 149:   
 150:              #endregion
 151:   
 152:   
 153:              #region Prolog
 154:              
 155:              // *
 156:              // *  -- LAPACK routine (version 3.1) --
 157:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 158:              // *     November 2006
 159:              // *
 160:              // *     .. Scalar Arguments ..
 161:              // *     ..
 162:              // *     .. Array Arguments ..
 163:              // *     ..
 164:              // *
 165:              // *  Purpose
 166:              // *  =======
 167:              // *
 168:              // *  DGBTF2 computes an LU factorization of a real m-by-n band matrix A
 169:              // *  using partial pivoting with row interchanges.
 170:              // *
 171:              // *  This is the unblocked version of the algorithm, calling Level 2 BLAS.
 172:              // *
 173:              // *  Arguments
 174:              // *  =========
 175:              // *
 176:              // *  M       (input) INTEGER
 177:              // *          The number of rows of the matrix A.  M >= 0.
 178:              // *
 179:              // *  N       (input) INTEGER
 180:              // *          The number of columns of the matrix A.  N >= 0.
 181:              // *
 182:              // *  KL      (input) INTEGER
 183:              // *          The number of subdiagonals within the band of A.  KL >= 0.
 184:              // *
 185:              // *  KU      (input) INTEGER
 186:              // *          The number of superdiagonals within the band of A.  KU >= 0.
 187:              // *
 188:              // *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 189:              // *          On entry, the matrix A in band storage, in rows KL+1 to
 190:              // *          2*KL+KU+1; rows 1 to KL of the array need not be set.
 191:              // *          The j-th column of A is stored in the j-th column of the
 192:              // *          array AB as follows:
 193:              // *          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
 194:              // *
 195:              // *          On exit, details of the factorization: U is stored as an
 196:              // *          upper triangular band matrix with KL+KU superdiagonals in
 197:              // *          rows 1 to KL+KU+1, and the multipliers used during the
 198:              // *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
 199:              // *          See below for further details.
 200:              // *
 201:              // *  LDAB    (input) INTEGER
 202:              // *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
 203:              // *
 204:              // *  IPIV    (output) INTEGER array, dimension (min(M,N))
 205:              // *          The pivot indices; for 1 <= i <= min(M,N), row i of the
 206:              // *          matrix was interchanged with row IPIV(i).
 207:              // *
 208:              // *  INFO    (output) INTEGER
 209:              // *          = 0: successful exit
 210:              // *          < 0: if INFO = -i, the i-th argument had an illegal value
 211:              // *          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
 212:              // *               has been completed, but the factor U is exactly
 213:              // *               singular, and division by zero will occur if it is used
 214:              // *               to solve a system of equations.
 215:              // *
 216:              // *  Further Details
 217:              // *  ===============
 218:              // *
 219:              // *  The band storage scheme is illustrated by the following example, when
 220:              // *  M = N = 6, KL = 2, KU = 1:
 221:              // *
 222:              // *  On entry:                       On exit:
 223:              // *
 224:              // *      *    *    *    +    +    +       *    *    *   u14  u25  u36
 225:              // *      *    *    +    +    +    +       *    *   u13  u24  u35  u46
 226:              // *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
 227:              // *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
 228:              // *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
 229:              // *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
 230:              // *
 231:              // *  Array elements marked * are not used by the routine; elements marked
 232:              // *  + need not be set on entry, but are required by the routine to store
 233:              // *  elements of U, because of fill-in resulting from the row
 234:              // *  interchanges.
 235:              // *
 236:              // *  =====================================================================
 237:              // *
 238:              // *     .. Parameters ..
 239:              // *     ..
 240:              // *     .. Local Scalars ..
 241:              // *     ..
 242:              // *     .. External Functions ..
 243:              // *     ..
 244:              // *     .. External Subroutines ..
 245:              // *     ..
 246:              // *     .. Intrinsic Functions ..
 247:              //      INTRINSIC          MAX, MIN;
 248:              // *     ..
 249:              // *     .. Executable Statements ..
 250:              // *
 251:              // *     KV is the number of superdiagonals in the factor U, allowing for
 252:              // *     fill-in.
 253:              // *
 254:   
 255:              #endregion
 256:   
 257:   
 258:              #region Body
 259:              
 260:              KV = KU + KL;
 261:              // *
 262:              // *     Test the input parameters.
 263:              // *
 264:              INFO = 0;
 265:              if (M < 0)
 266:              {
 267:                  INFO =  - 1;
 268:              }
 269:              else
 270:              {
 271:                  if (N < 0)
 272:                  {
 273:                      INFO =  - 2;
 274:                  }
 275:                  else
 276:                  {
 277:                      if (KL < 0)
 278:                      {
 279:                          INFO =  - 3;
 280:                      }
 281:                      else
 282:                      {
 283:                          if (KU < 0)
 284:                          {
 285:                              INFO =  - 4;
 286:                          }
 287:                          else
 288:                          {
 289:                              if (LDAB < KL + KV + 1)
 290:                              {
 291:                                  INFO =  - 6;
 292:                              }
 293:                          }
 294:                      }
 295:                  }
 296:              }
 297:              if (INFO != 0)
 298:              {
 299:                  this._xerbla.Run("DGBTF2",  - INFO);
 300:                  return;
 301:              }
 302:              // *
 303:              // *     Quick return if possible
 304:              // *
 305:              if (M == 0 || N == 0) return;
 306:              // *
 307:              // *     Gaussian elimination with partial pivoting
 308:              // *
 309:              // *     Set fill-in elements in columns KU+2 to KV to zero.
 310:              // *
 311:              for (J = KU + 2; J <= Math.Min(KV, N); J++)
 312:              {
 313:                  for (I = KV - J + 2; I <= KL; I++)
 314:                  {
 315:                      AB[I+J * LDAB + o_ab] = ZERO;
 316:                  }
 317:              }
 318:              // *
 319:              // *     JU is the index of the last column affected by the current stage
 320:              // *     of the factorization.
 321:              // *
 322:              JU = 1;
 323:              // *
 324:              for (J = 1; J <= Math.Min(M, N); J++)
 325:              {
 326:                  // *
 327:                  // *        Set fill-in elements in column J+KV to zero.
 328:                  // *
 329:                  if (J + KV <= N)
 330:                  {
 331:                      for (I = 1; I <= KL; I++)
 332:                      {
 333:                          AB[I+(J + KV) * LDAB + o_ab] = ZERO;
 334:                      }
 335:                  }
 336:                  // *
 337:                  // *        Find pivot and test for singularity. KM is the number of
 338:                  // *        subdiagonal elements in the current column.
 339:                  // *
 340:                  KM = Math.Min(KL, M - J);
 341:                  JP = this._idamax.Run(KM + 1, AB, KV + 1+J * LDAB + o_ab, 1);
 342:                  IPIV[J + o_ipiv] = JP + J - 1;
 343:                  if (AB[KV + JP+J * LDAB + o_ab] != ZERO)
 344:                  {
 345:                      JU = Math.Max(JU, Math.Min(J + KU + JP - 1, N));
 346:                      // *
 347:                      // *           Apply interchange to columns J to JU.
 348:                      // *
 349:                      if (JP != 1) this._dswap.Run(JU - J + 1, ref AB, KV + JP+J * LDAB + o_ab, LDAB - 1, ref AB, KV + 1+J * LDAB + o_ab, LDAB - 1);
 350:                      // *
 351:                      if (KM > 0)
 352:                      {
 353:                          // *
 354:                          // *              Compute multipliers.
 355:                          // *
 356:                          this._dscal.Run(KM, ONE / AB[KV + 1+J * LDAB + o_ab], ref AB, KV + 2+J * LDAB + o_ab, 1);
 357:                          // *
 358:                          // *              Update trailing submatrix within the band.
 359:                          // *
 360:                          if (JU > J)
 361:                          {
 362:                              this._dger.Run(KM, JU - J,  - ONE, AB, KV + 2+J * LDAB + o_ab, 1, AB, KV+(J + 1) * LDAB + o_ab
 363:                                             , LDAB - 1, ref AB, KV + 1+(J + 1) * LDAB + o_ab, LDAB - 1);
 364:                          }
 365:                      }
 366:                  }
 367:                  else
 368:                  {
 369:                      // *
 370:                      // *           If pivot is zero, set INFO to the index of the pivot
 371:                      // *           unless a zero pivot has already been found.
 372:                      // *
 373:                      if (INFO == 0) INFO = J;
 374:                  }
 375:              }
 376:              return;
 377:              // *
 378:              // *     End of DGBTF2
 379:              // *
 380:   
 381:              #endregion
 382:   
 383:          }
 384:      }
 385:  }