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:      /// DGBSV computes the solution to a real system of linear equations
  27:      /// A * X = B, where A is a band matrix of order N with KL subdiagonals
  28:      /// and KU superdiagonals, and X and B are N-by-NRHS matrices.
  29:      /// 
  30:      /// The LU decomposition with partial pivoting and row interchanges is
  31:      /// used to factor A as A = L * U, where L is a product of permutation
  32:      /// and unit lower triangular matrices with KL subdiagonals, and U is
  33:      /// upper triangular with KL+KU superdiagonals.  The factored form of A
  34:      /// is then used to solve the system of equations A * X = B.
  35:      /// 
  36:      ///</summary>
  37:      public class DGBSV
  38:      {
  39:      
  40:   
  41:          #region Dependencies
  42:          
  43:          DGBTRF _dgbtrf; DGBTRS _dgbtrs; XERBLA _xerbla; 
  44:   
  45:          #endregion
  46:   
  47:          public DGBSV(DGBTRF dgbtrf, DGBTRS dgbtrs, XERBLA xerbla)
  48:          {
  49:      
  50:   
  51:              #region Set Dependencies
  52:              
  53:              this._dgbtrf = dgbtrf; this._dgbtrs = dgbtrs; this._xerbla = xerbla; 
  54:   
  55:              #endregion
  56:   
  57:          }
  58:      
  59:          public DGBSV()
  60:          {
  61:      
  62:   
  63:              #region Dependencies (Initialization)
  64:              
  65:              IDAMAX idamax = new IDAMAX();
  66:              IEEECK ieeeck = new IEEECK();
  67:              IPARMQ iparmq = new IPARMQ();
  68:              DCOPY dcopy = new DCOPY();
  69:              XERBLA xerbla = new XERBLA();
  70:              DSCAL dscal = new DSCAL();
  71:              DSWAP dswap = new DSWAP();
  72:              LSAME lsame = new LSAME();
  73:              DLASWP dlaswp = new DLASWP();
  74:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  75:              DGER dger = new DGER(xerbla);
  76:              DGBTF2 dgbtf2 = new DGBTF2(idamax, dger, dscal, dswap, xerbla);
  77:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  78:              DTRSM dtrsm = new DTRSM(lsame, xerbla);
  79:              DGBTRF dgbtrf = new DGBTRF(idamax, ilaenv, dcopy, dgbtf2, dgemm, dger, dlaswp, dscal, dswap, dtrsm
  80:                                         , xerbla);
  81:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  82:              DTBSV dtbsv = new DTBSV(lsame, xerbla);
  83:              DGBTRS dgbtrs = new DGBTRS(lsame, dgemv, dger, dswap, dtbsv, xerbla);
  84:   
  85:              #endregion
  86:   
  87:   
  88:              #region Set Dependencies
  89:              
  90:              this._dgbtrf = dgbtrf; this._dgbtrs = dgbtrs; this._xerbla = xerbla; 
  91:   
  92:              #endregion
  93:   
  94:          }
  95:          /// <summary>
  96:          /// Purpose
  97:          /// =======
  98:          /// 
  99:          /// DGBSV computes the solution to a real system of linear equations
 100:          /// A * X = B, where A is a band matrix of order N with KL subdiagonals
 101:          /// and KU superdiagonals, and X and B are N-by-NRHS matrices.
 102:          /// 
 103:          /// The LU decomposition with partial pivoting and row interchanges is
 104:          /// used to factor A as A = L * U, where L is a product of permutation
 105:          /// and unit lower triangular matrices with KL subdiagonals, and U is
 106:          /// upper triangular with KL+KU superdiagonals.  The factored form of A
 107:          /// is then used to solve the system of equations A * X = B.
 108:          /// 
 109:          ///</summary>
 110:          /// <param name="N">
 111:          /// (input) INTEGER
 112:          /// The number of linear equations, i.e., the order of the
 113:          /// matrix A.  N .GE. 0.
 114:          ///</param>
 115:          /// <param name="KL">
 116:          /// (input) INTEGER
 117:          /// The number of subdiagonals within the band of A.  KL .GE. 0.
 118:          ///</param>
 119:          /// <param name="KU">
 120:          /// (input) INTEGER
 121:          /// The number of superdiagonals within the band of A.  KU .GE. 0.
 122:          ///</param>
 123:          /// <param name="NRHS">
 124:          /// (input) INTEGER
 125:          /// The number of right hand sides, i.e., the number of columns
 126:          /// of the matrix B.  NRHS .GE. 0.
 127:          ///</param>
 128:          /// <param name="AB">
 129:          /// (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 130:          /// On entry, the matrix A in band storage, in rows KL+1 to
 131:          /// 2*KL+KU+1; rows 1 to KL of the array need not be set.
 132:          /// The j-th column of A is stored in the j-th column of the
 133:          /// array AB as follows:
 134:          /// AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU).LE.i.LE.min(N,j+KL)
 135:          /// On exit, details of the factorization: U is stored as an
 136:          /// upper triangular band matrix with KL+KU superdiagonals in
 137:          /// rows 1 to KL+KU+1, and the multipliers used during the
 138:          /// factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
 139:          /// See below for further details.
 140:          ///</param>
 141:          /// <param name="LDAB">
 142:          /// (input) INTEGER
 143:          /// The leading dimension of the array AB.  LDAB .GE. 2*KL+KU+1.
 144:          ///</param>
 145:          /// <param name="IPIV">
 146:          /// (output) INTEGER array, dimension (N)
 147:          /// The pivot indices that define the permutation matrix P;
 148:          /// row i of the matrix was interchanged with row IPIV(i).
 149:          ///</param>
 150:          /// <param name="B">
 151:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 152:          /// On entry, the N-by-NRHS right hand side matrix B.
 153:          /// On exit, if INFO = 0, the N-by-NRHS solution matrix X.
 154:          ///</param>
 155:          /// <param name="LDB">
 156:          /// (input) INTEGER
 157:          /// The leading dimension of the array B.  LDB .GE. max(1,N).
 158:          ///</param>
 159:          /// <param name="INFO">
 160:          /// (output) INTEGER
 161:          /// = 0:  successful exit
 162:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 163:          /// .GT. 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
 164:          /// has been completed, but the factor U is exactly
 165:          /// singular, and the solution has not been computed.
 166:          ///</param>
 167:          public void Run(int N, int KL, int KU, int NRHS, ref double[] AB, int offset_ab, int LDAB
 168:                           , ref int[] IPIV, int offset_ipiv, ref double[] B, int offset_b, int LDB, ref int INFO)
 169:          {
 170:   
 171:              #region Array Index Correction
 172:              
 173:               int o_ab = -1 - LDAB + offset_ab;  int o_ipiv = -1 + offset_ipiv;  int o_b = -1 - LDB + offset_b; 
 174:   
 175:              #endregion
 176:   
 177:   
 178:              #region Prolog
 179:              
 180:              // *
 181:              // *  -- LAPACK driver routine (version 3.1) --
 182:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 183:              // *     November 2006
 184:              // *
 185:              // *     .. Scalar Arguments ..
 186:              // *     ..
 187:              // *     .. Array Arguments ..
 188:              // *     ..
 189:              // *
 190:              // *  Purpose
 191:              // *  =======
 192:              // *
 193:              // *  DGBSV computes the solution to a real system of linear equations
 194:              // *  A * X = B, where A is a band matrix of order N with KL subdiagonals
 195:              // *  and KU superdiagonals, and X and B are N-by-NRHS matrices.
 196:              // *
 197:              // *  The LU decomposition with partial pivoting and row interchanges is
 198:              // *  used to factor A as A = L * U, where L is a product of permutation
 199:              // *  and unit lower triangular matrices with KL subdiagonals, and U is
 200:              // *  upper triangular with KL+KU superdiagonals.  The factored form of A
 201:              // *  is then used to solve the system of equations A * X = B.
 202:              // *
 203:              // *  Arguments
 204:              // *  =========
 205:              // *
 206:              // *  N       (input) INTEGER
 207:              // *          The number of linear equations, i.e., the order of the
 208:              // *          matrix A.  N >= 0.
 209:              // *
 210:              // *  KL      (input) INTEGER
 211:              // *          The number of subdiagonals within the band of A.  KL >= 0.
 212:              // *
 213:              // *  KU      (input) INTEGER
 214:              // *          The number of superdiagonals within the band of A.  KU >= 0.
 215:              // *
 216:              // *  NRHS    (input) INTEGER
 217:              // *          The number of right hand sides, i.e., the number of columns
 218:              // *          of the matrix B.  NRHS >= 0.
 219:              // *
 220:              // *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 221:              // *          On entry, the matrix A in band storage, in rows KL+1 to
 222:              // *          2*KL+KU+1; rows 1 to KL of the array need not be set.
 223:              // *          The j-th column of A is stored in the j-th column of the
 224:              // *          array AB as follows:
 225:              // *          AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL)
 226:              // *          On exit, details of the factorization: U is stored as an
 227:              // *          upper triangular band matrix with KL+KU superdiagonals in
 228:              // *          rows 1 to KL+KU+1, and the multipliers used during the
 229:              // *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
 230:              // *          See below for further details.
 231:              // *
 232:              // *  LDAB    (input) INTEGER
 233:              // *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
 234:              // *
 235:              // *  IPIV    (output) INTEGER array, dimension (N)
 236:              // *          The pivot indices that define the permutation matrix P;
 237:              // *          row i of the matrix was interchanged with row IPIV(i).
 238:              // *
 239:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
 240:              // *          On entry, the N-by-NRHS right hand side matrix B.
 241:              // *          On exit, if INFO = 0, the N-by-NRHS solution matrix X.
 242:              // *
 243:              // *  LDB     (input) INTEGER
 244:              // *          The leading dimension of the array B.  LDB >= max(1,N).
 245:              // *
 246:              // *  INFO    (output) INTEGER
 247:              // *          = 0:  successful exit
 248:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 249:              // *          > 0:  if INFO = i, U(i,i) is exactly zero.  The factorization
 250:              // *                has been completed, but the factor U is exactly
 251:              // *                singular, and the solution has not been computed.
 252:              // *
 253:              // *  Further Details
 254:              // *  ===============
 255:              // *
 256:              // *  The band storage scheme is illustrated by the following example, when
 257:              // *  M = N = 6, KL = 2, KU = 1:
 258:              // *
 259:              // *  On entry:                       On exit:
 260:              // *
 261:              // *      *    *    *    +    +    +       *    *    *   u14  u25  u36
 262:              // *      *    *    +    +    +    +       *    *   u13  u24  u35  u46
 263:              // *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
 264:              // *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
 265:              // *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
 266:              // *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
 267:              // *
 268:              // *  Array elements marked * are not used by the routine; elements marked
 269:              // *  + need not be set on entry, but are required by the routine to store
 270:              // *  elements of U because of fill-in resulting from the row interchanges.
 271:              // *
 272:              // *  =====================================================================
 273:              // *
 274:              // *     .. External Subroutines ..
 275:              // *     ..
 276:              // *     .. Intrinsic Functions ..
 277:              //      INTRINSIC          MAX;
 278:              // *     ..
 279:              // *     .. Executable Statements ..
 280:              // *
 281:              // *     Test the input parameters.
 282:              // *
 283:   
 284:              #endregion
 285:   
 286:   
 287:              #region Body
 288:              
 289:              INFO = 0;
 290:              if (N < 0)
 291:              {
 292:                  INFO =  - 1;
 293:              }
 294:              else
 295:              {
 296:                  if (KL < 0)
 297:                  {
 298:                      INFO =  - 2;
 299:                  }
 300:                  else
 301:                  {
 302:                      if (KU < 0)
 303:                      {
 304:                          INFO =  - 3;
 305:                      }
 306:                      else
 307:                      {
 308:                          if (NRHS < 0)
 309:                          {
 310:                              INFO =  - 4;
 311:                          }
 312:                          else
 313:                          {
 314:                              if (LDAB < 2 * KL + KU + 1)
 315:                              {
 316:                                  INFO =  - 6;
 317:                              }
 318:                              else
 319:                              {
 320:                                  if (LDB < Math.Max(N, 1))
 321:                                  {
 322:                                      INFO =  - 9;
 323:                                  }
 324:                              }
 325:                          }
 326:                      }
 327:                  }
 328:              }
 329:              if (INFO != 0)
 330:              {
 331:                  this._xerbla.Run("DGBSV ",  - INFO);
 332:                  return;
 333:              }
 334:              // *
 335:              // *     Compute the LU factorization of the band matrix A.
 336:              // *
 337:              this._dgbtrf.Run(N, N, KL, KU, ref AB, offset_ab, LDAB
 338:                               , ref IPIV, offset_ipiv, ref INFO);
 339:              if (INFO == 0)
 340:              {
 341:                  // *
 342:                  // *        Solve the system A*X = B, overwriting B with X.
 343:                  // *
 344:                  this._dgbtrs.Run("No transpose", N, KL, KU, NRHS, AB, offset_ab
 345:                                   , LDAB, IPIV, offset_ipiv, ref B, offset_b, LDB, ref INFO);
 346:              }
 347:              return;
 348:              // *
 349:              // *     End of DGBSV
 350:              // *
 351:   
 352:              #endregion
 353:   
 354:          }
 355:      }
 356:  }