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:      /// DSBTRD reduces a real symmetric band matrix A to symmetric
  27:      /// tridiagonal form T by an orthogonal similarity transformation:
  28:      /// Q**T * A * Q = T.
  29:      /// 
  30:      ///</summary>
  31:      public class DSBTRD
  32:      {
  33:      
  34:   
  35:          #region Dependencies
  36:          
  37:          DLAR2V _dlar2v; DLARGV _dlargv; DLARTG _dlartg; DLARTV _dlartv; DLASET _dlaset; DROT _drot; XERBLA _xerbla; LSAME _lsame; 
  38:   
  39:          #endregion
  40:   
  41:   
  42:          #region Fields
  43:          
  44:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool INITQ = false; bool UPPER = false; bool WANTQ = false; 
  45:          int I = 0;int I2 = 0; int IBL = 0; int INCA = 0; int INCX = 0; int IQAEND = 0; int IQB = 0; int IQEND = 0; int J = 0; 
  46:          int J1 = 0;int J1END = 0; int J1INC = 0; int J2 = 0; int JEND = 0; int JIN = 0; int JINC = 0; int K = 0; int KD1 = 0; 
  47:          int KDM1 = 0;int KDN = 0; int L = 0; int LAST = 0; int LEND = 0; int NQ = 0; int NR = 0; int NRT = 0; double TEMP = 0; 
  48:   
  49:          #endregion
  50:   
  51:          public DSBTRD(DLAR2V dlar2v, DLARGV dlargv, DLARTG dlartg, DLARTV dlartv, DLASET dlaset, DROT drot, XERBLA xerbla, LSAME lsame)
  52:          {
  53:      
  54:   
  55:              #region Set Dependencies
  56:              
  57:              this._dlar2v = dlar2v; this._dlargv = dlargv; this._dlartg = dlartg; this._dlartv = dlartv; this._dlaset = dlaset; 
  58:              this._drot = drot;this._xerbla = xerbla; this._lsame = lsame; 
  59:   
  60:              #endregion
  61:   
  62:          }
  63:      
  64:          public DSBTRD()
  65:          {
  66:      
  67:   
  68:              #region Dependencies (Initialization)
  69:              
  70:              DLAR2V dlar2v = new DLAR2V();
  71:              DLARGV dlargv = new DLARGV();
  72:              LSAME lsame = new LSAME();
  73:              DLAMC3 dlamc3 = new DLAMC3();
  74:              DLARTV dlartv = new DLARTV();
  75:              DROT drot = new DROT();
  76:              XERBLA xerbla = new XERBLA();
  77:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  78:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  79:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  80:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  81:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  82:              DLARTG dlartg = new DLARTG(dlamch);
  83:              DLASET dlaset = new DLASET(lsame);
  84:   
  85:              #endregion
  86:   
  87:   
  88:              #region Set Dependencies
  89:              
  90:              this._dlar2v = dlar2v; this._dlargv = dlargv; this._dlartg = dlartg; this._dlartv = dlartv; this._dlaset = dlaset; 
  91:              this._drot = drot;this._xerbla = xerbla; this._lsame = lsame; 
  92:   
  93:              #endregion
  94:   
  95:          }
  96:          /// <summary>
  97:          /// Purpose
  98:          /// =======
  99:          /// 
 100:          /// DSBTRD reduces a real symmetric band matrix A to symmetric
 101:          /// tridiagonal form T by an orthogonal similarity transformation:
 102:          /// Q**T * A * Q = T.
 103:          /// 
 104:          ///</summary>
 105:          /// <param name="VECT">
 106:          /// (input) CHARACTER*1
 107:          /// = 'N':  do not form Q;
 108:          /// = 'V':  form Q;
 109:          /// = 'U':  update a matrix X, by forming X*Q.
 110:          ///</param>
 111:          /// <param name="UPLO">
 112:          /// (input) CHARACTER*1
 113:          /// = 'U':  Upper triangle of A is stored;
 114:          /// = 'L':  Lower triangle of A is stored.
 115:          ///</param>
 116:          /// <param name="N">
 117:          /// (input) INTEGER
 118:          /// The order of the matrix A.  N .GE. 0.
 119:          ///</param>
 120:          /// <param name="KD">
 121:          /// (input) INTEGER
 122:          /// The number of superdiagonals of the matrix A if UPLO = 'U',
 123:          /// or the number of subdiagonals if UPLO = 'L'.  KD .GE. 0.
 124:          ///</param>
 125:          /// <param name="AB">
 126:          /// (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 127:          /// On entry, the upper or lower triangle of the symmetric band
 128:          /// matrix A, stored in the first KD+1 rows of the array.  The
 129:          /// j-th column of A is stored in the j-th column of the array AB
 130:          /// as follows:
 131:          /// if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd).LE.i.LE.j;
 132:          /// if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j.LE.i.LE.min(n,j+kd).
 133:          /// On exit, the diagonal elements of AB are overwritten by the
 134:          /// diagonal elements of the tridiagonal matrix T; if KD .GT. 0, the
 135:          /// elements on the first superdiagonal (if UPLO = 'U') or the
 136:          /// first subdiagonal (if UPLO = 'L') are overwritten by the
 137:          /// off-diagonal elements of T; the rest of AB is overwritten by
 138:          /// values generated during the reduction.
 139:          ///</param>
 140:          /// <param name="LDAB">
 141:          /// (input) INTEGER
 142:          /// The leading dimension of the array AB.  LDAB .GE. KD+1.
 143:          ///</param>
 144:          /// <param name="D">
 145:          /// (output) DOUBLE PRECISION array, dimension (N)
 146:          /// The diagonal elements of the tridiagonal matrix T.
 147:          ///</param>
 148:          /// <param name="E">
 149:          /// (output) DOUBLE PRECISION array, dimension (N-1)
 150:          /// The off-diagonal elements of the tridiagonal matrix T:
 151:          /// E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
 152:          ///</param>
 153:          /// <param name="Q">
 154:          /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
 155:          /// On entry, if VECT = 'U', then Q must contain an N-by-N
 156:          /// matrix X; if VECT = 'N' or 'V', then Q need not be set.
 157:          /// 
 158:          /// On exit:
 159:          /// if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
 160:          /// if VECT = 'U', Q contains the product X*Q;
 161:          /// if VECT = 'N', the array Q is not referenced.
 162:          ///</param>
 163:          /// <param name="LDQ">
 164:          /// (input) INTEGER
 165:          /// The leading dimension of the array Q.
 166:          /// LDQ .GE. 1, and LDQ .GE. N if VECT = 'V' or 'U'.
 167:          ///</param>
 168:          /// <param name="WORK">
 169:          /// (workspace) DOUBLE PRECISION array, dimension (N)
 170:          ///</param>
 171:          /// <param name="INFO">
 172:          /// (output) INTEGER
 173:          /// = 0:  successful exit
 174:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 175:          ///</param>
 176:          public void Run(string VECT, string UPLO, int N, int KD, ref double[] AB, int offset_ab, int LDAB
 177:                           , ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] Q, int offset_q, int LDQ, ref double[] WORK, int offset_work, ref int INFO)
 178:          {
 179:   
 180:              #region Array Index Correction
 181:              
 182:               int o_ab = -1 - LDAB + offset_ab;  int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_q = -1 - LDQ + offset_q; 
 183:               int o_work = -1 + offset_work;
 184:   
 185:              #endregion
 186:   
 187:   
 188:              #region Strings
 189:              
 190:              VECT = VECT.Substring(0, 1);  UPLO = UPLO.Substring(0, 1);  
 191:   
 192:              #endregion
 193:   
 194:   
 195:              #region Prolog
 196:              
 197:              // *
 198:              // *  -- LAPACK routine (version 3.1) --
 199:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 200:              // *     November 2006
 201:              // *
 202:              // *     .. Scalar Arguments ..
 203:              // *     ..
 204:              // *     .. Array Arguments ..
 205:              // *     ..
 206:              // *
 207:              // *  Purpose
 208:              // *  =======
 209:              // *
 210:              // *  DSBTRD reduces a real symmetric band matrix A to symmetric
 211:              // *  tridiagonal form T by an orthogonal similarity transformation:
 212:              // *  Q**T * A * Q = T.
 213:              // *
 214:              // *  Arguments
 215:              // *  =========
 216:              // *
 217:              // *  VECT    (input) CHARACTER*1
 218:              // *          = 'N':  do not form Q;
 219:              // *          = 'V':  form Q;
 220:              // *          = 'U':  update a matrix X, by forming X*Q.
 221:              // *
 222:              // *  UPLO    (input) CHARACTER*1
 223:              // *          = 'U':  Upper triangle of A is stored;
 224:              // *          = 'L':  Lower triangle of A is stored.
 225:              // *
 226:              // *  N       (input) INTEGER
 227:              // *          The order of the matrix A.  N >= 0.
 228:              // *
 229:              // *  KD      (input) INTEGER
 230:              // *          The number of superdiagonals of the matrix A if UPLO = 'U',
 231:              // *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
 232:              // *
 233:              // *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
 234:              // *          On entry, the upper or lower triangle of the symmetric band
 235:              // *          matrix A, stored in the first KD+1 rows of the array.  The
 236:              // *          j-th column of A is stored in the j-th column of the array AB
 237:              // *          as follows:
 238:              // *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
 239:              // *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
 240:              // *          On exit, the diagonal elements of AB are overwritten by the
 241:              // *          diagonal elements of the tridiagonal matrix T; if KD > 0, the
 242:              // *          elements on the first superdiagonal (if UPLO = 'U') or the
 243:              // *          first subdiagonal (if UPLO = 'L') are overwritten by the
 244:              // *          off-diagonal elements of T; the rest of AB is overwritten by
 245:              // *          values generated during the reduction.
 246:              // *
 247:              // *  LDAB    (input) INTEGER
 248:              // *          The leading dimension of the array AB.  LDAB >= KD+1.
 249:              // *
 250:              // *  D       (output) DOUBLE PRECISION array, dimension (N)
 251:              // *          The diagonal elements of the tridiagonal matrix T.
 252:              // *
 253:              // *  E       (output) DOUBLE PRECISION array, dimension (N-1)
 254:              // *          The off-diagonal elements of the tridiagonal matrix T:
 255:              // *          E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
 256:              // *
 257:              // *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
 258:              // *          On entry, if VECT = 'U', then Q must contain an N-by-N
 259:              // *          matrix X; if VECT = 'N' or 'V', then Q need not be set.
 260:              // *
 261:              // *          On exit:
 262:              // *          if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
 263:              // *          if VECT = 'U', Q contains the product X*Q;
 264:              // *          if VECT = 'N', the array Q is not referenced.
 265:              // *
 266:              // *  LDQ     (input) INTEGER
 267:              // *          The leading dimension of the array Q.
 268:              // *          LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
 269:              // *
 270:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
 271:              // *
 272:              // *  INFO    (output) INTEGER
 273:              // *          = 0:  successful exit
 274:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 275:              // *
 276:              // *  Further Details
 277:              // *  ===============
 278:              // *
 279:              // *  Modified by Linda Kaufman, Bell Labs.
 280:              // *
 281:              // *  =====================================================================
 282:              // *
 283:              // *     .. Parameters ..
 284:              // *     ..
 285:              // *     .. Local Scalars ..
 286:              // *     ..
 287:              // *     .. External Subroutines ..
 288:              // *     ..
 289:              // *     .. Intrinsic Functions ..
 290:              //      INTRINSIC          MAX, MIN;
 291:              // *     ..
 292:              // *     .. External Functions ..
 293:              // *     ..
 294:              // *     .. Executable Statements ..
 295:              // *
 296:              // *     Test the input parameters
 297:              // *
 298:   
 299:              #endregion
 300:   
 301:   
 302:              #region Body
 303:              
 304:              INITQ = this._lsame.Run(VECT, "V");
 305:              WANTQ = INITQ || this._lsame.Run(VECT, "U");
 306:              UPPER = this._lsame.Run(UPLO, "U");
 307:              KD1 = KD + 1;
 308:              KDM1 = KD - 1;
 309:              INCX = LDAB - 1;
 310:              IQEND = 1;
 311:              // *
 312:              INFO = 0;
 313:              if (!WANTQ && !this._lsame.Run(VECT, "N"))
 314:              {
 315:                  INFO =  - 1;
 316:              }
 317:              else
 318:              {
 319:                  if (!UPPER && !this._lsame.Run(UPLO, "L"))
 320:                  {
 321:                      INFO =  - 2;
 322:                  }
 323:                  else
 324:                  {
 325:                      if (N < 0)
 326:                      {
 327:                          INFO =  - 3;
 328:                      }
 329:                      else
 330:                      {
 331:                          if (KD < 0)
 332:                          {
 333:                              INFO =  - 4;
 334:                          }
 335:                          else
 336:                          {
 337:                              if (LDAB < KD1)
 338:                              {
 339:                                  INFO =  - 6;
 340:                              }
 341:                              else
 342:                              {
 343:                                  if (LDQ < Math.Max(1, N) && WANTQ)
 344:                                  {
 345:                                      INFO =  - 10;
 346:                                  }
 347:                              }
 348:                          }
 349:                      }
 350:                  }
 351:              }
 352:              if (INFO != 0)
 353:              {
 354:                  this._xerbla.Run("DSBTRD",  - INFO);
 355:                  return;
 356:              }
 357:              // *
 358:              // *     Quick return if possible
 359:              // *
 360:              if (N == 0) return;
 361:              // *
 362:              // *     Initialize Q to the unit matrix, if needed
 363:              // *
 364:              if (INITQ)
 365:              {
 366:                  this._dlaset.Run("Full", N, N, ZERO, ONE, ref Q, offset_q
 367:                                   , LDQ);
 368:              }
 369:              // *
 370:              // *     Wherever possible, plane rotations are generated and applied in
 371:              // *     vector operations of length NR over the index set J1:J2:KD1.
 372:              // *
 373:              // *     The cosines and sines of the plane rotations are stored in the
 374:              // *     arrays D and WORK.
 375:              // *
 376:              INCA = KD1 * LDAB;
 377:              KDN = Math.Min(N - 1, KD);
 378:              if (UPPER)
 379:              {
 380:                  // *
 381:                  if (KD > 1)
 382:                  {
 383:                      // *
 384:                      // *           Reduce to tridiagonal form, working with upper triangle
 385:                      // *
 386:                      NR = 0;
 387:                      J1 = KDN + 2;
 388:                      J2 = 1;
 389:                      // *
 390:                      for (I = 1; I <= N - 2; I++)
 391:                      {
 392:                          // *
 393:                          // *              Reduce i-th row of matrix to tridiagonal form
 394:                          // *
 395:                          for (K = KDN + 1; K >= 2; K +=  - 1)
 396:                          {
 397:                              J1 = J1 + KDN;
 398:                              J2 = J2 + KDN;
 399:                              // *
 400:                              if (NR > 0)
 401:                              {
 402:                                  // *
 403:                                  // *                    generate plane rotations to annihilate nonzero
 404:                                  // *                    elements which have been created outside the band
 405:                                  // *
 406:                                  this._dlargv.Run(NR, ref AB, 1+(J1 - 1) * LDAB + o_ab, INCA, ref WORK, J1 + o_work, KD1, ref D, J1 + o_d
 407:                                                   , KD1);
 408:                                  // *
 409:                                  // *                    apply rotations from the right
 410:                                  // *
 411:                                  // *
 412:                                  // *                    Dependent on the the number of diagonals either
 413:                                  // *                    DLARTV or DROT is used
 414:                                  // *
 415:                                  if (NR >= 2 * KD - 1)
 416:                                  {
 417:                                      for (L = 1; L <= KD - 1; L++)
 418:                                      {
 419:                                          this._dlartv.Run(NR, ref AB, L + 1+(J1 - 1) * LDAB + o_ab, INCA, ref AB, L+J1 * LDAB + o_ab, INCA, D, J1 + o_d
 420:                                                           , WORK, J1 + o_work, KD1);
 421:                                      }
 422:                                      // *
 423:                                  }
 424:                                  else
 425:                                  {
 426:                                      JEND = J1 + (NR - 1) * KD1;
 427:                                      for (JINC = J1; (KD1 >= 0) ? (JINC <= JEND) : (JINC >= JEND); JINC += KD1)
 428:                                      {
 429:                                          this._drot.Run(KDM1, ref AB, 2+(JINC - 1) * LDAB + o_ab, 1, ref AB, 1+JINC * LDAB + o_ab, 1, D[JINC + o_d]
 430:                                                         , WORK[JINC + o_work]);
 431:                                      }
 432:                                  }
 433:                              }
 434:                              // *
 435:                              // *
 436:                              if (K > 2)
 437:                              {
 438:                                  if (K <= N - I + 1)
 439:                                  {
 440:                                      // *
 441:                                      // *                       generate plane rotation to annihilate a(i,i+k-1)
 442:                                      // *                       within the band
 443:                                      // *
 444:                                      this._dlartg.Run(AB[KD - K + 3+(I + K - 2) * LDAB + o_ab], AB[KD - K + 2+(I + K - 1) * LDAB + o_ab], ref D[I + K - 1 + o_d], ref WORK[I + K - 1 + o_work], ref TEMP);
 445:                                      AB[KD - K + 3+(I + K - 2) * LDAB + o_ab] = TEMP;
 446:                                      // *
 447:                                      // *                       apply rotation from the right
 448:                                      // *
 449:                                      this._drot.Run(K - 3, ref AB, KD - K + 4+(I + K - 2) * LDAB + o_ab, 1, ref AB, KD - K + 3+(I + K - 1) * LDAB + o_ab, 1, D[I + K - 1 + o_d]
 450:                                                     , WORK[I + K - 1 + o_work]);
 451:                                  }
 452:                                  NR = NR + 1;
 453:                                  J1 = J1 - KDN - 1;
 454:                              }
 455:                              // *
 456:                              // *                 apply plane rotations from both sides to diagonal
 457:                              // *                 blocks
 458:                              // *
 459:                              if (NR > 0)
 460:                              {
 461:                                  this._dlar2v.Run(NR, ref AB, KD1+(J1 - 1) * LDAB + o_ab, ref AB, KD1+J1 * LDAB + o_ab, ref AB, KD+J1 * LDAB + o_ab, INCA, D, J1 + o_d
 462:                                                   , WORK, J1 + o_work, KD1);
 463:                              }
 464:                              // *
 465:                              // *                 apply plane rotations from the left
 466:                              // *
 467:                              if (NR > 0)
 468:                              {
 469:                                  if (2 * KD - 1 < NR)
 470:                                  {
 471:                                      // *
 472:                                      // *                    Dependent on the the number of diagonals either
 473:                                      // *                    DLARTV or DROT is used
 474:                                      // *
 475:                                      for (L = 1; L <= KD - 1; L++)
 476:                                      {
 477:                                          if (J2 + L > N)
 478:                                          {
 479:                                              NRT = NR - 1;
 480:                                          }
 481:                                          else
 482:                                          {
 483:                                              NRT = NR;
 484:                                          }
 485:                                          if (NRT > 0)
 486:                                          {
 487:                                              this._dlartv.Run(NRT, ref AB, KD - L+(J1 + L) * LDAB + o_ab, INCA, ref AB, KD - L + 1+(J1 + L) * LDAB + o_ab, INCA, D, J1 + o_d
 488:                                                               , WORK, J1 + o_work, KD1);
 489:                                          }
 490:                                      }
 491:                                  }
 492:                                  else
 493:                                  {
 494:                                      J1END = J1 + KD1 * (NR - 2);
 495:                                      if (J1END >= J1)
 496:                                      {
 497:                                          for (JIN = J1; (KD1 >= 0) ? (JIN <= J1END) : (JIN >= J1END); JIN += KD1)
 498:                                          {
 499:                                              this._drot.Run(KD - 1, ref AB, KD - 1+(JIN + 1) * LDAB + o_ab, INCX, ref AB, KD+(JIN + 1) * LDAB + o_ab, INCX, D[JIN + o_d]
 500:                                                             , WORK[JIN + o_work]);
 501:                                          }
 502:                                      }
 503:                                      LEND = Math.Min(KDM1, N - J2);
 504:                                      LAST = J1END + KD1;
 505:                                      if (LEND > 0)
 506:                                      {
 507:                                          this._drot.Run(LEND, ref AB, KD - 1+(LAST + 1) * LDAB + o_ab, INCX, ref AB, KD+(LAST + 1) * LDAB + o_ab, INCX, D[LAST + o_d]
 508:                                                         , WORK[LAST + o_work]);
 509:                                      }
 510:                                  }
 511:                              }
 512:                              // *
 513:                              if (WANTQ)
 514:                              {
 515:                                  // *
 516:                                  // *                    accumulate product of plane rotations in Q
 517:                                  // *
 518:                                  if (INITQ)
 519:                                  {
 520:                                      // *
 521:                                      // *                 take advantage of the fact that Q was
 522:                                      // *                 initially the Identity matrix
 523:                                      // *
 524:                                      IQEND = Math.Max(IQEND, J2);
 525:                                      I2 = Math.Max(0, K - 3);
 526:                                      IQAEND = 1 + I * KD;
 527:                                      if (K == 2) IQAEND = IQAEND + KD;
 528:                                      IQAEND = Math.Min(IQAEND, IQEND);
 529:                                      for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
 530:                                      {
 531:                                          IBL = I - I2 / KDM1;
 532:                                          I2 = I2 + 1;
 533:                                          IQB = Math.Max(1, J - IBL);
 534:                                          NQ = 1 + IQAEND - IQB;
 535:                                          IQAEND = Math.Min(IQAEND + KD, IQEND);
 536:                                          this._drot.Run(NQ, ref Q, IQB+(J - 1) * LDQ + o_q, 1, ref Q, IQB+J * LDQ + o_q, 1, D[J + o_d]
 537:                                                         , WORK[J + o_work]);
 538:                                      }
 539:                                  }
 540:                                  else
 541:                                  {
 542:                                      // *
 543:                                      for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
 544:                                      {
 545:                                          this._drot.Run(N, ref Q, 1+(J - 1) * LDQ + o_q, 1, ref Q, 1+J * LDQ + o_q, 1, D[J + o_d]
 546:                                                         , WORK[J + o_work]);
 547:                                      }
 548:                                  }
 549:                                  // *
 550:                              }
 551:                              // *
 552:                              if (J2 + KDN > N)
 553:                              {
 554:                                  // *
 555:                                  // *                    adjust J2 to keep within the bounds of the matrix
 556:                                  // *
 557:                                  NR = NR - 1;
 558:                                  J2 = J2 - KDN - 1;
 559:                              }
 560:                              // *
 561:                              for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
 562:                              {
 563:                                  // *
 564:                                  // *                    create nonzero element a(j-1,j+kd) outside the band
 565:                                  // *                    and store it in WORK
 566:                                  // *
 567:                                  WORK[J + KD + o_work] = WORK[J + o_work] * AB[1+(J + KD) * LDAB + o_ab];
 568:                                  AB[1+(J + KD) * LDAB + o_ab] = D[J + o_d] * AB[1+(J + KD) * LDAB + o_ab];
 569:                              }
 570:                          }
 571:                      }
 572:                  }
 573:                  // *
 574:                  if (KD > 0)
 575:                  {
 576:                      // *
 577:                      // *           copy off-diagonal elements to E
 578:                      // *
 579:                      for (I = 1; I <= N - 1; I++)
 580:                      {
 581:                          E[I + o_e] = AB[KD+(I + 1) * LDAB + o_ab];
 582:                      }
 583:                  }
 584:                  else
 585:                  {
 586:                      // *
 587:                      // *           set E to zero if original matrix was diagonal
 588:                      // *
 589:                      for (I = 1; I <= N - 1; I++)
 590:                      {
 591:                          E[I + o_e] = ZERO;
 592:                      }
 593:                  }
 594:                  // *
 595:                  // *        copy diagonal elements to D
 596:                  // *
 597:                  for (I = 1; I <= N; I++)
 598:                  {
 599:                      D[I + o_d] = AB[KD1+I * LDAB + o_ab];
 600:                  }
 601:                  // *
 602:              }
 603:              else
 604:              {
 605:                  // *
 606:                  if (KD > 1)
 607:                  {
 608:                      // *
 609:                      // *           Reduce to tridiagonal form, working with lower triangle
 610:                      // *
 611:                      NR = 0;
 612:                      J1 = KDN + 2;
 613:                      J2 = 1;
 614:                      // *
 615:                      for (I = 1; I <= N - 2; I++)
 616:                      {
 617:                          // *
 618:                          // *              Reduce i-th column of matrix to tridiagonal form
 619:                          // *
 620:                          for (K = KDN + 1; K >= 2; K +=  - 1)
 621:                          {
 622:                              J1 = J1 + KDN;
 623:                              J2 = J2 + KDN;
 624:                              // *
 625:                              if (NR > 0)
 626:                              {
 627:                                  // *
 628:                                  // *                    generate plane rotations to annihilate nonzero
 629:                                  // *                    elements which have been created outside the band
 630:                                  // *
 631:                                  this._dlargv.Run(NR, ref AB, KD1+(J1 - KD1) * LDAB + o_ab, INCA, ref WORK, J1 + o_work, KD1, ref D, J1 + o_d
 632:                                                   , KD1);
 633:                                  // *
 634:                                  // *                    apply plane rotations from one side
 635:                                  // *
 636:                                  // *
 637:                                  // *                    Dependent on the the number of diagonals either
 638:                                  // *                    DLARTV or DROT is used
 639:                                  // *
 640:                                  if (NR > 2 * KD - 1)
 641:                                  {
 642:                                      for (L = 1; L <= KD - 1; L++)
 643:                                      {
 644:                                          this._dlartv.Run(NR, ref AB, KD1 - L+(J1 - KD1 + L) * LDAB + o_ab, INCA, ref AB, KD1 - L + 1+(J1 - KD1 + L) * LDAB + o_ab, INCA, D, J1 + o_d
 645:                                                           , WORK, J1 + o_work, KD1);
 646:                                      }
 647:                                  }
 648:                                  else
 649:                                  {
 650:                                      JEND = J1 + KD1 * (NR - 1);
 651:                                      for (JINC = J1; (KD1 >= 0) ? (JINC <= JEND) : (JINC >= JEND); JINC += KD1)
 652:                                      {
 653:                                          this._drot.Run(KDM1, ref AB, KD+(JINC - KD) * LDAB + o_ab, INCX, ref AB, KD1+(JINC - KD) * LDAB + o_ab, INCX, D[JINC + o_d]
 654:                                                         , WORK[JINC + o_work]);
 655:                                      }
 656:                                  }
 657:                                  // *
 658:                              }
 659:                              // *
 660:                              if (K > 2)
 661:                              {
 662:                                  if (K <= N - I + 1)
 663:                                  {
 664:                                      // *
 665:                                      // *                       generate plane rotation to annihilate a(i+k-1,i)
 666:                                      // *                       within the band
 667:                                      // *
 668:                                      this._dlartg.Run(AB[K - 1+I * LDAB + o_ab], AB[K+I * LDAB + o_ab], ref D[I + K - 1 + o_d], ref WORK[I + K - 1 + o_work], ref TEMP);
 669:                                      AB[K - 1+I * LDAB + o_ab] = TEMP;
 670:                                      // *
 671:                                      // *                       apply rotation from the left
 672:                                      // *
 673:                                      this._drot.Run(K - 3, ref AB, K - 2+(I + 1) * LDAB + o_ab, LDAB - 1, ref AB, K - 1+(I + 1) * LDAB + o_ab, LDAB - 1, D[I + K - 1 + o_d]
 674:                                                     , WORK[I + K - 1 + o_work]);
 675:                                  }
 676:                                  NR = NR + 1;
 677:                                  J1 = J1 - KDN - 1;
 678:                              }
 679:                              // *
 680:                              // *                 apply plane rotations from both sides to diagonal
 681:                              // *                 blocks
 682:                              // *
 683:                              if (NR > 0)
 684:                              {
 685:                                  this._dlar2v.Run(NR, ref AB, 1+(J1 - 1) * LDAB + o_ab, ref AB, 1+J1 * LDAB + o_ab, ref AB, 2+(J1 - 1) * LDAB + o_ab, INCA, D, J1 + o_d
 686:                                                   , WORK, J1 + o_work, KD1);
 687:                              }
 688:                              // *
 689:                              // *                 apply plane rotations from the right
 690:                              // *
 691:                              // *
 692:                              // *                    Dependent on the the number of diagonals either
 693:                              // *                    DLARTV or DROT is used
 694:                              // *
 695:                              if (NR > 0)
 696:                              {
 697:                                  if (NR > 2 * KD - 1)
 698:                                  {
 699:                                      for (L = 1; L <= KD - 1; L++)
 700:                                      {
 701:                                          if (J2 + L > N)
 702:                                          {
 703:                                              NRT = NR - 1;
 704:                                          }
 705:                                          else
 706:                                          {
 707:                                              NRT = NR;
 708:                                          }
 709:                                          if (NRT > 0)
 710:                                          {
 711:                                              this._dlartv.Run(NRT, ref AB, L + 2+(J1 - 1) * LDAB + o_ab, INCA, ref AB, L + 1+J1 * LDAB + o_ab, INCA, D, J1 + o_d
 712:                                                               , WORK, J1 + o_work, KD1);
 713:                                          }
 714:                                      }
 715:                                  }
 716:                                  else
 717:                                  {
 718:                                      J1END = J1 + KD1 * (NR - 2);
 719:                                      if (J1END >= J1)
 720:                                      {
 721:                                          for (J1INC = J1; (KD1 >= 0) ? (J1INC <= J1END) : (J1INC >= J1END); J1INC += KD1)
 722:                                          {
 723:                                              this._drot.Run(KDM1, ref AB, 3+(J1INC - 1) * LDAB + o_ab, 1, ref AB, 2+J1INC * LDAB + o_ab, 1, D[J1INC + o_d]
 724:                                                             , WORK[J1INC + o_work]);
 725:                                          }
 726:                                      }
 727:                                      LEND = Math.Min(KDM1, N - J2);
 728:                                      LAST = J1END + KD1;
 729:                                      if (LEND > 0)
 730:                                      {
 731:                                          this._drot.Run(LEND, ref AB, 3+(LAST - 1) * LDAB + o_ab, 1, ref AB, 2+LAST * LDAB + o_ab, 1, D[LAST + o_d]
 732:                                                         , WORK[LAST + o_work]);
 733:                                      }
 734:                                  }
 735:                              }
 736:                              // *
 737:                              // *
 738:                              // *
 739:                              if (WANTQ)
 740:                              {
 741:                                  // *
 742:                                  // *                    accumulate product of plane rotations in Q
 743:                                  // *
 744:                                  if (INITQ)
 745:                                  {
 746:                                      // *
 747:                                      // *                 take advantage of the fact that Q was
 748:                                      // *                 initially the Identity matrix
 749:                                      // *
 750:                                      IQEND = Math.Max(IQEND, J2);
 751:                                      I2 = Math.Max(0, K - 3);
 752:                                      IQAEND = 1 + I * KD;
 753:                                      if (K == 2) IQAEND = IQAEND + KD;
 754:                                      IQAEND = Math.Min(IQAEND, IQEND);
 755:                                      for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
 756:                                      {
 757:                                          IBL = I - I2 / KDM1;
 758:                                          I2 = I2 + 1;
 759:                                          IQB = Math.Max(1, J - IBL);
 760:                                          NQ = 1 + IQAEND - IQB;
 761:                                          IQAEND = Math.Min(IQAEND + KD, IQEND);
 762:                                          this._drot.Run(NQ, ref Q, IQB+(J - 1) * LDQ + o_q, 1, ref Q, IQB+J * LDQ + o_q, 1, D[J + o_d]
 763:                                                         , WORK[J + o_work]);
 764:                                      }
 765:                                  }
 766:                                  else
 767:                                  {
 768:                                      // *
 769:                                      for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
 770:                                      {
 771:                                          this._drot.Run(N, ref Q, 1+(J - 1) * LDQ + o_q, 1, ref Q, 1+J * LDQ + o_q, 1, D[J + o_d]
 772:                                                         , WORK[J + o_work]);
 773:                                      }
 774:                                  }
 775:                              }
 776:                              // *
 777:                              if (J2 + KDN > N)
 778:                              {
 779:                                  // *
 780:                                  // *                    adjust J2 to keep within the bounds of the matrix
 781:                                  // *
 782:                                  NR = NR - 1;
 783:                                  J2 = J2 - KDN - 1;
 784:                              }
 785:                              // *
 786:                              for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
 787:                              {
 788:                                  // *
 789:                                  // *                    create nonzero element a(j+kd,j-1) outside the
 790:                                  // *                    band and store it in WORK
 791:                                  // *
 792:                                  WORK[J + KD + o_work] = WORK[J + o_work] * AB[KD1+J * LDAB + o_ab];
 793:                                  AB[KD1+J * LDAB + o_ab] = D[J + o_d] * AB[KD1+J * LDAB + o_ab];
 794:                              }
 795:                          }
 796:                      }
 797:                  }
 798:                  // *
 799:                  if (KD > 0)
 800:                  {
 801:                      // *
 802:                      // *           copy off-diagonal elements to E
 803:                      // *
 804:                      for (I = 1; I <= N - 1; I++)
 805:                      {
 806:                          E[I + o_e] = AB[2+I * LDAB + o_ab];
 807:                      }
 808:                  }
 809:                  else
 810:                  {
 811:                      // *
 812:                      // *           set E to zero if original matrix was diagonal
 813:                      // *
 814:                      for (I = 1; I <= N - 1; I++)
 815:                      {
 816:                          E[I + o_e] = ZERO;
 817:                      }
 818:                  }
 819:                  // *
 820:                  // *        copy diagonal elements to D
 821:                  // *
 822:                  for (I = 1; I <= N; I++)
 823:                  {
 824:                      D[I + o_d] = AB[1+I * LDAB + o_ab];
 825:                  }
 826:              }
 827:              // *
 828:              return;
 829:              // *
 830:              // *     End of DSBTRD
 831:              // *
 832:   
 833:              #endregion
 834:   
 835:          }
 836:      }
 837:  }