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