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