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 deprecated driver routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// This routine is deprecated and has been replaced by routine DGEQP3.
  27:      /// 
  28:      /// DGEQPF computes a QR factorization with column pivoting of a
  29:      /// real M-by-N matrix A: A*P = Q*R.
  30:      /// 
  31:      ///</summary>
  32:      public class DGEQPF
  33:      {
  34:      
  35:   
  36:          #region Dependencies
  37:          
  38:          DGEQR2 _dgeqr2; DLARF _dlarf; DLARFG _dlarfg; DORM2R _dorm2r; DSWAP _dswap; XERBLA _xerbla; IDAMAX _idamax; 
  39:          DLAMCH _dlamch;DNRM2 _dnrm2; 
  40:   
  41:          #endregion
  42:   
  43:   
  44:          #region Fields
  45:          
  46:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; int ITEMP = 0; int J = 0; int MA = 0; int MN = 0; 
  47:          int PVT = 0;double AII = 0; double TEMP = 0; double TEMP2 = 0; double TOL3Z = 0; 
  48:   
  49:          #endregion
  50:   
  51:          public DGEQPF(DGEQR2 dgeqr2, DLARF dlarf, DLARFG dlarfg, DORM2R dorm2r, DSWAP dswap, XERBLA xerbla, IDAMAX idamax, DLAMCH dlamch, DNRM2 dnrm2)
  52:          {
  53:      
  54:   
  55:              #region Set Dependencies
  56:              
  57:              this._dgeqr2 = dgeqr2; this._dlarf = dlarf; this._dlarfg = dlarfg; this._dorm2r = dorm2r; this._dswap = dswap; 
  58:              this._xerbla = xerbla;this._idamax = idamax; this._dlamch = dlamch; this._dnrm2 = dnrm2; 
  59:   
  60:              #endregion
  61:   
  62:          }
  63:      
  64:          public DGEQPF()
  65:          {
  66:      
  67:   
  68:              #region Dependencies (Initialization)
  69:              
  70:              LSAME lsame = new LSAME();
  71:              XERBLA xerbla = new XERBLA();
  72:              DLAMC3 dlamc3 = new DLAMC3();
  73:              DLAPY2 dlapy2 = new DLAPY2();
  74:              DNRM2 dnrm2 = new DNRM2();
  75:              DSCAL dscal = new DSCAL();
  76:              DSWAP dswap = new DSWAP();
  77:              IDAMAX idamax = new IDAMAX();
  78:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  79:              DGER dger = new DGER(xerbla);
  80:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  81:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  82:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  83:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  84:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  85:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  86:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  87:              DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
  88:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
  89:   
  90:              #endregion
  91:   
  92:   
  93:              #region Set Dependencies
  94:              
  95:              this._dgeqr2 = dgeqr2; this._dlarf = dlarf; this._dlarfg = dlarfg; this._dorm2r = dorm2r; this._dswap = dswap; 
  96:              this._xerbla = xerbla;this._idamax = idamax; this._dlamch = dlamch; this._dnrm2 = dnrm2; 
  97:   
  98:              #endregion
  99:   
 100:          }
 101:          /// <summary>
 102:          /// Purpose
 103:          /// =======
 104:          /// 
 105:          /// This routine is deprecated and has been replaced by routine DGEQP3.
 106:          /// 
 107:          /// DGEQPF computes a QR factorization with column pivoting of a
 108:          /// real M-by-N matrix A: A*P = Q*R.
 109:          /// 
 110:          ///</summary>
 111:          /// <param name="M">
 112:          /// (input) INTEGER
 113:          /// The number of rows of the matrix A. M .GE. 0.
 114:          ///</param>
 115:          /// <param name="N">
 116:          /// (input) INTEGER
 117:          /// The number of columns of the matrix A. N .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 the array contains the
 123:          /// min(M,N)-by-N upper triangular matrix R; the elements
 124:          /// below the diagonal, together with the array TAU,
 125:          /// represent the orthogonal matrix Q as a product of
 126:          /// min(m,n) elementary reflectors.
 127:          ///</param>
 128:          /// <param name="LDA">
 129:          /// (input) INTEGER
 130:          /// The leading dimension of the array A. LDA .GE. max(1,M).
 131:          ///</param>
 132:          /// <param name="JPVT">
 133:          /// (input/output) INTEGER array, dimension (N)
 134:          /// On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
 135:          /// to the front of A*P (a leading column); if JPVT(i) = 0,
 136:          /// the i-th column of A is a free column.
 137:          /// On exit, if JPVT(i) = k, then the i-th column of A*P
 138:          /// was the k-th column of A.
 139:          ///</param>
 140:          /// <param name="TAU">
 141:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 142:          /// The scalar factors of the elementary reflectors.
 143:          ///</param>
 144:          /// <param name="WORK">
 145:          /// (workspace) DOUBLE PRECISION array, dimension (3*N)
 146:          ///</param>
 147:          /// <param name="INFO">
 148:          /// (output) INTEGER
 149:          /// = 0:  successful exit
 150:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 151:          ///</param>
 152:          public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref int[] JPVT, int offset_jpvt, ref double[] TAU, int offset_tau
 153:                           , ref double[] WORK, int offset_work, ref int INFO)
 154:          {
 155:   
 156:              #region Array Index Correction
 157:              
 158:               int o_a = -1 - LDA + offset_a;  int o_jpvt = -1 + offset_jpvt;  int o_tau = -1 + offset_tau; 
 159:               int o_work = -1 + offset_work;
 160:   
 161:              #endregion
 162:   
 163:   
 164:              #region Prolog
 165:              
 166:              // *
 167:              // *  -- LAPACK deprecated driver routine (version 3.1) --
 168:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 169:              // *     November 2006
 170:              // *
 171:              // *     .. Scalar Arguments ..
 172:              // *     ..
 173:              // *     .. Array Arguments ..
 174:              // *     ..
 175:              // *
 176:              // *  Purpose
 177:              // *  =======
 178:              // *
 179:              // *  This routine is deprecated and has been replaced by routine DGEQP3.
 180:              // *
 181:              // *  DGEQPF computes a QR factorization with column pivoting of a
 182:              // *  real M-by-N matrix A: A*P = Q*R.
 183:              // *
 184:              // *  Arguments
 185:              // *  =========
 186:              // *
 187:              // *  M       (input) INTEGER
 188:              // *          The number of rows of the matrix A. M >= 0.
 189:              // *
 190:              // *  N       (input) INTEGER
 191:              // *          The number of columns of the matrix A. N >= 0
 192:              // *
 193:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 194:              // *          On entry, the M-by-N matrix A.
 195:              // *          On exit, the upper triangle of the array contains the
 196:              // *          min(M,N)-by-N upper triangular matrix R; the elements
 197:              // *          below the diagonal, together with the array TAU,
 198:              // *          represent the orthogonal matrix Q as a product of
 199:              // *          min(m,n) elementary reflectors.
 200:              // *
 201:              // *  LDA     (input) INTEGER
 202:              // *          The leading dimension of the array A. LDA >= max(1,M).
 203:              // *
 204:              // *  JPVT    (input/output) INTEGER array, dimension (N)
 205:              // *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
 206:              // *          to the front of A*P (a leading column); if JPVT(i) = 0,
 207:              // *          the i-th column of A is a free column.
 208:              // *          On exit, if JPVT(i) = k, then the i-th column of A*P
 209:              // *          was the k-th column of A.
 210:              // *
 211:              // *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
 212:              // *          The scalar factors of the elementary reflectors.
 213:              // *
 214:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
 215:              // *
 216:              // *  INFO    (output) INTEGER
 217:              // *          = 0:  successful exit
 218:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 219:              // *
 220:              // *  Further Details
 221:              // *  ===============
 222:              // *
 223:              // *  The matrix Q is represented as a product of elementary reflectors
 224:              // *
 225:              // *     Q = H(1) H(2) . . . H(n)
 226:              // *
 227:              // *  Each H(i) has the form
 228:              // *
 229:              // *     H = I - tau * v * v'
 230:              // *
 231:              // *  where tau is a real scalar, and v is a real vector with
 232:              // *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
 233:              // *
 234:              // *  The matrix P is represented in jpvt as follows: If
 235:              // *     jpvt(j) = i
 236:              // *  then the jth column of P is the ith canonical unit vector.
 237:              // *
 238:              // *  Partial column norm updating strategy modified by
 239:              // *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
 240:              // *    University of Zagreb, Croatia.
 241:              // *    June 2006.
 242:              // *  For more details see LAPACK Working Note 176.
 243:              // *
 244:              // *  =====================================================================
 245:              // *
 246:              // *     .. Parameters ..
 247:              // *     ..
 248:              // *     .. Local Scalars ..
 249:              // *     ..
 250:              // *     .. External Subroutines ..
 251:              // *     ..
 252:              // *     .. Intrinsic Functions ..
 253:              //      INTRINSIC          ABS, MAX, MIN, SQRT;
 254:              // *     ..
 255:              // *     .. External Functions ..
 256:              // *     ..
 257:              // *     .. Executable Statements ..
 258:              // *
 259:              // *     Test the input arguments
 260:              // *
 261:   
 262:              #endregion
 263:   
 264:   
 265:              #region Body
 266:              
 267:              INFO = 0;
 268:              if (M < 0)
 269:              {
 270:                  INFO =  - 1;
 271:              }
 272:              else
 273:              {
 274:                  if (N < 0)
 275:                  {
 276:                      INFO =  - 2;
 277:                  }
 278:                  else
 279:                  {
 280:                      if (LDA < Math.Max(1, M))
 281:                      {
 282:                          INFO =  - 4;
 283:                      }
 284:                  }
 285:              }
 286:              if (INFO != 0)
 287:              {
 288:                  this._xerbla.Run("DGEQPF",  - INFO);
 289:                  return;
 290:              }
 291:              // *
 292:              MN = Math.Min(M, N);
 293:              TOL3Z = Math.Sqrt(this._dlamch.Run("Epsilon"));
 294:              // *
 295:              // *     Move initial columns up front
 296:              // *
 297:              ITEMP = 1;
 298:              for (I = 1; I <= N; I++)
 299:              {
 300:                  if (JPVT[I + o_jpvt] != 0)
 301:                  {
 302:                      if (I != ITEMP)
 303:                      {
 304:                          this._dswap.Run(M, ref A, 1+I * LDA + o_a, 1, ref A, 1+ITEMP * LDA + o_a, 1);
 305:                          JPVT[I + o_jpvt] = JPVT[ITEMP + o_jpvt];
 306:                          JPVT[ITEMP + o_jpvt] = I;
 307:                      }
 308:                      else
 309:                      {
 310:                          JPVT[I + o_jpvt] = I;
 311:                      }
 312:                      ITEMP = ITEMP + 1;
 313:                  }
 314:                  else
 315:                  {
 316:                      JPVT[I + o_jpvt] = I;
 317:                  }
 318:              }
 319:              ITEMP = ITEMP - 1;
 320:              // *
 321:              // *     Compute the QR factorization and update remaining columns
 322:              // *
 323:              if (ITEMP > 0)
 324:              {
 325:                  MA = Math.Min(ITEMP, M);
 326:                  this._dgeqr2.Run(M, MA, ref A, offset_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
 327:                                   , ref INFO);
 328:                  if (MA < N)
 329:                  {
 330:                      this._dorm2r.Run("Left", "Transpose", M, N - MA, MA, ref A, offset_a
 331:                                       , LDA, TAU, offset_tau, ref A, 1+(MA + 1) * LDA + o_a, LDA, ref WORK, offset_work, ref INFO);
 332:                  }
 333:              }
 334:              // *
 335:              if (ITEMP < MN)
 336:              {
 337:                  // *
 338:                  // *        Initialize partial column norms. The first n elements of
 339:                  // *        work store the exact column norms.
 340:                  // *
 341:                  for (I = ITEMP + 1; I <= N; I++)
 342:                  {
 343:                      WORK[I + o_work] = this._dnrm2.Run(M - ITEMP, A, ITEMP + 1+I * LDA + o_a, 1);
 344:                      WORK[N + I + o_work] = WORK[I + o_work];
 345:                  }
 346:                  // *
 347:                  // *        Compute factorization
 348:                  // *
 349:                  for (I = ITEMP + 1; I <= MN; I++)
 350:                  {
 351:                      // *
 352:                      // *           Determine ith pivot column and swap if necessary
 353:                      // *
 354:                      PVT = (I - 1) + this._idamax.Run(N - I + 1, WORK, I + o_work, 1);
 355:                      // *
 356:                      if (PVT != I)
 357:                      {
 358:                          this._dswap.Run(M, ref A, 1+PVT * LDA + o_a, 1, ref A, 1+I * LDA + o_a, 1);
 359:                          ITEMP = JPVT[PVT + o_jpvt];
 360:                          JPVT[PVT + o_jpvt] = JPVT[I + o_jpvt];
 361:                          JPVT[I + o_jpvt] = ITEMP;
 362:                          WORK[PVT + o_work] = WORK[I + o_work];
 363:                          WORK[N + PVT + o_work] = WORK[N + I + o_work];
 364:                      }
 365:                      // *
 366:                      // *           Generate elementary reflector H(i)
 367:                      // *
 368:                      if (I < M)
 369:                      {
 370:                          this._dlarfg.Run(M - I + 1, ref A[I+I * LDA + o_a], ref A, I + 1+I * LDA + o_a, 1, ref TAU[I + o_tau]);
 371:                      }
 372:                      else
 373:                      {
 374:                          this._dlarfg.Run(1, ref A[M+M * LDA + o_a], ref A, M+M * LDA + o_a, 1, ref TAU[M + o_tau]);
 375:                      }
 376:                      // *
 377:                      if (I < N)
 378:                      {
 379:                          // *
 380:                          // *              Apply H(i) to A(i:m,i+1:n) from the left
 381:                          // *
 382:                          AII = A[I+I * LDA + o_a];
 383:                          A[I+I * LDA + o_a] = ONE;
 384:                          this._dlarf.Run("LEFT", M - I + 1, N - I, A, I+I * LDA + o_a, 1, TAU[I + o_tau]
 385:                                          , ref A, I+(I + 1) * LDA + o_a, LDA, ref WORK, 2 * N + 1 + o_work);
 386:                          A[I+I * LDA + o_a] = AII;
 387:                      }
 388:                      // *
 389:                      // *           Update partial column norms
 390:                      // *
 391:                      for (J = I + 1; J <= N; J++)
 392:                      {
 393:                          if (WORK[J + o_work] != ZERO)
 394:                          {
 395:                              // *
 396:                              // *                 NOTE: The following 4 lines follow from the analysis in
 397:                              // *                 Lapack Working Note 176.
 398:                              // *                 
 399:                              TEMP = Math.Abs(A[I+J * LDA + o_a]) / WORK[J + o_work];
 400:                              TEMP = Math.Max(ZERO, (ONE + TEMP) * (ONE - TEMP));
 401:                              TEMP2 = TEMP * Math.Pow(WORK[J + o_work] / WORK[N + J + o_work],2);
 402:                              if (TEMP2 <= TOL3Z)
 403:                              {
 404:                                  if (M - I > 0)
 405:                                  {
 406:                                      WORK[J + o_work] = this._dnrm2.Run(M - I, A, I + 1+J * LDA + o_a, 1);
 407:                                      WORK[N + J + o_work] = WORK[J + o_work];
 408:                                  }
 409:                                  else
 410:                                  {
 411:                                      WORK[J + o_work] = ZERO;
 412:                                      WORK[N + J + o_work] = ZERO;
 413:                                  }
 414:                              }
 415:                              else
 416:                              {
 417:                                  WORK[J + o_work] = WORK[J + o_work] * Math.Sqrt(TEMP);
 418:                              }
 419:                          }
 420:                      }
 421:                      // *
 422:                  }
 423:              }
 424:              return;
 425:              // *
 426:              // *     End of DGEQPF
 427:              // *
 428:   
 429:              #endregion
 430:   
 431:          }
 432:      }
 433:  }