`   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:      /// DLAQP2 computes a QR factorization with column pivoting of`
`  27:      /// the block A(OFFSET+1:M,1:N).`
`  28:      /// The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.`
`  29:      /// `
`  30:      ///</summary>`
`  31:      public class DLAQP2`
`  32:      {`
`  33:      `
`  34:   `
`  35:          #region Dependencies`
`  36:          `
`  37:          DLARF _dlarf; DLARFG _dlarfg; DSWAP _dswap; IDAMAX _idamax; DLAMCH _dlamch; DNRM2 _dnrm2; `
`  38:   `
`  39:          #endregion`
`  40:   `
`  41:   `
`  42:          #region Fields`
`  43:          `
`  44:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; int ITEMP = 0; int J = 0; int MN = 0; int OFFPI = 0; `
`  45:          int PVT = 0;double AII = 0; double TEMP = 0; double TEMP2 = 0; double TOL3Z = 0; `
`  46:   `
`  47:          #endregion`
`  48:   `
`  49:          public DLAQP2(DLARF dlarf, DLARFG dlarfg, DSWAP dswap, IDAMAX idamax, DLAMCH dlamch, DNRM2 dnrm2)`
`  50:          {`
`  51:      `
`  52:   `
`  53:              #region Set Dependencies`
`  54:              `
`  55:              this._dlarf = dlarf; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax; this._dlamch = dlamch; `
`  56:              this._dnrm2 = dnrm2;`
`  57:   `
`  58:              #endregion`
`  59:   `
`  60:          }`
`  61:      `
`  62:          public DLAQP2()`
`  63:          {`
`  64:      `
`  65:   `
`  66:              #region Dependencies (Initialization)`
`  67:              `
`  68:              LSAME lsame = new LSAME();`
`  69:              XERBLA xerbla = new XERBLA();`
`  70:              DLAMC3 dlamc3 = new DLAMC3();`
`  71:              DLAPY2 dlapy2 = new DLAPY2();`
`  72:              DNRM2 dnrm2 = new DNRM2();`
`  73:              DSCAL dscal = new DSCAL();`
`  74:              DSWAP dswap = new DSWAP();`
`  75:              IDAMAX idamax = new IDAMAX();`
`  76:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
`  77:              DGER dger = new DGER(xerbla);`
`  78:              DLARF dlarf = new DLARF(dgemv, dger, lsame);`
`  79:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  80:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  81:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  82:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  83:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  84:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
`  85:   `
`  86:              #endregion`
`  87:   `
`  88:   `
`  89:              #region Set Dependencies`
`  90:              `
`  91:              this._dlarf = dlarf; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax; this._dlamch = dlamch; `
`  92:              this._dnrm2 = dnrm2;`
`  93:   `
`  94:              #endregion`
`  95:   `
`  96:          }`
`  97:          /// <summary>`
`  98:          /// Purpose`
`  99:          /// =======`
` 100:          /// `
` 101:          /// DLAQP2 computes a QR factorization with column pivoting of`
` 102:          /// the block A(OFFSET+1:M,1:N).`
` 103:          /// The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.`
` 104:          /// `
` 105:          ///</summary>`
` 106:          /// <param name="M">`
` 107:          /// (input) INTEGER`
` 108:          /// The number of rows of the matrix A. M .GE. 0.`
` 109:          ///</param>`
` 110:          /// <param name="N">`
` 111:          /// (input) INTEGER`
` 112:          /// The number of columns of the matrix A. N .GE. 0.`
` 113:          ///</param>`
` 114:          /// <param name="OFFSET">`
` 115:          /// (input) INTEGER`
` 116:          /// The number of rows of the matrix A that must be pivoted`
` 117:          /// but no factorized. OFFSET .GE. 0.`
` 118:          ///</param>`
` 119:          /// <param name="A">`
` 120:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 121:          /// On entry, the M-by-N matrix A.`
` 122:          /// On exit, the upper triangle of block A(OFFSET+1:M,1:N) is `
` 123:          /// the triangular factor obtained; the elements in block`
` 124:          /// A(OFFSET+1:M,1:N) below the diagonal, together with the`
` 125:          /// array TAU, represent the orthogonal matrix Q as a product of`
` 126:          /// elementary reflectors. Block A(1:OFFSET,1:N) has been`
` 127:          /// accordingly pivoted, but no factorized.`
` 128:          ///</param>`
` 129:          /// <param name="LDA">`
` 130:          /// (input) INTEGER`
` 131:          /// The leading dimension of the array A. LDA .GE. max(1,M).`
` 132:          ///</param>`
` 133:          /// <param name="JPVT">`
` 134:          /// (input/output) INTEGER array, dimension (N)`
` 135:          /// On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted`
` 136:          /// to the front of A*P (a leading column); if JPVT(i) = 0,`
` 137:          /// the i-th column of A is a free column.`
` 138:          /// On exit, if JPVT(i) = k, then the i-th column of A*P`
` 139:          /// was the k-th column of A.`
` 140:          ///</param>`
` 141:          /// <param name="TAU">`
` 142:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))`
` 143:          /// The scalar factors of the elementary reflectors.`
` 144:          ///</param>`
` 145:          /// <param name="VN1">`
` 146:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 147:          /// The vector with the partial column norms.`
` 148:          ///</param>`
` 149:          /// <param name="VN2">`
` 150:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 151:          /// The vector with the exact column norms.`
` 152:          ///</param>`
` 153:          /// <param name="WORK">`
` 154:          /// (workspace) DOUBLE PRECISION array, dimension (N)`
` 155:          ///</param>`
` 156:          public void Run(int M, int N, int OFFSET, ref double[] A, int offset_a, int LDA, ref int[] JPVT, int offset_jpvt`
` 157:                           , ref double[] TAU, int offset_tau, ref double[] VN1, int offset_vn1, ref double[] VN2, int offset_vn2, ref double[] WORK, int offset_work)`
` 158:          {`
` 159:   `
` 160:              #region Array Index Correction`
` 161:              `
` 162:               int o_a = -1 - LDA + offset_a;  int o_jpvt = -1 + offset_jpvt;  int o_tau = -1 + offset_tau; `
` 163:               int o_vn1 = -1 + offset_vn1; int o_vn2 = -1 + offset_vn2;  int o_work = -1 + offset_work; `
` 164:   `
` 165:              #endregion`
` 166:   `
` 167:   `
` 168:              #region Prolog`
` 169:              `
` 170:              // *`
` 171:              // *  -- LAPACK auxiliary routine (version 3.1) --`
` 172:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 173:              // *     November 2006`
` 174:              // *`
` 175:              // *     .. Scalar Arguments ..`
` 176:              // *     ..`
` 177:              // *     .. Array Arguments ..`
` 178:              // *     ..`
` 179:              // *`
` 180:              // *  Purpose`
` 181:              // *  =======`
` 182:              // *`
` 183:              // *  DLAQP2 computes a QR factorization with column pivoting of`
` 184:              // *  the block A(OFFSET+1:M,1:N).`
` 185:              // *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.`
` 186:              // *`
` 187:              // *  Arguments`
` 188:              // *  =========`
` 189:              // *`
` 190:              // *  M       (input) INTEGER`
` 191:              // *          The number of rows of the matrix A. M >= 0.`
` 192:              // *`
` 193:              // *  N       (input) INTEGER`
` 194:              // *          The number of columns of the matrix A. N >= 0.`
` 195:              // *`
` 196:              // *  OFFSET  (input) INTEGER`
` 197:              // *          The number of rows of the matrix A that must be pivoted`
` 198:              // *          but no factorized. OFFSET >= 0.`
` 199:              // *`
` 200:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 201:              // *          On entry, the M-by-N matrix A.`
` 202:              // *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is `
` 203:              // *          the triangular factor obtained; the elements in block`
` 204:              // *          A(OFFSET+1:M,1:N) below the diagonal, together with the`
` 205:              // *          array TAU, represent the orthogonal matrix Q as a product of`
` 206:              // *          elementary reflectors. Block A(1:OFFSET,1:N) has been`
` 207:              // *          accordingly pivoted, but no factorized.`
` 208:              // *`
` 209:              // *  LDA     (input) INTEGER`
` 210:              // *          The leading dimension of the array A. LDA >= max(1,M).`
` 211:              // *`
` 212:              // *  JPVT    (input/output) INTEGER array, dimension (N)`
` 213:              // *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted`
` 214:              // *          to the front of A*P (a leading column); if JPVT(i) = 0,`
` 215:              // *          the i-th column of A is a free column.`
` 216:              // *          On exit, if JPVT(i) = k, then the i-th column of A*P`
` 217:              // *          was the k-th column of A.`
` 218:              // *`
` 219:              // *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))`
` 220:              // *          The scalar factors of the elementary reflectors.`
` 221:              // *`
` 222:              // *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)`
` 223:              // *          The vector with the partial column norms.`
` 224:              // *`
` 225:              // *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)`
` 226:              // *          The vector with the exact column norms.`
` 227:              // *`
` 228:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)`
` 229:              // *`
` 230:              // *  Further Details`
` 231:              // *  ===============`
` 232:              // *`
` 233:              // *  Based on contributions by`
` 234:              // *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain`
` 235:              // *    X. Sun, Computer Science Dept., Duke University, USA`
` 236:              // *`
` 237:              // *  Partial column norm updating strategy modified by`
` 238:              // *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,`
` 239:              // *    University of Zagreb, Croatia.`
` 240:              // *    June 2006.`
` 241:              // *  For more details see LAPACK Working Note 176.`
` 242:              // *  =====================================================================`
` 243:              // *`
` 244:              // *     .. Parameters ..`
` 245:              // *     ..`
` 246:              // *     .. Local Scalars ..`
` 247:              // *     ..`
` 248:              // *     .. External Subroutines ..`
` 249:              // *     ..`
` 250:              // *     .. Intrinsic Functions ..`
` 251:              //      INTRINSIC          ABS, MAX, MIN, SQRT;`
` 252:              // *     ..`
` 253:              // *     .. External Functions ..`
` 254:              // *     ..`
` 255:              // *     .. Executable Statements ..`
` 256:              // *`
` 257:   `
` 258:              #endregion`
` 259:   `
` 260:   `
` 261:              #region Body`
` 262:              `
` 263:              MN = Math.Min(M - OFFSET, N);`
` 264:              TOL3Z = Math.Sqrt(this._dlamch.Run("Epsilon"));`
` 265:              // *`
` 266:              // *     Compute factorization.`
` 267:              // *`
` 268:              for (I = 1; I <= MN; I++)`
` 269:              {`
` 270:                  // *`
` 271:                  OFFPI = OFFSET + I;`
` 272:                  // *`
` 273:                  // *        Determine ith pivot column and swap if necessary.`
` 274:                  // *`
` 275:                  PVT = (I - 1) + this._idamax.Run(N - I + 1, VN1, I + o_vn1, 1);`
` 276:                  // *`
` 277:                  if (PVT != I)`
` 278:                  {`
` 279:                      this._dswap.Run(M, ref A, 1+PVT * LDA + o_a, 1, ref A, 1+I * LDA + o_a, 1);`
` 280:                      ITEMP = JPVT[PVT + o_jpvt];`
` 281:                      JPVT[PVT + o_jpvt] = JPVT[I + o_jpvt];`
` 282:                      JPVT[I + o_jpvt] = ITEMP;`
` 283:                      VN1[PVT + o_vn1] = VN1[I + o_vn1];`
` 284:                      VN2[PVT + o_vn2] = VN2[I + o_vn2];`
` 285:                  }`
` 286:                  // *`
` 287:                  // *        Generate elementary reflector H(i).`
` 288:                  // *`
` 289:                  if (OFFPI < M)`
` 290:                  {`
` 291:                      this._dlarfg.Run(M - OFFPI + 1, ref A[OFFPI+I * LDA + o_a], ref A, OFFPI + 1+I * LDA + o_a, 1, ref TAU[I + o_tau]);`
` 292:                  }`
` 293:                  else`
` 294:                  {`
` 295:                      this._dlarfg.Run(1, ref A[M+I * LDA + o_a], ref A, M+I * LDA + o_a, 1, ref TAU[I + o_tau]);`
` 296:                  }`
` 297:                  // *`
` 298:                  if (I < N)`
` 299:                  {`
` 300:                      // *`
` 301:                      // *           Apply H(i)' to A(offset+i:m,i+1:n) from the left.`
` 302:                      // *`
` 303:                      AII = A[OFFPI+I * LDA + o_a];`
` 304:                      A[OFFPI+I * LDA + o_a] = ONE;`
` 305:                      this._dlarf.Run("Left", M - OFFPI + 1, N - I, A, OFFPI+I * LDA + o_a, 1, TAU[I + o_tau]`
` 306:                                      , ref A, OFFPI+(I + 1) * LDA + o_a, LDA, ref WORK, 1 + o_work);`
` 307:                      A[OFFPI+I * LDA + o_a] = AII;`
` 308:                  }`
` 309:                  // *`
` 310:                  // *        Update partial column norms.`
` 311:                  // *`
` 312:                  for (J = I + 1; J <= N; J++)`
` 313:                  {`
` 314:                      if (VN1[J + o_vn1] != ZERO)`
` 315:                      {`
` 316:                          // *`
` 317:                          // *              NOTE: The following 4 lines follow from the analysis in`
` 318:                          // *              Lapack Working Note 176.`
` 319:                          // *`
` 320:                          TEMP = ONE - Math.Pow(Math.Abs(A[OFFPI+J * LDA + o_a]) / VN1[J + o_vn1],2);`
` 321:                          TEMP = Math.Max(TEMP, ZERO);`
` 322:                          TEMP2 = TEMP * Math.Pow(VN1[J + o_vn1] / VN2[J + o_vn2],2);`
` 323:                          if (TEMP2 <= TOL3Z)`
` 324:                          {`
` 325:                              if (OFFPI < M)`
` 326:                              {`
` 327:                                  VN1[J + o_vn1] = this._dnrm2.Run(M - OFFPI, A, OFFPI + 1+J * LDA + o_a, 1);`
` 328:                                  VN2[J + o_vn2] = VN1[J + o_vn1];`
` 329:                              }`
` 330:                              else`
` 331:                              {`
` 332:                                  VN1[J + o_vn1] = ZERO;`
` 333:                                  VN2[J + o_vn2] = ZERO;`
` 334:                              }`
` 335:                          }`
` 336:                          else`
` 337:                          {`
` 338:                              VN1[J + o_vn1] = VN1[J + o_vn1] * Math.Sqrt(TEMP);`
` 339:                          }`
` 340:                      }`
` 341:                  }`
` 342:                  // *`
` 343:              }`
` 344:              // *`
` 345:              return;`
` 346:              // *`
` 347:              // *     End of DLAQP2`
` 348:              // *`
` 349:   `
` 350:              #endregion`
` 351:   `
` 352:          }`
` 353:      }`
` 354:  }`