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:      /// DGBTRF 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 blocked version of the algorithm, calling Level 3 BLAS.
  30:      /// 
  31:      ///</summary>
  32:      public class DGBTRF
  33:      {
  34:      
  35:   
  36:          #region Dependencies
  37:          
  38:          IDAMAX _idamax; ILAENV _ilaenv; DCOPY _dcopy; DGBTF2 _dgbtf2; DGEMM _dgemm; DGER _dger; DLASWP _dlaswp; DSCAL _dscal; 
  39:          DSWAP _dswap;DTRSM _dtrsm; XERBLA _xerbla; 
  40:   
  41:          #endregion
  42:   
  43:   
  44:          #region Fields
  45:          
  46:          const double ONE = 1.0E+0; const double ZERO = 0.0E+0; const int NBMAX = 64; const int LDWORK = NBMAX + 1; int I = 0; 
  47:          int I2 = 0;int I3 = 0; int II = 0; int IP = 0; int J = 0; int J2 = 0; int J3 = 0; int JB = 0; int JJ = 0; int JM = 0; 
  48:          int JP = 0;int JU = 0; int K2 = 0; int KM = 0; int KV = 0; int NB = 0; int NW = 0; double TEMP = 0; 
  49:          double[] WORK13 = new double[LDWORK * NBMAX]; int offset_work13 = 0; int o_work13 = -1 - LDWORK;
  50:          double[] WORK31 = new double[LDWORK * NBMAX]; int offset_work31 = 0; int o_work31 = -1 - LDWORK;
  51:   
  52:          #endregion
  53:   
  54:          public DGBTRF(IDAMAX idamax, ILAENV ilaenv, DCOPY dcopy, DGBTF2 dgbtf2, DGEMM dgemm, DGER dger, DLASWP dlaswp, DSCAL dscal, DSWAP dswap, DTRSM dtrsm
  55:                        , XERBLA xerbla)
  56:          {
  57:      
  58:   
  59:              #region Set Dependencies
  60:              
  61:              this._idamax = idamax; this._ilaenv = ilaenv; this._dcopy = dcopy; this._dgbtf2 = dgbtf2; this._dgemm = dgemm; 
  62:              this._dger = dger;this._dlaswp = dlaswp; this._dscal = dscal; this._dswap = dswap; this._dtrsm = dtrsm; 
  63:              this._xerbla = xerbla;
  64:   
  65:              #endregion
  66:   
  67:          }
  68:      
  69:          public DGBTRF()
  70:          {
  71:      
  72:   
  73:              #region Dependencies (Initialization)
  74:              
  75:              IDAMAX idamax = new IDAMAX();
  76:              IEEECK ieeeck = new IEEECK();
  77:              IPARMQ iparmq = new IPARMQ();
  78:              DCOPY dcopy = new DCOPY();
  79:              XERBLA xerbla = new XERBLA();
  80:              DSCAL dscal = new DSCAL();
  81:              DSWAP dswap = new DSWAP();
  82:              LSAME lsame = new LSAME();
  83:              DLASWP dlaswp = new DLASWP();
  84:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  85:              DGER dger = new DGER(xerbla);
  86:              DGBTF2 dgbtf2 = new DGBTF2(idamax, dger, dscal, dswap, xerbla);
  87:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  88:              DTRSM dtrsm = new DTRSM(lsame, xerbla);
  89:   
  90:              #endregion
  91:   
  92:   
  93:              #region Set Dependencies
  94:              
  95:              this._idamax = idamax; this._ilaenv = ilaenv; this._dcopy = dcopy; this._dgbtf2 = dgbtf2; this._dgemm = dgemm; 
  96:              this._dger = dger;this._dlaswp = dlaswp; this._dscal = dscal; this._dswap = dswap; this._dtrsm = dtrsm; 
  97:              this._xerbla = xerbla;
  98:   
  99:              #endregion
 100:   
 101:          }
 102:          /// <summary>
 103:          /// Purpose
 104:          /// =======
 105:          /// 
 106:          /// DGBTRF computes an LU factorization of a real m-by-n band matrix A
 107:          /// using partial pivoting with row interchanges.
 108:          /// 
 109:          /// This is the blocked version of the algorithm, calling Level 3 BLAS.
 110:          /// 
 111:          ///</summary>
 112:          /// <param name="M">
 113:          /// (input) INTEGER
 114:          /// The number of rows of the matrix A.  M .GE. 0.
 115:          ///</param>
 116:          /// <param name="N">
 117:          /// (input) INTEGER
 118:          /// The number of columns of the matrix A.  N .GE. 0.
 119:          ///</param>
 120:          /// <param name="KL">
 121:          /// (input) INTEGER
 122:          /// The number of subdiagonals within the band of A.  KL .GE. 0.
 123:          ///</param>
 124:          /// <param name="KU">
 125:          /// (input) INTEGER
 126:          /// The number of superdiagonals within the band of A.  KU .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(m,j+kl)
 135:          /// 
 136:          /// On exit, details of the factorization: U is stored as an
 137:          /// upper triangular band matrix with KL+KU superdiagonals in
 138:          /// rows 1 to KL+KU+1, and the multipliers used during the
 139:          /// factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
 140:          /// See below for further details.
 141:          ///</param>
 142:          /// <param name="LDAB">
 143:          /// (input) INTEGER
 144:          /// The leading dimension of the array AB.  LDAB .GE. 2*KL+KU+1.
 145:          ///</param>
 146:          /// <param name="IPIV">
 147:          /// (output) INTEGER array, dimension (min(M,N))
 148:          /// The pivot indices; for 1 .LE. i .LE. min(M,N), row i of the
 149:          /// matrix was interchanged with row IPIV(i).
 150:          ///</param>
 151:          /// <param name="INFO">
 152:          /// (output) INTEGER
 153:          /// = 0: successful exit
 154:          /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
 155:          /// .GT. 0: if INFO = +i, U(i,i) is exactly zero. The factorization
 156:          /// has been completed, but the factor U is exactly
 157:          /// singular, and division by zero will occur if it is used
 158:          /// to solve a system of equations.
 159:          ///</param>
 160:          public void Run(int M, int N, int KL, int KU, ref double[] AB, int offset_ab, int LDAB
 161:                           , ref int[] IPIV, int offset_ipiv, ref int INFO)
 162:          {
 163:   
 164:              #region Array Index Correction
 165:              
 166:               int o_ab = -1 - LDAB + offset_ab;  int o_ipiv = -1 + offset_ipiv; 
 167:   
 168:              #endregion
 169:   
 170:   
 171:              #region Prolog
 172:              
 173:              // *
 174:              // *  -- LAPACK routine (version 3.1) --
 175:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 176:              // *     November 2006
 177:              // *
 178:              // *     .. Scalar Arguments ..
 179:              // *     ..
 180:              // *     .. Array Arguments ..
 181:              // *     ..
 182:              // *
 183:              // *  Purpose
 184:              // *  =======
 185:              // *
 186:              // *  DGBTRF computes an LU factorization of a real m-by-n band matrix A
 187:              // *  using partial pivoting with row interchanges.
 188:              // *
 189:              // *  This is the blocked version of the algorithm, calling Level 3 BLAS.
 190:              // *
 191:              // *  Arguments
 192:              // *  =========
 193:              // *
 194:              // *  M       (input) INTEGER
 195:              // *          The number of rows of the matrix A.  M >= 0.
 196:              // *
 197:              // *  N       (input) INTEGER
 198:              // *          The number of columns of the matrix A.  N >= 0.
 199:              // *
 200:              // *  KL      (input) INTEGER
 201:              // *          The number of subdiagonals within the band of A.  KL >= 0.
 202:              // *
 203:              // *  KU      (input) INTEGER
 204:              // *          The number of superdiagonals within the band of A.  KU >= 0.
 205:              // *
 206:              // *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 207:              // *          On entry, the matrix A in band storage, in rows KL+1 to
 208:              // *          2*KL+KU+1; rows 1 to KL of the array need not be set.
 209:              // *          The j-th column of A is stored in the j-th column of the
 210:              // *          array AB as follows:
 211:              // *          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
 212:              // *
 213:              // *          On exit, details of the factorization: U is stored as an
 214:              // *          upper triangular band matrix with KL+KU superdiagonals in
 215:              // *          rows 1 to KL+KU+1, and the multipliers used during the
 216:              // *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
 217:              // *          See below for further details.
 218:              // *
 219:              // *  LDAB    (input) INTEGER
 220:              // *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
 221:              // *
 222:              // *  IPIV    (output) INTEGER array, dimension (min(M,N))
 223:              // *          The pivot indices; for 1 <= i <= min(M,N), row i of the
 224:              // *          matrix was interchanged with row IPIV(i).
 225:              // *
 226:              // *  INFO    (output) INTEGER
 227:              // *          = 0: successful exit
 228:              // *          < 0: if INFO = -i, the i-th argument had an illegal value
 229:              // *          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
 230:              // *               has been completed, but the factor U is exactly
 231:              // *               singular, and division by zero will occur if it is used
 232:              // *               to solve a system of equations.
 233:              // *
 234:              // *  Further Details
 235:              // *  ===============
 236:              // *
 237:              // *  The band storage scheme is illustrated by the following example, when
 238:              // *  M = N = 6, KL = 2, KU = 1:
 239:              // *
 240:              // *  On entry:                       On exit:
 241:              // *
 242:              // *      *    *    *    +    +    +       *    *    *   u14  u25  u36
 243:              // *      *    *    +    +    +    +       *    *   u13  u24  u35  u46
 244:              // *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
 245:              // *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
 246:              // *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
 247:              // *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
 248:              // *
 249:              // *  Array elements marked * are not used by the routine; elements marked
 250:              // *  + need not be set on entry, but are required by the routine to store
 251:              // *  elements of U because of fill-in resulting from the row interchanges.
 252:              // *
 253:              // *  =====================================================================
 254:              // *
 255:              // *     .. Parameters ..
 256:              // *     ..
 257:              // *     .. Local Scalars ..
 258:              // *     ..
 259:              // *     .. Local Arrays ..
 260:              // *     ..
 261:              // *     .. External Functions ..
 262:              // *     ..
 263:              // *     .. External Subroutines ..
 264:              // *     ..
 265:              // *     .. Intrinsic Functions ..
 266:              //      INTRINSIC          MAX, MIN;
 267:              // *     ..
 268:              // *     .. Executable Statements ..
 269:              // *
 270:              // *     KV is the number of superdiagonals in the factor U, allowing for
 271:              // *     fill-in
 272:              // *
 273:   
 274:              #endregion
 275:   
 276:   
 277:              #region Body
 278:              
 279:              KV = KU + KL;
 280:              // *
 281:              // *     Test the input parameters.
 282:              // *
 283:              INFO = 0;
 284:              if (M < 0)
 285:              {
 286:                  INFO =  - 1;
 287:              }
 288:              else
 289:              {
 290:                  if (N < 0)
 291:                  {
 292:                      INFO =  - 2;
 293:                  }
 294:                  else
 295:                  {
 296:                      if (KL < 0)
 297:                      {
 298:                          INFO =  - 3;
 299:                      }
 300:                      else
 301:                      {
 302:                          if (KU < 0)
 303:                          {
 304:                              INFO =  - 4;
 305:                          }
 306:                          else
 307:                          {
 308:                              if (LDAB < KL + KV + 1)
 309:                              {
 310:                                  INFO =  - 6;
 311:                              }
 312:                          }
 313:                      }
 314:                  }
 315:              }
 316:              if (INFO != 0)
 317:              {
 318:                  this._xerbla.Run("DGBTRF",  - INFO);
 319:                  return;
 320:              }
 321:              // *
 322:              // *     Quick return if possible
 323:              // *
 324:              if (M == 0 || N == 0) return;
 325:              // *
 326:              // *     Determine the block size for this environment
 327:              // *
 328:              NB = this._ilaenv.Run(1, "DGBTRF", " ", M, N, KL, KU);
 329:              // *
 330:              // *     The block size must not exceed the limit set by the size of the
 331:              // *     local arrays WORK13 and WORK31.
 332:              // *
 333:              NB = Math.Min(NB, NBMAX);
 334:              // *
 335:              if (NB <= 1 || NB > KL)
 336:              {
 337:                  // *
 338:                  // *        Use unblocked code
 339:                  // *
 340:                  this._dgbtf2.Run(M, N, KL, KU, ref AB, offset_ab, LDAB
 341:                                   , ref IPIV, offset_ipiv, ref INFO);
 342:              }
 343:              else
 344:              {
 345:                  // *
 346:                  // *        Use blocked code
 347:                  // *
 348:                  // *        Zero the superdiagonal elements of the work array WORK13
 349:                  // *
 350:                  for (J = 1; J <= NB; J++)
 351:                  {
 352:                      for (I = 1; I <= J - 1; I++)
 353:                      {
 354:                          WORK13[I+J * LDWORK + o_work13] = ZERO;
 355:                      }
 356:                  }
 357:                  // *
 358:                  // *        Zero the subdiagonal elements of the work array WORK31
 359:                  // *
 360:                  for (J = 1; J <= NB; J++)
 361:                  {
 362:                      for (I = J + 1; I <= NB; I++)
 363:                      {
 364:                          WORK31[I+J * LDWORK + o_work31] = ZERO;
 365:                      }
 366:                  }
 367:                  // *
 368:                  // *        Gaussian elimination with partial pivoting
 369:                  // *
 370:                  // *        Set fill-in elements in columns KU+2 to KV to zero
 371:                  // *
 372:                  for (J = KU + 2; J <= Math.Min(KV, N); J++)
 373:                  {
 374:                      for (I = KV - J + 2; I <= KL; I++)
 375:                      {
 376:                          AB[I+J * LDAB + o_ab] = ZERO;
 377:                      }
 378:                  }
 379:                  // *
 380:                  // *        JU is the index of the last column affected by the current
 381:                  // *        stage of the factorization
 382:                  // *
 383:                  JU = 1;
 384:                  // *
 385:                  for (J = 1; (NB >= 0) ? (J <= Math.Min(M, N)) : (J >= Math.Min(M, N)); J += NB)
 386:                  {
 387:                      JB = Math.Min(NB, Math.Min(M, N) - J + 1);
 388:                      // *
 389:                      // *           The active part of the matrix is partitioned
 390:                      // *
 391:                      // *              A11   A12   A13
 392:                      // *              A21   A22   A23
 393:                      // *              A31   A32   A33
 394:                      // *
 395:                      // *           Here A11, A21 and A31 denote the current block of JB columns
 396:                      // *           which is about to be factorized. The number of rows in the
 397:                      // *           partitioning are JB, I2, I3 respectively, and the numbers
 398:                      // *           of columns are JB, J2, J3. The superdiagonal elements of A13
 399:                      // *           and the subdiagonal elements of A31 lie outside the band.
 400:                      // *
 401:                      I2 = Math.Min(KL - JB, M - J - JB + 1);
 402:                      I3 = Math.Min(JB, M - J - KL + 1);
 403:                      // *
 404:                      // *           J2 and J3 are computed after JU has been updated.
 405:                      // *
 406:                      // *           Factorize the current block of JB columns
 407:                      // *
 408:                      for (JJ = J; JJ <= J + JB - 1; JJ++)
 409:                      {
 410:                          // *
 411:                          // *              Set fill-in elements in column JJ+KV to zero
 412:                          // *
 413:                          if (JJ + KV <= N)
 414:                          {
 415:                              for (I = 1; I <= KL; I++)
 416:                              {
 417:                                  AB[I+(JJ + KV) * LDAB + o_ab] = ZERO;
 418:                              }
 419:                          }
 420:                          // *
 421:                          // *              Find pivot and test for singularity. KM is the number of
 422:                          // *              subdiagonal elements in the current column.
 423:                          // *
 424:                          KM = Math.Min(KL, M - JJ);
 425:                          JP = this._idamax.Run(KM + 1, AB, KV + 1+JJ * LDAB + o_ab, 1);
 426:                          IPIV[JJ + o_ipiv] = JP + JJ - J;
 427:                          if (AB[KV + JP+JJ * LDAB + o_ab] != ZERO)
 428:                          {
 429:                              JU = Math.Max(JU, Math.Min(JJ + KU + JP - 1, N));
 430:                              if (JP != 1)
 431:                              {
 432:                                  // *
 433:                                  // *                    Apply interchange to columns J to J+JB-1
 434:                                  // *
 435:                                  if (JP + JJ - 1 < J + KL)
 436:                                  {
 437:                                      // *
 438:                                      this._dswap.Run(JB, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref AB, KV + JP + JJ - J+J * LDAB + o_ab, LDAB - 1);
 439:                                  }
 440:                                  else
 441:                                  {
 442:                                      // *
 443:                                      // *                       The interchange affects columns J to JJ-1 of A31
 444:                                      // *                       which are stored in the work array WORK31
 445:                                      // *
 446:                                      this._dswap.Run(JJ - J, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref WORK31, JP + JJ - J - KL+1 * LDWORK + o_work31, LDWORK);
 447:                                      this._dswap.Run(J + JB - JJ, ref AB, KV + 1+JJ * LDAB + o_ab, LDAB - 1, ref AB, KV + JP+JJ * LDAB + o_ab, LDAB - 1);
 448:                                  }
 449:                              }
 450:                              // *
 451:                              // *                 Compute multipliers
 452:                              // *
 453:                              this._dscal.Run(KM, ONE / AB[KV + 1+JJ * LDAB + o_ab], ref AB, KV + 2+JJ * LDAB + o_ab, 1);
 454:                              // *
 455:                              // *                 Update trailing submatrix within the band and within
 456:                              // *                 the current block. JM is the index of the last column
 457:                              // *                 which needs to be updated.
 458:                              // *
 459:                              JM = Math.Min(JU, J + JB - 1);
 460:                              if (JM > JJ)
 461:                              {
 462:                                  this._dger.Run(KM, JM - JJ,  - ONE, AB, KV + 2+JJ * LDAB + o_ab, 1, AB, KV+(JJ + 1) * LDAB + o_ab
 463:                                                 , LDAB - 1, ref AB, KV + 1+(JJ + 1) * LDAB + o_ab, LDAB - 1);
 464:                              }
 465:                          }
 466:                          else
 467:                          {
 468:                              // *
 469:                              // *                 If pivot is zero, set INFO to the index of the pivot
 470:                              // *                 unless a zero pivot has already been found.
 471:                              // *
 472:                              if (INFO == 0) INFO = JJ;
 473:                          }
 474:                          // *
 475:                          // *              Copy current column of A31 into the work array WORK31
 476:                          // *
 477:                          NW = Math.Min(JJ - J + 1, I3);
 478:                          if (NW > 0) this._dcopy.Run(NW, AB, KV + KL + 1 - JJ + J+JJ * LDAB + o_ab, 1, ref WORK31, 1+(JJ - J + 1) * LDWORK + o_work31, 1);
 479:                      }
 480:                      if (J + JB <= N)
 481:                      {
 482:                          // *
 483:                          // *              Apply the row interchanges to the other blocks.
 484:                          // *
 485:                          J2 = Math.Min(JU - J + 1, KV) - JB;
 486:                          J3 = Math.Max(0, JU - J - KV + 1);
 487:                          // *
 488:                          // *              Use DLASWP to apply the row interchanges to A12, A22, and
 489:                          // *              A32.
 490:                          // *
 491:                          this._dlaswp.Run(J2, ref AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1, 1, JB, IPIV, J + o_ipiv
 492:                                           , 1);
 493:                          // *
 494:                          // *              Adjust the pivot indices.
 495:                          // *
 496:                          for (I = J; I <= J + JB - 1; I++)
 497:                          {
 498:                              IPIV[I + o_ipiv] = IPIV[I + o_ipiv] + J - 1;
 499:                          }
 500:                          // *
 501:                          // *              Apply the row interchanges to A13, A23, and A33
 502:                          // *              columnwise.
 503:                          // *
 504:                          K2 = J - 1 + JB + J2;
 505:                          for (I = 1; I <= J3; I++)
 506:                          {
 507:                              JJ = K2 + I;
 508:                              for (II = J + I - 1; II <= J + JB - 1; II++)
 509:                              {
 510:                                  IP = IPIV[II + o_ipiv];
 511:                                  if (IP != II)
 512:                                  {
 513:                                      TEMP = AB[KV + 1 + II - JJ+JJ * LDAB + o_ab];
 514:                                      AB[KV + 1 + II - JJ+JJ * LDAB + o_ab] = AB[KV + 1 + IP - JJ+JJ * LDAB + o_ab];
 515:                                      AB[KV + 1 + IP - JJ+JJ * LDAB + o_ab] = TEMP;
 516:                                  }
 517:                              }
 518:                          }
 519:                          // *
 520:                          // *              Update the relevant part of the trailing submatrix
 521:                          // *
 522:                          if (J2 > 0)
 523:                          {
 524:                              // *
 525:                              // *                 Update A12
 526:                              // *
 527:                              this._dtrsm.Run("Left", "Lower", "No transpose", "Unit", JB, J2
 528:                                              , ONE, AB, KV + 1+J * LDAB + o_ab, LDAB - 1, ref AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1);
 529:                              // *
 530:                              if (I2 > 0)
 531:                              {
 532:                                  // *
 533:                                  // *                    Update A22
 534:                                  // *
 535:                                  this._dgemm.Run("No transpose", "No transpose", I2, J2, JB,  - ONE
 536:                                                  , AB, KV + 1 + JB+J * LDAB + o_ab, LDAB - 1, AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1, ONE, ref AB, KV + 1+(J + JB) * LDAB + o_ab
 537:                                                  , LDAB - 1);
 538:                              }
 539:                              // *
 540:                              if (I3 > 0)
 541:                              {
 542:                                  // *
 543:                                  // *                    Update A32
 544:                                  // *
 545:                                  this._dgemm.Run("No transpose", "No transpose", I3, J2, JB,  - ONE
 546:                                                  , WORK31, offset_work31, LDWORK, AB, KV + 1 - JB+(J + JB) * LDAB + o_ab, LDAB - 1, ONE, ref AB, KV + KL + 1 - JB+(J + JB) * LDAB + o_ab
 547:                                                  , LDAB - 1);
 548:                              }
 549:                          }
 550:                          // *
 551:                          if (J3 > 0)
 552:                          {
 553:                              // *
 554:                              // *                 Copy the lower triangle of A13 into the work array
 555:                              // *                 WORK13
 556:                              // *
 557:                              for (JJ = 1; JJ <= J3; JJ++)
 558:                              {
 559:                                  for (II = JJ; II <= JB; II++)
 560:                                  {
 561:                                      WORK13[II+JJ * LDWORK + o_work13] = AB[II - JJ + 1+(JJ + J + KV - 1) * LDAB + o_ab];
 562:                                  }
 563:                              }
 564:                              // *
 565:                              // *                 Update A13 in the work array
 566:                              // *
 567:                              this._dtrsm.Run("Left", "Lower", "No transpose", "Unit", JB, J3
 568:                                              , ONE, AB, KV + 1+J * LDAB + o_ab, LDAB - 1, ref WORK13, offset_work13, LDWORK);
 569:                              // *
 570:                              if (I2 > 0)
 571:                              {
 572:                                  // *
 573:                                  // *                    Update A23
 574:                                  // *
 575:                                  this._dgemm.Run("No transpose", "No transpose", I2, J3, JB,  - ONE
 576:                                                  , AB, KV + 1 + JB+J * LDAB + o_ab, LDAB - 1, WORK13, offset_work13, LDWORK, ONE, ref AB, 1 + JB+(J + KV) * LDAB + o_ab
 577:                                                  , LDAB - 1);
 578:                              }
 579:                              // *
 580:                              if (I3 > 0)
 581:                              {
 582:                                  // *
 583:                                  // *                    Update A33
 584:                                  // *
 585:                                  this._dgemm.Run("No transpose", "No transpose", I3, J3, JB,  - ONE
 586:                                                  , WORK31, offset_work31, LDWORK, WORK13, offset_work13, LDWORK, ONE, ref AB, 1 + KL+(J + KV) * LDAB + o_ab
 587:                                                  , LDAB - 1);
 588:                              }
 589:                              // *
 590:                              // *                 Copy the lower triangle of A13 back into place
 591:                              // *
 592:                              for (JJ = 1; JJ <= J3; JJ++)
 593:                              {
 594:                                  for (II = JJ; II <= JB; II++)
 595:                                  {
 596:                                      AB[II - JJ + 1+(JJ + J + KV - 1) * LDAB + o_ab] = WORK13[II+JJ * LDWORK + o_work13];
 597:                                  }
 598:                              }
 599:                          }
 600:                      }
 601:                      else
 602:                      {
 603:                          // *
 604:                          // *              Adjust the pivot indices.
 605:                          // *
 606:                          for (I = J; I <= J + JB - 1; I++)
 607:                          {
 608:                              IPIV[I + o_ipiv] = IPIV[I + o_ipiv] + J - 1;
 609:                          }
 610:                      }
 611:                      // *
 612:                      // *           Partially undo the interchanges in the current block to
 613:                      // *           restore the upper triangular form of A31 and copy the upper
 614:                      // *           triangle of A31 back into place
 615:                      // *
 616:                      for (JJ = J + JB - 1; JJ >= J; JJ +=  - 1)
 617:                      {
 618:                          JP = IPIV[JJ + o_ipiv] - JJ + 1;
 619:                          if (JP != 1)
 620:                          {
 621:                              // *
 622:                              // *                 Apply interchange to columns J to JJ-1
 623:                              // *
 624:                              if (JP + JJ - 1 < J + KL)
 625:                              {
 626:                                  // *
 627:                                  // *                    The interchange does not affect A31
 628:                                  // *
 629:                                  this._dswap.Run(JJ - J, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref AB, KV + JP + JJ - J+J * LDAB + o_ab, LDAB - 1);
 630:                              }
 631:                              else
 632:                              {
 633:                                  // *
 634:                                  // *                    The interchange does affect A31
 635:                                  // *
 636:                                  this._dswap.Run(JJ - J, ref AB, KV + 1 + JJ - J+J * LDAB + o_ab, LDAB - 1, ref WORK31, JP + JJ - J - KL+1 * LDWORK + o_work31, LDWORK);
 637:                              }
 638:                          }
 639:                          // *
 640:                          // *              Copy the current column of A31 back into place
 641:                          // *
 642:                          NW = Math.Min(I3, JJ - J + 1);
 643:                          if (NW > 0) this._dcopy.Run(NW, WORK31, 1+(JJ - J + 1) * LDWORK + o_work31, 1, ref AB, KV + KL + 1 - JJ + J+JJ * LDAB + o_ab, 1);
 644:                      }
 645:                  }
 646:              }
 647:              // *
 648:              return;
 649:              // *
 650:              // *     End of DGBTRF
 651:              // *
 652:   
 653:              #endregion
 654:   
 655:          }
 656:      }
 657:  }