`   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 auxiliary routine (version 3.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// November 2006`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)`
`  27:      /// matrix A so that elements below the k-th subdiagonal are zero. The`
`  28:      /// reduction is performed by an orthogonal similarity transformation`
`  29:      /// Q' * A * Q. The routine returns the matrices V and T which determine`
`  30:      /// Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.`
`  31:      /// `
`  32:      /// This is an auxiliary routine called by DGEHRD.`
`  33:      /// `
`  34:      ///</summary>`
`  35:      public class DLAHR2`
`  36:      {`
`  37:      `
`  38:   `
`  39:          #region Dependencies`
`  40:          `
`  41:          DAXPY _daxpy; DCOPY _dcopy; DGEMM _dgemm; DGEMV _dgemv; DLACPY _dlacpy; DLARFG _dlarfg; DSCAL _dscal; DTRMM _dtrmm; `
`  42:          DTRMV _dtrmv;`
`  43:   `
`  44:          #endregion`
`  45:   `
`  46:   `
`  47:          #region Fields`
`  48:          `
`  49:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; double EI = 0; `
`  50:   `
`  51:          #endregion`
`  52:   `
`  53:          public DLAHR2(DAXPY daxpy, DCOPY dcopy, DGEMM dgemm, DGEMV dgemv, DLACPY dlacpy, DLARFG dlarfg, DSCAL dscal, DTRMM dtrmm, DTRMV dtrmv)`
`  54:          {`
`  55:      `
`  56:   `
`  57:              #region Set Dependencies`
`  58:              `
`  59:              this._daxpy = daxpy; this._dcopy = dcopy; this._dgemm = dgemm; this._dgemv = dgemv; this._dlacpy = dlacpy; `
`  60:              this._dlarfg = dlarfg;this._dscal = dscal; this._dtrmm = dtrmm; this._dtrmv = dtrmv; `
`  61:   `
`  62:              #endregion`
`  63:   `
`  64:          }`
`  65:      `
`  66:          public DLAHR2()`
`  67:          {`
`  68:      `
`  69:   `
`  70:              #region Dependencies (Initialization)`
`  71:              `
`  72:              DAXPY daxpy = new DAXPY();`
`  73:              DCOPY dcopy = new DCOPY();`
`  74:              LSAME lsame = new LSAME();`
`  75:              XERBLA xerbla = new XERBLA();`
`  76:              DLAMC3 dlamc3 = new DLAMC3();`
`  77:              DLAPY2 dlapy2 = new DLAPY2();`
`  78:              DNRM2 dnrm2 = new DNRM2();`
`  79:              DSCAL dscal = new DSCAL();`
`  80:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
`  81:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
`  82:              DLACPY dlacpy = new DLACPY(lsame);`
`  83:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  84:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  85:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  86:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  87:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  88:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
`  89:              DTRMM dtrmm = new DTRMM(lsame, xerbla);`
`  90:              DTRMV dtrmv = new DTRMV(lsame, xerbla);`
`  91:   `
`  92:              #endregion`
`  93:   `
`  94:   `
`  95:              #region Set Dependencies`
`  96:              `
`  97:              this._daxpy = daxpy; this._dcopy = dcopy; this._dgemm = dgemm; this._dgemv = dgemv; this._dlacpy = dlacpy; `
`  98:              this._dlarfg = dlarfg;this._dscal = dscal; this._dtrmm = dtrmm; this._dtrmv = dtrmv; `
`  99:   `
` 100:              #endregion`
` 101:   `
` 102:          }`
` 103:          /// <summary>`
` 104:          /// Purpose`
` 105:          /// =======`
` 106:          /// `
` 107:          /// DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)`
` 108:          /// matrix A so that elements below the k-th subdiagonal are zero. The`
` 109:          /// reduction is performed by an orthogonal similarity transformation`
` 110:          /// Q' * A * Q. The routine returns the matrices V and T which determine`
` 111:          /// Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.`
` 112:          /// `
` 113:          /// This is an auxiliary routine called by DGEHRD.`
` 114:          /// `
` 115:          ///</summary>`
` 116:          /// <param name="N">`
` 117:          /// (input) INTEGER`
` 118:          /// The order of the matrix A.`
` 119:          ///</param>`
` 120:          /// <param name="K">`
` 121:          /// (input) INTEGER`
` 122:          /// The offset for the reduction. Elements below the k-th`
` 123:          /// subdiagonal in the first NB columns are reduced to zero.`
` 124:          /// K .LT. N.`
` 125:          ///</param>`
` 126:          /// <param name="NB">`
` 127:          /// (input) INTEGER`
` 128:          /// The number of columns to be reduced.`
` 129:          ///</param>`
` 130:          /// <param name="A">`
` 131:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1)`
` 132:          /// On entry, the n-by-(n-k+1) general matrix A.`
` 133:          /// On exit, the elements on and above the k-th subdiagonal in`
` 134:          /// the first NB columns are overwritten with the corresponding`
` 135:          /// elements of the reduced matrix; the elements below the k-th`
` 136:          /// subdiagonal, with the array TAU, represent the matrix Q as a`
` 137:          /// product of elementary reflectors. The other columns of A are`
` 138:          /// unchanged. See Further Details.`
` 139:          ///</param>`
` 140:          /// <param name="LDA">`
` 141:          /// (input) INTEGER`
` 142:          /// The leading dimension of the array A.  LDA .GE. max(1,N).`
` 143:          ///</param>`
` 144:          /// <param name="TAU">`
` 145:          /// (output) DOUBLE PRECISION array, dimension (NB)`
` 146:          /// The scalar factors of the elementary reflectors. See Further`
` 147:          /// Details.`
` 148:          ///</param>`
` 149:          /// <param name="T">`
` 150:          /// (output) DOUBLE PRECISION array, dimension (LDT,NB)`
` 151:          /// The upper triangular matrix T.`
` 152:          ///</param>`
` 153:          /// <param name="LDT">`
` 154:          /// (input) INTEGER`
` 155:          /// The leading dimension of the array T.  LDT .GE. NB.`
` 156:          ///</param>`
` 157:          /// <param name="Y">`
` 158:          /// (output) DOUBLE PRECISION array, dimension (LDY,NB)`
` 159:          /// The n-by-nb matrix Y.`
` 160:          ///</param>`
` 161:          /// <param name="LDY">`
` 162:          /// (input) INTEGER`
` 163:          /// The leading dimension of the array Y. LDY .GE. N.`
` 164:          ///</param>`
` 165:          public void Run(int N, int K, int NB, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau`
` 166:                           , ref double[] T, int offset_t, int LDT, ref double[] Y, int offset_y, int LDY)`
` 167:          {`
` 168:   `
` 169:              #region Array Index Correction`
` 170:              `
` 171:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_t = -1 - LDT + offset_t; `
` 172:               int o_y = -1 - LDY + offset_y;`
` 173:   `
` 174:              #endregion`
` 175:   `
` 176:   `
` 177:              #region Prolog`
` 178:              `
` 179:              // *`
` 180:              // *  -- LAPACK auxiliary routine (version 3.1) --`
` 181:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 182:              // *     November 2006`
` 183:              // *`
` 184:              // *     .. Scalar Arguments ..`
` 185:              // *     ..`
` 186:              // *     .. Array Arguments ..`
` 187:              // *     ..`
` 188:              // *`
` 189:              // *  Purpose`
` 190:              // *  =======`
` 191:              // *`
` 192:              // *  DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1)`
` 193:              // *  matrix A so that elements below the k-th subdiagonal are zero. The`
` 194:              // *  reduction is performed by an orthogonal similarity transformation`
` 195:              // *  Q' * A * Q. The routine returns the matrices V and T which determine`
` 196:              // *  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.`
` 197:              // *`
` 198:              // *  This is an auxiliary routine called by DGEHRD.`
` 199:              // *`
` 200:              // *  Arguments`
` 201:              // *  =========`
` 202:              // *`
` 203:              // *  N       (input) INTEGER`
` 204:              // *          The order of the matrix A.`
` 205:              // *`
` 206:              // *  K       (input) INTEGER`
` 207:              // *          The offset for the reduction. Elements below the k-th`
` 208:              // *          subdiagonal in the first NB columns are reduced to zero.`
` 209:              // *          K < N.`
` 210:              // *`
` 211:              // *  NB      (input) INTEGER`
` 212:              // *          The number of columns to be reduced.`
` 213:              // *`
` 214:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1)`
` 215:              // *          On entry, the n-by-(n-k+1) general matrix A.`
` 216:              // *          On exit, the elements on and above the k-th subdiagonal in`
` 217:              // *          the first NB columns are overwritten with the corresponding`
` 218:              // *          elements of the reduced matrix; the elements below the k-th`
` 219:              // *          subdiagonal, with the array TAU, represent the matrix Q as a`
` 220:              // *          product of elementary reflectors. The other columns of A are`
` 221:              // *          unchanged. See Further Details.`
` 222:              // *`
` 223:              // *  LDA     (input) INTEGER`
` 224:              // *          The leading dimension of the array A.  LDA >= max(1,N).`
` 225:              // *`
` 226:              // *  TAU     (output) DOUBLE PRECISION array, dimension (NB)`
` 227:              // *          The scalar factors of the elementary reflectors. See Further`
` 228:              // *          Details.`
` 229:              // *`
` 230:              // *  T       (output) DOUBLE PRECISION array, dimension (LDT,NB)`
` 231:              // *          The upper triangular matrix T.`
` 232:              // *`
` 233:              // *  LDT     (input) INTEGER`
` 234:              // *          The leading dimension of the array T.  LDT >= NB.`
` 235:              // *`
` 236:              // *  Y       (output) DOUBLE PRECISION array, dimension (LDY,NB)`
` 237:              // *          The n-by-nb matrix Y.`
` 238:              // *`
` 239:              // *  LDY     (input) INTEGER`
` 240:              // *          The leading dimension of the array Y. LDY >= N.`
` 241:              // *`
` 242:              // *  Further Details`
` 243:              // *  ===============`
` 244:              // *`
` 245:              // *  The matrix Q is represented as a product of nb elementary reflectors`
` 246:              // *`
` 247:              // *     Q = H(1) H(2) . . . H(nb).`
` 248:              // *`
` 249:              // *  Each H(i) has the form`
` 250:              // *`
` 251:              // *     H(i) = I - tau * v * v'`
` 252:              // *`
` 253:              // *  where tau is a real scalar, and v is a real vector with`
` 254:              // *  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in`
` 255:              // *  A(i+k+1:n,i), and tau in TAU(i).`
` 256:              // *`
` 257:              // *  The elements of the vectors v together form the (n-k+1)-by-nb matrix`
` 258:              // *  V which is needed, with T and Y, to apply the transformation to the`
` 259:              // *  unreduced part of the matrix, using an update of the form:`
` 260:              // *  A := (I - V*T*V') * (A - Y*V').`
` 261:              // *`
` 262:              // *  The contents of A on exit are illustrated by the following example`
` 263:              // *  with n = 7, k = 3 and nb = 2:`
` 264:              // *`
` 265:              // *     ( a   a   a   a   a )`
` 266:              // *     ( a   a   a   a   a )`
` 267:              // *     ( a   a   a   a   a )`
` 268:              // *     ( h   h   a   a   a )`
` 269:              // *     ( v1  h   a   a   a )`
` 270:              // *     ( v1  v2  a   a   a )`
` 271:              // *     ( v1  v2  a   a   a )`
` 272:              // *`
` 273:              // *  where a denotes an element of the original matrix A, h denotes a`
` 274:              // *  modified element of the upper Hessenberg matrix H, and vi denotes an`
` 275:              // *  element of the vector defining H(i).`
` 276:              // *`
` 277:              // *  This file is a slight modification of LAPACK-3.0's DLAHRD`
` 278:              // *  incorporating improvements proposed by Quintana-Orti and Van de`
` 279:              // *  Gejin. Note that the entries of A(1:K,2:NB) differ from those`
` 280:              // *  returned by the original LAPACK routine. This function is`
` 281:              // *  not backward compatible with LAPACK3.0.`
` 282:              // *`
` 283:              // *  =====================================================================`
` 284:              // *`
` 285:              // *     .. Parameters ..`
` 286:              // *     ..`
` 287:              // *     .. Local Scalars ..`
` 288:              // *     ..`
` 289:              // *     .. External Subroutines ..`
` 290:              // *     ..`
` 291:              // *     .. Intrinsic Functions ..`
` 292:              //      INTRINSIC          MIN;`
` 293:              // *     ..`
` 294:              // *     .. Executable Statements ..`
` 295:              // *`
` 296:              // *     Quick return if possible`
` 297:              // *`
` 298:   `
` 299:              #endregion`
` 300:   `
` 301:   `
` 302:              #region Body`
` 303:              `
` 304:              if (N <= 1) return;`
` 305:              // *`
` 306:              for (I = 1; I <= NB; I++)`
` 307:              {`
` 308:                  if (I > 1)`
` 309:                  {`
` 310:                      // *`
` 311:                      // *           Update A(K+1:N,I)`
` 312:                      // *`
` 313:                      // *           Update I-th column of A - Y * V'`
` 314:                      // *`
` 315:                      this._dgemv.Run("NO TRANSPOSE", N - K, I - 1,  - ONE, Y, K + 1+1 * LDY + o_y, LDY`
` 316:                                      , A, K + I - 1+1 * LDA + o_a, LDA, ONE, ref A, K + 1+I * LDA + o_a, 1);`
` 317:                      // *`
` 318:                      // *           Apply I - V * T' * V' to this column (call it b) from the`
` 319:                      // *           left, using the last column of T as workspace`
` 320:                      // *`
` 321:                      // *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)`
` 322:                      // *                    ( V2 )             ( b2 )`
` 323:                      // *`
` 324:                      // *           where V1 is unit lower triangular`
` 325:                      // *`
` 326:                      // *           w := V1' * b1`
` 327:                      // *`
` 328:                      this._dcopy.Run(I - 1, A, K + 1+I * LDA + o_a, 1, ref T, 1+NB * LDT + o_t, 1);`
` 329:                      this._dtrmv.Run("Lower", "Transpose", "UNIT", I - 1, A, K + 1+1 * LDA + o_a, LDA`
` 330:                                      , ref T, 1+NB * LDT + o_t, 1);`
` 331:                      // *`
` 332:                      // *           w := w + V2'*b2`
` 333:                      // *`
` 334:                      this._dgemv.Run("Transpose", N - K - I + 1, I - 1, ONE, A, K + I+1 * LDA + o_a, LDA`
` 335:                                      , A, K + I+I * LDA + o_a, 1, ONE, ref T, 1+NB * LDT + o_t, 1);`
` 336:                      // *`
` 337:                      // *           w := T'*w`
` 338:                      // *`
` 339:                      this._dtrmv.Run("Upper", "Transpose", "NON-UNIT", I - 1, T, offset_t, LDT`
` 340:                                      , ref T, 1+NB * LDT + o_t, 1);`
` 341:                      // *`
` 342:                      // *           b2 := b2 - V2*w`
` 343:                      // *`
` 344:                      this._dgemv.Run("NO TRANSPOSE", N - K - I + 1, I - 1,  - ONE, A, K + I+1 * LDA + o_a, LDA`
` 345:                                      , T, 1+NB * LDT + o_t, 1, ONE, ref A, K + I+I * LDA + o_a, 1);`
` 346:                      // *`
` 347:                      // *           b1 := b1 - V1*w`
` 348:                      // *`
` 349:                      this._dtrmv.Run("Lower", "NO TRANSPOSE", "UNIT", I - 1, A, K + 1+1 * LDA + o_a, LDA`
` 350:                                      , ref T, 1+NB * LDT + o_t, 1);`
` 351:                      this._daxpy.Run(I - 1,  - ONE, T, 1+NB * LDT + o_t, 1, ref A, K + 1+I * LDA + o_a, 1);`
` 352:                      // *`
` 353:                      A[K + I - 1+(I - 1) * LDA + o_a] = EI;`
` 354:                  }`
` 355:                  // *`
` 356:                  // *        Generate the elementary reflector H(I) to annihilate`
` 357:                  // *        A(K+I+1:N,I)`
` 358:                  // *`
` 359:                  this._dlarfg.Run(N - K - I + 1, ref A[K + I+I * LDA + o_a], ref A, Math.Min(K + I + 1, N)+I * LDA + o_a, 1, ref TAU[I + o_tau]);`
` 360:                  EI = A[K + I+I * LDA + o_a];`
` 361:                  A[K + I+I * LDA + o_a] = ONE;`
` 362:                  // *`
` 363:                  // *        Compute  Y(K+1:N,I)`
` 364:                  // *`
` 365:                  this._dgemv.Run("NO TRANSPOSE", N - K, N - K - I + 1, ONE, A, K + 1+(I + 1) * LDA + o_a, LDA`
` 366:                                  , A, K + I+I * LDA + o_a, 1, ZERO, ref Y, K + 1+I * LDY + o_y, 1);`
` 367:                  this._dgemv.Run("Transpose", N - K - I + 1, I - 1, ONE, A, K + I+1 * LDA + o_a, LDA`
` 368:                                  , A, K + I+I * LDA + o_a, 1, ZERO, ref T, 1+I * LDT + o_t, 1);`
` 369:                  this._dgemv.Run("NO TRANSPOSE", N - K, I - 1,  - ONE, Y, K + 1+1 * LDY + o_y, LDY`
` 370:                                  , T, 1+I * LDT + o_t, 1, ONE, ref Y, K + 1+I * LDY + o_y, 1);`
` 371:                  this._dscal.Run(N - K, TAU[I + o_tau], ref Y, K + 1+I * LDY + o_y, 1);`
` 372:                  // *`
` 373:                  // *        Compute T(1:I,I)`
` 374:                  // *`
` 375:                  this._dscal.Run(I - 1,  - TAU[I + o_tau], ref T, 1+I * LDT + o_t, 1);`
` 376:                  this._dtrmv.Run("Upper", "No Transpose", "NON-UNIT", I - 1, T, offset_t, LDT`
` 377:                                  , ref T, 1+I * LDT + o_t, 1);`
` 378:                  T[I+I * LDT + o_t] = TAU[I + o_tau];`
` 379:                  // *`
` 380:              }`
` 381:              A[K + NB+NB * LDA + o_a] = EI;`
` 382:              // *`
` 383:              // *     Compute Y(1:K,1:NB)`
` 384:              // *`
` 385:              this._dlacpy.Run("ALL", K, NB, A, 1+2 * LDA + o_a, LDA, ref Y, offset_y`
` 386:                               , LDY);`
` 387:              this._dtrmm.Run("RIGHT", "Lower", "NO TRANSPOSE", "UNIT", K, NB`
` 388:                              , ONE, A, K + 1+1 * LDA + o_a, LDA, ref Y, offset_y, LDY);`
` 389:              if (N > K + NB)`
` 390:              {`
` 391:                  this._dgemm.Run("NO TRANSPOSE", "NO TRANSPOSE", K, NB, N - K - NB, ONE`
` 392:                                  , A, 1+(2 + NB) * LDA + o_a, LDA, A, K + 1 + NB+1 * LDA + o_a, LDA, ONE, ref Y, offset_y`
` 393:                                  , LDY);`
` 394:              }`
` 395:              this._dtrmm.Run("RIGHT", "Upper", "NO TRANSPOSE", "NON-UNIT", K, NB`
` 396:                              , ONE, T, offset_t, LDT, ref Y, offset_y, LDY);`
` 397:              // *`
` 398:              return;`
` 399:              // *`
` 400:              // *     End of DLAHR2`
` 401:              // *`
` 402:   `
` 403:              #endregion`
` 404:   `
` 405:          }`
` 406:      }`
` 407:  }`