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 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:  }