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:      /// DORGQL generates an M-by-N real matrix Q with orthonormal columns,
27:      /// which is defined as the last N columns of a product of K elementary
28:      /// reflectors of order M
29:      ///
30:      /// Q  =  H(k) . . . H(2) H(1)
31:      ///
32:      /// as returned by DGEQLF.
33:      ///
34:      ///</summary>
35:      public class DORGQL
36:      {
37:
38:
39:          #region Dependencies
40:
41:          DLARFB _dlarfb; DLARFT _dlarft; DORG2L _dorg2l; XERBLA _xerbla; ILAENV _ilaenv;
42:
43:          #endregion
44:
45:
46:          #region Fields
47:
48:          const double ZERO = 0.0E+0; bool LQUERY = false; int I = 0; int IB = 0; int IINFO = 0; int IWS = 0; int J = 0; int KK = 0;
49:          int L = 0;int LDWORK = 0; int LWKOPT = 0; int NB = 0; int NBMIN = 0; int NX = 0;
50:
51:          #endregion
52:
53:          public DORGQL(DLARFB dlarfb, DLARFT dlarft, DORG2L dorg2l, XERBLA xerbla, ILAENV ilaenv)
54:          {
55:
56:
57:              #region Set Dependencies
58:
59:              this._dlarfb = dlarfb; this._dlarft = dlarft; this._dorg2l = dorg2l; this._xerbla = xerbla; this._ilaenv = ilaenv;
60:
61:              #endregion
62:
63:          }
64:
65:          public DORGQL()
66:          {
67:
68:
69:              #region Dependencies (Initialization)
70:
71:              LSAME lsame = new LSAME();
72:              DCOPY dcopy = new DCOPY();
73:              XERBLA xerbla = new XERBLA();
74:              DSCAL dscal = new DSCAL();
75:              IEEECK ieeeck = new IEEECK();
76:              IPARMQ iparmq = new IPARMQ();
77:              DGEMM dgemm = new DGEMM(lsame, xerbla);
78:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
79:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
80:              DGEMV dgemv = new DGEMV(lsame, xerbla);
81:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
82:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
83:              DGER dger = new DGER(xerbla);
84:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
85:              DORG2L dorg2l = new DORG2L(dlarf, dscal, xerbla);
86:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
87:
88:              #endregion
89:
90:
91:              #region Set Dependencies
92:
93:              this._dlarfb = dlarfb; this._dlarft = dlarft; this._dorg2l = dorg2l; this._xerbla = xerbla; this._ilaenv = ilaenv;
94:
95:              #endregion
96:
97:          }
98:          /// <summary>
99:          /// Purpose
100:          /// =======
101:          ///
102:          /// DORGQL generates an M-by-N real matrix Q with orthonormal columns,
103:          /// which is defined as the last N columns of a product of K elementary
104:          /// reflectors of order M
105:          ///
106:          /// Q  =  H(k) . . . H(2) H(1)
107:          ///
108:          /// as returned by DGEQLF.
109:          ///
110:          ///</summary>
111:          /// <param name="M">
112:          /// (input) INTEGER
113:          /// The number of rows of the matrix Q. M .GE. 0.
114:          ///</param>
115:          /// <param name="N">
116:          /// (input) INTEGER
117:          /// The number of columns of the matrix Q. M .GE. N .GE. 0.
118:          ///</param>
119:          /// <param name="K">
120:          /// (input) INTEGER
121:          /// The number of elementary reflectors whose product defines the
122:          /// matrix Q. N .GE. K .GE. 0.
123:          ///</param>
124:          /// <param name="A">
125:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
126:          /// On entry, the (n-k+i)-th column must contain the vector which
127:          /// defines the elementary reflector H(i), for i = 1,2,...,k, as
128:          /// returned by DGEQLF in the last k columns of its array
129:          /// argument A.
130:          /// On exit, the M-by-N matrix Q.
131:          ///</param>
132:          /// <param name="LDA">
133:          /// (input) INTEGER
134:          /// The first dimension of the array A. LDA .GE. max(1,M).
135:          ///</param>
136:          /// <param name="TAU">
137:          /// (input) DOUBLE PRECISION array, dimension (K)
138:          /// TAU(i) must contain the scalar factor of the elementary
139:          /// reflector H(i), as returned by DGEQLF.
140:          ///</param>
141:          /// <param name="WORK">
142:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
143:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
144:          ///</param>
145:          /// <param name="LWORK">
146:          /// (input) INTEGER
147:          /// The dimension of the array WORK. LWORK .GE. max(1,N).
148:          /// For optimum performance LWORK .GE. N*NB, where NB is the
149:          /// optimal blocksize.
150:          ///
151:          /// If LWORK = -1, then a workspace query is assumed; the routine
152:          /// only calculates the optimal size of the WORK array, returns
153:          /// this value as the first entry of the WORK array, and no error
154:          /// message related to LWORK is issued by XERBLA.
155:          ///</param>
156:          /// <param name="INFO">
157:          /// (output) INTEGER
158:          /// = 0:  successful exit
159:          /// .LT. 0:  if INFO = -i, the i-th argument has an illegal value
160:          ///</param>
161:          public void Run(int M, int N, int K, ref double[] A, int offset_a, int LDA, double[] TAU, int offset_tau
162:                           , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
163:          {
164:
165:              #region Array Index Correction
166:
167:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_work = -1 + offset_work;
168:
169:              #endregion
170:
171:
172:              #region Prolog
173:
174:              // *
175:              // *  -- LAPACK routine (version 3.1) --
176:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
177:              // *     November 2006
178:              // *
179:              // *     .. Scalar Arguments ..
180:              // *     ..
181:              // *     .. Array Arguments ..
182:              // *     ..
183:              // *
184:              // *  Purpose
185:              // *  =======
186:              // *
187:              // *  DORGQL generates an M-by-N real matrix Q with orthonormal columns,
188:              // *  which is defined as the last N columns of a product of K elementary
189:              // *  reflectors of order M
190:              // *
191:              // *        Q  =  H(k) . . . H(2) H(1)
192:              // *
193:              // *  as returned by DGEQLF.
194:              // *
195:              // *  Arguments
196:              // *  =========
197:              // *
198:              // *  M       (input) INTEGER
199:              // *          The number of rows of the matrix Q. M >= 0.
200:              // *
201:              // *  N       (input) INTEGER
202:              // *          The number of columns of the matrix Q. M >= N >= 0.
203:              // *
204:              // *  K       (input) INTEGER
205:              // *          The number of elementary reflectors whose product defines the
206:              // *          matrix Q. N >= K >= 0.
207:              // *
208:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
209:              // *          On entry, the (n-k+i)-th column must contain the vector which
210:              // *          defines the elementary reflector H(i), for i = 1,2,...,k, as
211:              // *          returned by DGEQLF in the last k columns of its array
212:              // *          argument A.
213:              // *          On exit, the M-by-N matrix Q.
214:              // *
215:              // *  LDA     (input) INTEGER
216:              // *          The first dimension of the array A. LDA >= max(1,M).
217:              // *
218:              // *  TAU     (input) DOUBLE PRECISION array, dimension (K)
219:              // *          TAU(i) must contain the scalar factor of the elementary
220:              // *          reflector H(i), as returned by DGEQLF.
221:              // *
222:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
223:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
224:              // *
225:              // *  LWORK   (input) INTEGER
226:              // *          The dimension of the array WORK. LWORK >= max(1,N).
227:              // *          For optimum performance LWORK >= N*NB, where NB is the
228:              // *          optimal blocksize.
229:              // *
230:              // *          If LWORK = -1, then a workspace query is assumed; the routine
231:              // *          only calculates the optimal size of the WORK array, returns
232:              // *          this value as the first entry of the WORK array, and no error
233:              // *          message related to LWORK is issued by XERBLA.
234:              // *
235:              // *  INFO    (output) INTEGER
236:              // *          = 0:  successful exit
237:              // *          < 0:  if INFO = -i, the i-th argument has an illegal value
238:              // *
239:              // *  =====================================================================
240:              // *
241:              // *     .. Parameters ..
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:              LQUERY = (LWORK ==  - 1);
264:              if (M < 0)
265:              {
266:                  INFO =  - 1;
267:              }
268:              else
269:              {
270:                  if (N < 0 || N > M)
271:                  {
272:                      INFO =  - 2;
273:                  }
274:                  else
275:                  {
276:                      if (K < 0 || K > N)
277:                      {
278:                          INFO =  - 3;
279:                      }
280:                      else
281:                      {
282:                          if (LDA < Math.Max(1, M))
283:                          {
284:                              INFO =  - 5;
285:                          }
286:                      }
287:                  }
288:              }
289:              // *
290:              if (INFO == 0)
291:              {
292:                  if (N == 0)
293:                  {
294:                      LWKOPT = 1;
295:                  }
296:                  else
297:                  {
298:                      NB = this._ilaenv.Run(1, "DORGQL", " ", M, N, K,  - 1);
299:                      LWKOPT = N * NB;
300:                  }
301:                  WORK[1 + o_work] = LWKOPT;
302:                  // *
303:                  if (LWORK < Math.Max(1, N) && !LQUERY)
304:                  {
305:                      INFO =  - 8;
306:                  }
307:              }
308:              // *
309:              if (INFO != 0)
310:              {
311:                  this._xerbla.Run("DORGQL",  - INFO);
312:                  return;
313:              }
314:              else
315:              {
316:                  if (LQUERY)
317:                  {
318:                      return;
319:                  }
320:              }
321:              // *
322:              // *     Quick return if possible
323:              // *
324:              if (N <= 0)
325:              {
326:                  return;
327:              }
328:              // *
329:              NBMIN = 2;
330:              NX = 0;
331:              IWS = N;
332:              if (NB > 1 && NB < K)
333:              {
334:                  // *
335:                  // *        Determine when to cross over from blocked to unblocked code.
336:                  // *
337:                  NX = Math.Max(0, this._ilaenv.Run(3, "DORGQL", " ", M, N, K,  - 1));
338:                  if (NX < K)
339:                  {
340:                      // *
341:                      // *           Determine if workspace is large enough for blocked code.
342:                      // *
343:                      LDWORK = N;
344:                      IWS = LDWORK * NB;
345:                      if (LWORK < IWS)
346:                      {
347:                          // *
348:                          // *              Not enough workspace to use optimal NB:  reduce NB and
349:                          // *              determine the minimum value of NB.
350:                          // *
351:                          NB = LWORK / LDWORK;
352:                          NBMIN = Math.Max(2, this._ilaenv.Run(2, "DORGQL", " ", M, N, K,  - 1));
353:                      }
354:                  }
355:              }
356:              // *
357:              if (NB >= NBMIN && NB < K && NX < K)
358:              {
359:                  // *
360:                  // *        Use blocked code after the first block.
361:                  // *        The last kk columns are handled by the block method.
362:                  // *
363:                  KK = Math.Min(K, ((K - NX + NB - 1) / NB) * NB);
364:                  // *
365:                  // *        Set A(m-kk+1:m,1:n-kk) to zero.
366:                  // *
367:                  for (J = 1; J <= N - KK; J++)
368:                  {
369:                      for (I = M - KK + 1; I <= M; I++)
370:                      {
371:                          A[I+J * LDA + o_a] = ZERO;
372:                      }
373:                  }
374:              }
375:              else
376:              {
377:                  KK = 0;
378:              }
379:              // *
380:              // *     Use unblocked code for the first or only block.
381:              // *
382:              this._dorg2l.Run(M - KK, N - KK, K - KK, ref A, offset_a, LDA, TAU, offset_tau
383:                               , ref WORK, offset_work, ref IINFO);
384:              // *
385:              if (KK > 0)
386:              {
387:                  // *
388:                  // *        Use blocked code
389:                  // *
390:                  for (I = K - KK + 1; (NB >= 0) ? (I <= K) : (I >= K); I += NB)
391:                  {
392:                      IB = Math.Min(NB, K - I + 1);
393:                      if (N - K + I > 1)
394:                      {
395:                          // *
396:                          // *              Form the triangular factor of the block reflector
397:                          // *              H = H(i+ib-1) . . . H(i+1) H(i)
398:                          // *
399:                          this._dlarft.Run("Backward", "Columnwise", M - K + I + IB - 1, IB, ref A, 1+(N - K + I) * LDA + o_a, LDA
400:                                           , TAU, I + o_tau, ref WORK, offset_work, LDWORK);
401:                          // *
402:                          // *              Apply H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
403:                          // *
404:                          this._dlarfb.Run("Left", "No transpose", "Backward", "Columnwise", M - K + I + IB - 1, N - K + I - 1
405:                                           , IB, A, 1+(N - K + I) * LDA + o_a, LDA, WORK, offset_work, LDWORK, ref A, offset_a
406:                                           , LDA, ref WORK, IB + 1 + o_work, LDWORK);
407:                      }
408:                      // *
409:                      // *           Apply H to rows 1:m-k+i+ib-1 of current block
410:                      // *
411:                      this._dorg2l.Run(M - K + I + IB - 1, IB, IB, ref A, 1+(N - K + I) * LDA + o_a, LDA, TAU, I + o_tau
412:                                       , ref WORK, offset_work, ref IINFO);
413:                      // *
414:                      // *           Set rows m-k+i+ib:m of current block to zero
415:                      // *
416:                      for (J = N - K + I; J <= N - K + I + IB - 1; J++)
417:                      {
418:                          for (L = M - K + I + IB; L <= M; L++)
419:                          {
420:                              A[L+J * LDA + o_a] = ZERO;
421:                          }
422:                      }
423:                  }
424:              }
425:              // *
426:              WORK[1 + o_work] = IWS;
427:              return;
428:              // *
429:              // *     End of DORGQL
430:              // *
431:
432:              #endregion
433:
434:          }
435:      }
436:  }