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 routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DGELQF computes an LQ factorization of a real M-by-N matrix A:
  27:      /// A = L * Q.
  28:      /// 
  29:      ///</summary>
  30:      public class DGELQF
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          DGELQ2 _dgelq2; DLARFB _dlarfb; DLARFT _dlarft; XERBLA _xerbla; ILAENV _ilaenv; 
  37:   
  38:          #endregion
  39:   
  40:   
  41:          #region Fields
  42:          
  43:          bool LQUERY = false; int I = 0; int IB = 0; int IINFO = 0; int IWS = 0; int K = 0; int LDWORK = 0; int LWKOPT = 0; 
  44:          int NB = 0;int NBMIN = 0; int NX = 0; 
  45:   
  46:          #endregion
  47:   
  48:          public DGELQF(DGELQ2 dgelq2, DLARFB dlarfb, DLARFT dlarft, XERBLA xerbla, ILAENV ilaenv)
  49:          {
  50:      
  51:   
  52:              #region Set Dependencies
  53:              
  54:              this._dgelq2 = dgelq2; this._dlarfb = dlarfb; this._dlarft = dlarft; this._xerbla = xerbla; this._ilaenv = ilaenv; 
  55:   
  56:              #endregion
  57:   
  58:          }
  59:      
  60:          public DGELQF()
  61:          {
  62:      
  63:   
  64:              #region Dependencies (Initialization)
  65:              
  66:              LSAME lsame = new LSAME();
  67:              XERBLA xerbla = new XERBLA();
  68:              DLAMC3 dlamc3 = new DLAMC3();
  69:              DLAPY2 dlapy2 = new DLAPY2();
  70:              DNRM2 dnrm2 = new DNRM2();
  71:              DSCAL dscal = new DSCAL();
  72:              DCOPY dcopy = new DCOPY();
  73:              IEEECK ieeeck = new IEEECK();
  74:              IPARMQ iparmq = new IPARMQ();
  75:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  76:              DGER dger = new DGER(xerbla);
  77:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  78:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  79:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  80:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  81:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  82:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  83:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  84:              DGELQ2 dgelq2 = new DGELQ2(dlarf, dlarfg, xerbla);
  85:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  86:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  87:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
  88:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  89:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
  90:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  91:   
  92:              #endregion
  93:   
  94:   
  95:              #region Set Dependencies
  96:              
  97:              this._dgelq2 = dgelq2; this._dlarfb = dlarfb; this._dlarft = dlarft; this._xerbla = xerbla; this._ilaenv = ilaenv; 
  98:   
  99:              #endregion
 100:   
 101:          }
 102:          /// <summary>
 103:          /// Purpose
 104:          /// =======
 105:          /// 
 106:          /// DGELQF computes an LQ factorization of a real M-by-N matrix A:
 107:          /// A = L * Q.
 108:          /// 
 109:          ///</summary>
 110:          /// <param name="M">
 111:          /// (input) INTEGER
 112:          /// The number of rows of the matrix A.  M .GE. 0.
 113:          ///</param>
 114:          /// <param name="N">
 115:          /// (input) INTEGER
 116:          /// The number of columns of the matrix A.  N .GE. 0.
 117:          ///</param>
 118:          /// <param name="A">
 119:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 120:          /// On entry, the M-by-N matrix A.
 121:          /// On exit, the elements on and below the diagonal of the array
 122:          /// contain the m-by-min(m,n) lower trapezoidal matrix L (L is
 123:          /// lower triangular if m .LE. n); the elements above the diagonal,
 124:          /// with the array TAU, represent the orthogonal matrix Q as a
 125:          /// product of elementary reflectors (see Further Details).
 126:          ///</param>
 127:          /// <param name="LDA">
 128:          /// (input) INTEGER
 129:          /// The leading dimension of the array A.  LDA .GE. max(1,M).
 130:          ///</param>
 131:          /// <param name="TAU">
 132:          /// (output) DOUBLE PRECISION array, dimension (min(M,N))
 133:          /// The scalar factors of the elementary reflectors (see Further
 134:          /// Details).
 135:          ///</param>
 136:          /// <param name="WORK">
 137:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 138:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 139:          ///</param>
 140:          /// <param name="LWORK">
 141:          /// (input) INTEGER
 142:          /// The dimension of the array WORK.  LWORK .GE. max(1,M).
 143:          /// For optimum performance LWORK .GE. M*NB, where NB is the
 144:          /// optimal blocksize.
 145:          /// 
 146:          /// If LWORK = -1, then a workspace query is assumed; the routine
 147:          /// only calculates the optimal size of the WORK array, returns
 148:          /// this value as the first entry of the WORK array, and no error
 149:          /// message related to LWORK is issued by XERBLA.
 150:          ///</param>
 151:          /// <param name="INFO">
 152:          /// (output) INTEGER
 153:          /// = 0:  successful exit
 154:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 155:          ///</param>
 156:          public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau, ref double[] WORK, int offset_work
 157:                           , int LWORK, ref int INFO)
 158:          {
 159:   
 160:              #region Array Index Correction
 161:              
 162:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_work = -1 + offset_work; 
 163:   
 164:              #endregion
 165:   
 166:   
 167:              #region Prolog
 168:              
 169:              // *
 170:              // *  -- LAPACK routine (version 3.1) --
 171:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 172:              // *     November 2006
 173:              // *
 174:              // *     .. Scalar Arguments ..
 175:              // *     ..
 176:              // *     .. Array Arguments ..
 177:              // *     ..
 178:              // *
 179:              // *  Purpose
 180:              // *  =======
 181:              // *
 182:              // *  DGELQF computes an LQ factorization of a real M-by-N matrix A:
 183:              // *  A = L * Q.
 184:              // *
 185:              // *  Arguments
 186:              // *  =========
 187:              // *
 188:              // *  M       (input) INTEGER
 189:              // *          The number of rows of the matrix A.  M >= 0.
 190:              // *
 191:              // *  N       (input) INTEGER
 192:              // *          The number of columns of the matrix A.  N >= 0.
 193:              // *
 194:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 195:              // *          On entry, the M-by-N matrix A.
 196:              // *          On exit, the elements on and below the diagonal of the array
 197:              // *          contain the m-by-min(m,n) lower trapezoidal matrix L (L is
 198:              // *          lower triangular if m <= n); the elements above the diagonal,
 199:              // *          with the array TAU, represent the orthogonal matrix Q as a
 200:              // *          product of elementary reflectors (see Further Details).
 201:              // *
 202:              // *  LDA     (input) INTEGER
 203:              // *          The leading dimension of the array A.  LDA >= max(1,M).
 204:              // *
 205:              // *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
 206:              // *          The scalar factors of the elementary reflectors (see Further
 207:              // *          Details).
 208:              // *
 209:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 210:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 211:              // *
 212:              // *  LWORK   (input) INTEGER
 213:              // *          The dimension of the array WORK.  LWORK >= max(1,M).
 214:              // *          For optimum performance LWORK >= M*NB, where NB is the
 215:              // *          optimal blocksize.
 216:              // *
 217:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 218:              // *          only calculates the optimal size of the WORK array, returns
 219:              // *          this value as the first entry of the WORK array, and no error
 220:              // *          message related to LWORK is issued by XERBLA.
 221:              // *
 222:              // *  INFO    (output) INTEGER
 223:              // *          = 0:  successful exit
 224:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 225:              // *
 226:              // *  Further Details
 227:              // *  ===============
 228:              // *
 229:              // *  The matrix Q is represented as a product of elementary reflectors
 230:              // *
 231:              // *     Q = H(k) . . . H(2) H(1), where k = min(m,n).
 232:              // *
 233:              // *  Each H(i) has the form
 234:              // *
 235:              // *     H(i) = I - tau * v * v'
 236:              // *
 237:              // *  where tau is a real scalar, and v is a real vector with
 238:              // *  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n),
 239:              // *  and tau in TAU(i).
 240:              // *
 241:              // *  =====================================================================
 242:              // *
 243:              // *     .. Local Scalars ..
 244:              // *     ..
 245:              // *     .. External Subroutines ..
 246:              // *     ..
 247:              // *     .. Intrinsic Functions ..
 248:              //      INTRINSIC          MAX, MIN;
 249:              // *     ..
 250:              // *     .. External Functions ..
 251:              // *     ..
 252:              // *     .. Executable Statements ..
 253:              // *
 254:              // *     Test the input arguments
 255:              // *
 256:   
 257:              #endregion
 258:   
 259:   
 260:              #region Body
 261:              
 262:              INFO = 0;
 263:              NB = this._ilaenv.Run(1, "DGELQF", " ", M, N,  - 1,  - 1);
 264:              LWKOPT = M * NB;
 265:              WORK[1 + o_work] = LWKOPT;
 266:              LQUERY = (LWORK ==  - 1);
 267:              if (M < 0)
 268:              {
 269:                  INFO =  - 1;
 270:              }
 271:              else
 272:              {
 273:                  if (N < 0)
 274:                  {
 275:                      INFO =  - 2;
 276:                  }
 277:                  else
 278:                  {
 279:                      if (LDA < Math.Max(1, M))
 280:                      {
 281:                          INFO =  - 4;
 282:                      }
 283:                      else
 284:                      {
 285:                          if (LWORK < Math.Max(1, M) && !LQUERY)
 286:                          {
 287:                              INFO =  - 7;
 288:                          }
 289:                      }
 290:                  }
 291:              }
 292:              if (INFO != 0)
 293:              {
 294:                  this._xerbla.Run("DGELQF",  - INFO);
 295:                  return;
 296:              }
 297:              else
 298:              {
 299:                  if (LQUERY)
 300:                  {
 301:                      return;
 302:                  }
 303:              }
 304:              // *
 305:              // *     Quick return if possible
 306:              // *
 307:              K = Math.Min(M, N);
 308:              if (K == 0)
 309:              {
 310:                  WORK[1 + o_work] = 1;
 311:                  return;
 312:              }
 313:              // *
 314:              NBMIN = 2;
 315:              NX = 0;
 316:              IWS = M;
 317:              if (NB > 1 && NB < K)
 318:              {
 319:                  // *
 320:                  // *        Determine when to cross over from blocked to unblocked code.
 321:                  // *
 322:                  NX = Math.Max(0, this._ilaenv.Run(3, "DGELQF", " ", M, N,  - 1,  - 1));
 323:                  if (NX < K)
 324:                  {
 325:                      // *
 326:                      // *           Determine if workspace is large enough for blocked code.
 327:                      // *
 328:                      LDWORK = M;
 329:                      IWS = LDWORK * NB;
 330:                      if (LWORK < IWS)
 331:                      {
 332:                          // *
 333:                          // *              Not enough workspace to use optimal NB:  reduce NB and
 334:                          // *              determine the minimum value of NB.
 335:                          // *
 336:                          NB = LWORK / LDWORK;
 337:                          NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGELQF", " ", M, N,  - 1,  - 1));
 338:                      }
 339:                  }
 340:              }
 341:              // *
 342:              if (NB >= NBMIN && NB < K && NX < K)
 343:              {
 344:                  // *
 345:                  // *        Use blocked code initially
 346:                  // *
 347:                  for (I = 1; (NB >= 0) ? (I <= K - NX) : (I >= K - NX); I += NB)
 348:                  {
 349:                      IB = Math.Min(K - I + 1, NB);
 350:                      // *
 351:                      // *           Compute the LQ factorization of the current block
 352:                      // *           A(i:i+ib-1,i:n)
 353:                      // *
 354:                      this._dgelq2.Run(IB, N - I + 1, ref A, I+I * LDA + o_a, LDA, ref TAU, I + o_tau, ref WORK, offset_work
 355:                                       , ref IINFO);
 356:                      if (I + IB <= M)
 357:                      {
 358:                          // *
 359:                          // *              Form the triangular factor of the block reflector
 360:                          // *              H = H(i) H(i+1) . . . H(i+ib-1)
 361:                          // *
 362:                          this._dlarft.Run("Forward", "Rowwise", N - I + 1, IB, ref A, I+I * LDA + o_a, LDA
 363:                                           , TAU, I + o_tau, ref WORK, offset_work, LDWORK);
 364:                          // *
 365:                          // *              Apply H to A(i+ib:m,i:n) from the right
 366:                          // *
 367:                          this._dlarfb.Run("Right", "No transpose", "Forward", "Rowwise", M - I - IB + 1, N - I + 1
 368:                                           , IB, A, I+I * LDA + o_a, LDA, WORK, offset_work, LDWORK, ref A, I + IB+I * LDA + o_a
 369:                                           , LDA, ref WORK, IB + 1 + o_work, LDWORK);
 370:                      }
 371:                  }
 372:              }
 373:              else
 374:              {
 375:                  I = 1;
 376:              }
 377:              // *
 378:              // *     Use unblocked code to factor the last or only block.
 379:              // *
 380:              if (I <= K)
 381:              {
 382:                  this._dgelq2.Run(M - I + 1, N - I + 1, ref A, I+I * LDA + o_a, LDA, ref TAU, I + o_tau, ref WORK, offset_work
 383:                                   , ref IINFO);
 384:              }
 385:              // *
 386:              WORK[1 + o_work] = IWS;
 387:              return;
 388:              // *
 389:              // *     End of DGELQF
 390:              // *
 391:   
 392:              #endregion
 393:   
 394:          }
 395:      }
 396:  }