`   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:      /// DLAQPS computes a step of QR factorization with column pivoting`
`  27:      /// of a real M-by-N matrix A by using Blas-3.  It tries to factorize`
`  28:      /// NB columns from A starting from the row OFFSET+1, and updates all`
`  29:      /// of the matrix with Blas-3 xGEMM.`
`  30:      /// `
`  31:      /// In some cases, due to catastrophic cancellations, it cannot`
`  32:      /// factorize NB columns.  Hence, the actual number of factorized`
`  33:      /// columns is returned in KB.`
`  34:      /// `
`  35:      /// Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.`
`  36:      /// `
`  37:      ///</summary>`
`  38:      public class DLAQPS`
`  39:      {`
`  40:      `
`  41:   `
`  42:          #region Dependencies`
`  43:          `
`  44:          DGEMM _dgemm; DGEMV _dgemv; DLARFG _dlarfg; DSWAP _dswap; IDAMAX _idamax; DLAMCH _dlamch; DNRM2 _dnrm2; `
`  45:   `
`  46:          #endregion`
`  47:   `
`  48:   `
`  49:          #region Fields`
`  50:          `
`  51:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int ITEMP = 0; int J = 0; int K = 0; int LASTRK = 0; `
`  52:          int LSTICC = 0;int PVT = 0; int RK = 0; double AKK = 0; double TEMP = 0; double TEMP2 = 0; double TOL3Z = 0; `
`  53:   `
`  54:          #endregion`
`  55:   `
`  56:          public DLAQPS(DGEMM dgemm, DGEMV dgemv, DLARFG dlarfg, DSWAP dswap, IDAMAX idamax, DLAMCH dlamch, DNRM2 dnrm2)`
`  57:          {`
`  58:      `
`  59:   `
`  60:              #region Set Dependencies`
`  61:              `
`  62:              this._dgemm = dgemm; this._dgemv = dgemv; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax; `
`  63:              this._dlamch = dlamch;this._dnrm2 = dnrm2; `
`  64:   `
`  65:              #endregion`
`  66:   `
`  67:          }`
`  68:      `
`  69:          public DLAQPS()`
`  70:          {`
`  71:      `
`  72:   `
`  73:              #region Dependencies (Initialization)`
`  74:              `
`  75:              LSAME lsame = new LSAME();`
`  76:              XERBLA xerbla = new XERBLA();`
`  77:              DLAMC3 dlamc3 = new DLAMC3();`
`  78:              DLAPY2 dlapy2 = new DLAPY2();`
`  79:              DNRM2 dnrm2 = new DNRM2();`
`  80:              DSCAL dscal = new DSCAL();`
`  81:              DSWAP dswap = new DSWAP();`
`  82:              IDAMAX idamax = new IDAMAX();`
`  83:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
`  84:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
`  85:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  86:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  87:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  88:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  89:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  90:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);`
`  91:   `
`  92:              #endregion`
`  93:   `
`  94:   `
`  95:              #region Set Dependencies`
`  96:              `
`  97:              this._dgemm = dgemm; this._dgemv = dgemv; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax; `
`  98:              this._dlamch = dlamch;this._dnrm2 = dnrm2; `
`  99:   `
` 100:              #endregion`
` 101:   `
` 102:          }`
` 103:          /// <summary>`
` 104:          /// Purpose`
` 105:          /// =======`
` 106:          /// `
` 107:          /// DLAQPS computes a step of QR factorization with column pivoting`
` 108:          /// of a real M-by-N matrix A by using Blas-3.  It tries to factorize`
` 109:          /// NB columns from A starting from the row OFFSET+1, and updates all`
` 110:          /// of the matrix with Blas-3 xGEMM.`
` 111:          /// `
` 112:          /// In some cases, due to catastrophic cancellations, it cannot`
` 113:          /// factorize NB columns.  Hence, the actual number of factorized`
` 114:          /// columns is returned in KB.`
` 115:          /// `
` 116:          /// Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.`
` 117:          /// `
` 118:          ///</summary>`
` 119:          /// <param name="M">`
` 120:          /// (input) INTEGER`
` 121:          /// The number of rows of the matrix A. M .GE. 0.`
` 122:          ///</param>`
` 123:          /// <param name="N">`
` 124:          /// (input) INTEGER`
` 125:          /// The number of columns of the matrix A. N .GE. 0`
` 126:          ///</param>`
` 127:          /// <param name="OFFSET">`
` 128:          /// (input) INTEGER`
` 129:          /// The number of rows of A that have been factorized in`
` 130:          /// previous steps.`
` 131:          ///</param>`
` 132:          /// <param name="NB">`
` 133:          /// (input) INTEGER`
` 134:          /// The number of columns to factorize.`
` 135:          ///</param>`
` 136:          /// <param name="KB">`
` 137:          /// (output) INTEGER`
` 138:          /// The number of columns actually factorized.`
` 139:          ///</param>`
` 140:          /// <param name="A">`
` 141:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 142:          /// On entry, the M-by-N matrix A.`
` 143:          /// On exit, block A(OFFSET+1:M,1:KB) is the triangular`
` 144:          /// factor obtained and block A(1:OFFSET,1:N) has been`
` 145:          /// accordingly pivoted, but no factorized.`
` 146:          /// The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has`
` 147:          /// been updated.`
` 148:          ///</param>`
` 149:          /// <param name="LDA">`
` 150:          /// (input) INTEGER`
` 151:          /// The leading dimension of the array A. LDA .GE. max(1,M).`
` 152:          ///</param>`
` 153:          /// <param name="JPVT">`
` 154:          /// (input/output) INTEGER array, dimension (N)`
` 155:          /// JPVT(I) = K .LE.=.GT. Column K of the full matrix A has been`
` 156:          /// permuted into position I in AP.`
` 157:          ///</param>`
` 158:          /// <param name="TAU">`
` 159:          /// (output) DOUBLE PRECISION array, dimension (KB)`
` 160:          /// The scalar factors of the elementary reflectors.`
` 161:          ///</param>`
` 162:          /// <param name="VN1">`
` 163:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 164:          /// The vector with the partial column norms.`
` 165:          ///</param>`
` 166:          /// <param name="VN2">`
` 167:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 168:          /// The vector with the exact column norms.`
` 169:          ///</param>`
` 170:          /// <param name="AUXV">`
` 171:          /// (input/output) DOUBLE PRECISION array, dimension (NB)`
` 172:          /// Auxiliar vector.`
` 173:          ///</param>`
` 174:          /// <param name="F">`
` 175:          /// (input/output) DOUBLE PRECISION array, dimension (LDF,NB)`
` 176:          /// Matrix F' = L*Y'*A.`
` 177:          ///</param>`
` 178:          /// <param name="LDF">`
` 179:          /// (input) INTEGER`
` 180:          /// The leading dimension of the array F. LDF .GE. max(1,N).`
` 181:          ///</param>`
` 182:          public void Run(int M, int N, int OFFSET, int NB, ref int KB, ref double[] A, int offset_a`
` 183:                           , int LDA, ref int[] JPVT, int offset_jpvt, ref double[] TAU, int offset_tau, ref double[] VN1, int offset_vn1, ref double[] VN2, int offset_vn2, ref double[] AUXV, int offset_auxv`
` 184:                           , ref double[] F, int offset_f, int LDF)`
` 185:          {`
` 186:   `
` 187:              #region Array Index Correction`
` 188:              `
` 189:               int o_a = -1 - LDA + offset_a;  int o_jpvt = -1 + offset_jpvt;  int o_tau = -1 + offset_tau; `
` 190:               int o_vn1 = -1 + offset_vn1; int o_vn2 = -1 + offset_vn2;  int o_auxv = -1 + offset_auxv; `
` 191:               int o_f = -1 - LDF + offset_f;`
` 192:   `
` 193:              #endregion`
` 194:   `
` 195:   `
` 196:              #region Prolog`
` 197:              `
` 198:              // *`
` 199:              // *  -- LAPACK auxiliary routine (version 3.1) --`
` 200:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 201:              // *     November 2006`
` 202:              // *`
` 203:              // *     .. Scalar Arguments ..`
` 204:              // *     ..`
` 205:              // *     .. Array Arguments ..`
` 206:              // *     ..`
` 207:              // *`
` 208:              // *  Purpose`
` 209:              // *  =======`
` 210:              // *`
` 211:              // *  DLAQPS computes a step of QR factorization with column pivoting`
` 212:              // *  of a real M-by-N matrix A by using Blas-3.  It tries to factorize`
` 213:              // *  NB columns from A starting from the row OFFSET+1, and updates all`
` 214:              // *  of the matrix with Blas-3 xGEMM.`
` 215:              // *`
` 216:              // *  In some cases, due to catastrophic cancellations, it cannot`
` 217:              // *  factorize NB columns.  Hence, the actual number of factorized`
` 218:              // *  columns is returned in KB.`
` 219:              // *`
` 220:              // *  Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.`
` 221:              // *`
` 222:              // *  Arguments`
` 223:              // *  =========`
` 224:              // *`
` 225:              // *  M       (input) INTEGER`
` 226:              // *          The number of rows of the matrix A. M >= 0.`
` 227:              // *`
` 228:              // *  N       (input) INTEGER`
` 229:              // *          The number of columns of the matrix A. N >= 0`
` 230:              // *`
` 231:              // *  OFFSET  (input) INTEGER`
` 232:              // *          The number of rows of A that have been factorized in`
` 233:              // *          previous steps.`
` 234:              // *`
` 235:              // *  NB      (input) INTEGER`
` 236:              // *          The number of columns to factorize.`
` 237:              // *`
` 238:              // *  KB      (output) INTEGER`
` 239:              // *          The number of columns actually factorized.`
` 240:              // *`
` 241:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 242:              // *          On entry, the M-by-N matrix A.`
` 243:              // *          On exit, block A(OFFSET+1:M,1:KB) is the triangular`
` 244:              // *          factor obtained and block A(1:OFFSET,1:N) has been`
` 245:              // *          accordingly pivoted, but no factorized.`
` 246:              // *          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has`
` 247:              // *          been updated.`
` 248:              // *`
` 249:              // *  LDA     (input) INTEGER`
` 250:              // *          The leading dimension of the array A. LDA >= max(1,M).`
` 251:              // *`
` 252:              // *  JPVT    (input/output) INTEGER array, dimension (N)`
` 253:              // *          JPVT(I) = K <==> Column K of the full matrix A has been`
` 254:              // *          permuted into position I in AP.`
` 255:              // *`
` 256:              // *  TAU     (output) DOUBLE PRECISION array, dimension (KB)`
` 257:              // *          The scalar factors of the elementary reflectors.`
` 258:              // *`
` 259:              // *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)`
` 260:              // *          The vector with the partial column norms.`
` 261:              // *`
` 262:              // *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)`
` 263:              // *          The vector with the exact column norms.`
` 264:              // *`
` 265:              // *  AUXV    (input/output) DOUBLE PRECISION array, dimension (NB)`
` 266:              // *          Auxiliar vector.`
` 267:              // *`
` 268:              // *  F       (input/output) DOUBLE PRECISION array, dimension (LDF,NB)`
` 269:              // *          Matrix F' = L*Y'*A.`
` 270:              // *`
` 271:              // *  LDF     (input) INTEGER`
` 272:              // *          The leading dimension of the array F. LDF >= max(1,N).`
` 273:              // *`
` 274:              // *  Further Details`
` 275:              // *  ===============`
` 276:              // *`
` 277:              // *  Based on contributions by`
` 278:              // *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain`
` 279:              // *    X. Sun, Computer Science Dept., Duke University, USA`
` 280:              // *`
` 281:              // *  Partial column norm updating strategy modified by`
` 282:              // *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,`
` 283:              // *    University of Zagreb, Croatia.`
` 284:              // *    June 2006.`
` 285:              // *  For more details see LAPACK Working Note 176.`
` 286:              // *  =====================================================================`
` 287:              // *`
` 288:              // *     .. Parameters ..`
` 289:              // *     ..`
` 290:              // *     .. Local Scalars ..`
` 291:              // *     ..`
` 292:              // *     .. External Subroutines ..`
` 293:              // *     ..`
` 294:              // *     .. Intrinsic Functions ..`
` 295:              //      INTRINSIC          ABS, DBLE, MAX, MIN, NINT, SQRT;`
` 296:              // *     ..`
` 297:              // *     .. External Functions ..`
` 298:              // *     ..`
` 299:              // *     .. Executable Statements ..`
` 300:              // *`
` 301:   `
` 302:              #endregion`
` 303:   `
` 304:   `
` 305:              #region Body`
` 306:              `
` 307:              LASTRK = Math.Min(M, N + OFFSET);`
` 308:              LSTICC = 0;`
` 309:              K = 0;`
` 310:              TOL3Z = Math.Sqrt(this._dlamch.Run("Epsilon"));`
` 311:              // *`
` 312:              // *     Beginning of while loop.`
` 313:              // *`
` 314:          LABEL10:;`
` 315:              if ((K < NB) && (LSTICC == 0))`
` 316:              {`
` 317:                  K = K + 1;`
` 318:                  RK = OFFSET + K;`
` 319:                  // *`
` 320:                  // *        Determine ith pivot column and swap if necessary`
` 321:                  // *`
` 322:                  PVT = (K - 1) + this._idamax.Run(N - K + 1, VN1, K + o_vn1, 1);`
` 323:                  if (PVT != K)`
` 324:                  {`
` 325:                      this._dswap.Run(M, ref A, 1+PVT * LDA + o_a, 1, ref A, 1+K * LDA + o_a, 1);`
` 326:                      this._dswap.Run(K - 1, ref F, PVT+1 * LDF + o_f, LDF, ref F, K+1 * LDF + o_f, LDF);`
` 327:                      ITEMP = JPVT[PVT + o_jpvt];`
` 328:                      JPVT[PVT + o_jpvt] = JPVT[K + o_jpvt];`
` 329:                      JPVT[K + o_jpvt] = ITEMP;`
` 330:                      VN1[PVT + o_vn1] = VN1[K + o_vn1];`
` 331:                      VN2[PVT + o_vn2] = VN2[K + o_vn2];`
` 332:                  }`
` 333:                  // *`
` 334:                  // *        Apply previous Householder reflectors to column K:`
` 335:                  // *        A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'.`
` 336:                  // *`
` 337:                  if (K > 1)`
` 338:                  {`
` 339:                      this._dgemv.Run("No transpose", M - RK + 1, K - 1,  - ONE, A, RK+1 * LDA + o_a, LDA`
` 340:                                      , F, K+1 * LDF + o_f, LDF, ONE, ref A, RK+K * LDA + o_a, 1);`
` 341:                  }`
` 342:                  // *`
` 343:                  // *        Generate elementary reflector H(k).`
` 344:                  // *`
` 345:                  if (RK < M)`
` 346:                  {`
` 347:                      this._dlarfg.Run(M - RK + 1, ref A[RK+K * LDA + o_a], ref A, RK + 1+K * LDA + o_a, 1, ref TAU[K + o_tau]);`
` 348:                  }`
` 349:                  else`
` 350:                  {`
` 351:                      this._dlarfg.Run(1, ref A[RK+K * LDA + o_a], ref A, RK+K * LDA + o_a, 1, ref TAU[K + o_tau]);`
` 352:                  }`
` 353:                  // *`
` 354:                  AKK = A[RK+K * LDA + o_a];`
` 355:                  A[RK+K * LDA + o_a] = ONE;`
` 356:                  // *`
` 357:                  // *        Compute Kth column of F:`
` 358:                  // *`
` 359:                  // *        Compute  F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K).`
` 360:                  // *`
` 361:                  if (K < N)`
` 362:                  {`
` 363:                      this._dgemv.Run("Transpose", M - RK + 1, N - K, TAU[K + o_tau], A, RK+(K + 1) * LDA + o_a, LDA`
` 364:                                      , A, RK+K * LDA + o_a, 1, ZERO, ref F, K + 1+K * LDF + o_f, 1);`
` 365:                  }`
` 366:                  // *`
` 367:                  // *        Padding F(1:K,K) with zeros.`
` 368:                  // *`
` 369:                  for (J = 1; J <= K; J++)`
` 370:                  {`
` 371:                      F[J+K * LDF + o_f] = ZERO;`
` 372:                  }`
` 373:                  // *`
` 374:                  // *        Incremental updating of F:`
` 375:                  // *        F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)'`
` 376:                  // *                    *A(RK:M,K).`
` 377:                  // *`
` 378:                  if (K > 1)`
` 379:                  {`
` 380:                      this._dgemv.Run("Transpose", M - RK + 1, K - 1,  - TAU[K + o_tau], A, RK+1 * LDA + o_a, LDA`
` 381:                                      , A, RK+K * LDA + o_a, 1, ZERO, ref AUXV, 1 + o_auxv, 1);`
` 382:                      // *`
` 383:                      this._dgemv.Run("No transpose", N, K - 1, ONE, F, 1+1 * LDF + o_f, LDF`
` 384:                                      , AUXV, 1 + o_auxv, 1, ONE, ref F, 1+K * LDF + o_f, 1);`
` 385:                  }`
` 386:                  // *`
` 387:                  // *        Update the current row of A:`
` 388:                  // *        A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'.`
` 389:                  // *`
` 390:                  if (K < N)`
` 391:                  {`
` 392:                      this._dgemv.Run("No transpose", N - K, K,  - ONE, F, K + 1+1 * LDF + o_f, LDF`
` 393:                                      , A, RK+1 * LDA + o_a, LDA, ONE, ref A, RK+(K + 1) * LDA + o_a, LDA);`
` 394:                  }`
` 395:                  // *`
` 396:                  // *        Update partial column norms.`
` 397:                  // *`
` 398:                  if (RK < LASTRK)`
` 399:                  {`
` 400:                      for (J = K + 1; J <= N; J++)`
` 401:                      {`
` 402:                          if (VN1[J + o_vn1] != ZERO)`
` 403:                          {`
` 404:                              // *`
` 405:                              // *                 NOTE: The following 4 lines follow from the analysis in`
` 406:                              // *                 Lapack Working Note 176.`
` 407:                              // *`
` 408:                              TEMP = Math.Abs(A[RK+J * LDA + o_a]) / VN1[J + o_vn1];`
` 409:                              TEMP = Math.Max(ZERO, (ONE + TEMP) * (ONE - TEMP));`
` 410:                              TEMP2 = TEMP * Math.Pow(VN1[J + o_vn1] / VN2[J + o_vn2],2);`
` 411:                              if (TEMP2 <= TOL3Z)`
` 412:                              {`
` 413:                                  VN2[J + o_vn2] = Convert.ToDouble(LSTICC);`
` 414:                                  LSTICC = J;`
` 415:                              }`
` 416:                              else`
` 417:                              {`
` 418:                                  VN1[J + o_vn1] = VN1[J + o_vn1] * Math.Sqrt(TEMP);`
` 419:                              }`
` 420:                          }`
` 421:                      }`
` 422:                  }`
` 423:                  // *`
` 424:                  A[RK+K * LDA + o_a] = AKK;`
` 425:                  // *`
` 426:                  // *        End of while loop.`
` 427:                  // *`
` 428:                  goto LABEL10;`
` 429:              }`
` 430:              KB = K;`
` 431:              RK = OFFSET + KB;`
` 432:              // *`
` 433:              // *     Apply the block reflector to the rest of the matrix:`
` 434:              // *     A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -`
` 435:              // *                         A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'.`
` 436:              // *`
` 437:              if (KB < Math.Min(N, M - OFFSET))`
` 438:              {`
` 439:                  this._dgemm.Run("No transpose", "Transpose", M - RK, N - KB, KB,  - ONE`
` 440:                                  , A, RK + 1+1 * LDA + o_a, LDA, F, KB + 1+1 * LDF + o_f, LDF, ONE, ref A, RK + 1+(KB + 1) * LDA + o_a`
` 441:                                  , LDA);`
` 442:              }`
` 443:              // *`
` 444:              // *     Recomputation of difficult columns.`
` 445:              // *`
` 446:          LABEL40:;`
` 447:              if (LSTICC > 0)`
` 448:              {`
` 449:                  ITEMP = (int)Math.Round(VN2[LSTICC + o_vn2]);`
` 450:                  VN1[LSTICC + o_vn1] = this._dnrm2.Run(M - RK, A, RK + 1+LSTICC * LDA + o_a, 1);`
` 451:                  // *`
` 452:                  // *        NOTE: The computation of VN1( LSTICC ) relies on the fact that `
` 453:                  // *        SNRM2 does not fail on vectors with norm below the value of`
` 454:                  // *        SQRT(DLAMCH('S')) `
` 455:                  // *`
` 456:                  VN2[LSTICC + o_vn2] = VN1[LSTICC + o_vn1];`
` 457:                  LSTICC = ITEMP;`
` 458:                  goto LABEL40;`
` 459:              }`
` 460:              // *`
` 461:              return;`
` 462:              // *`
` 463:              // *     End of DLAQPS`
` 464:              // *`
` 465:   `
` 466:              #endregion`
` 467:   `
` 468:          }`
` 469:      }`
` 470:  }`