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: /// DORGQR generates an M-by-N real matrix Q with orthonormal columns,
27: /// which is defined as the first N columns of a product of K elementary
28: /// reflectors of order M
29: ///
30: /// Q = H(1) H(2) . . . H(k)
31: ///
32: /// as returned by DGEQRF.
33: ///
34: ///</summary>
35: public class DORGQR
36: {
37:
38:
39: #region Dependencies
40:
41: DLARFB _dlarfb; DLARFT _dlarft; DORG2R _dorg2r; 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 DORGQR(DLARFB dlarfb, DLARFT dlarft, DORG2R dorg2r, XERBLA xerbla, ILAENV ilaenv)
54: {
55:
56:
57: #region Set Dependencies
58:
59: this._dlarfb = dlarfb; this._dlarft = dlarft; this._dorg2r = dorg2r; this._xerbla = xerbla; this._ilaenv = ilaenv;
60:
61: #endregion
62:
63: }
64:
65: public DORGQR()
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: DORG2R dorg2r = new DORG2R(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._dorg2r = dorg2r; this._xerbla = xerbla; this._ilaenv = ilaenv;
94:
95: #endregion
96:
97: }
98: /// <summary>
99: /// Purpose
100: /// =======
101: ///
102: /// DORGQR generates an M-by-N real matrix Q with orthonormal columns,
103: /// which is defined as the first N columns of a product of K elementary
104: /// reflectors of order M
105: ///
106: /// Q = H(1) H(2) . . . H(k)
107: ///
108: /// as returned by DGEQRF.
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 i-th column must contain the vector which
127: /// defines the elementary reflector H(i), for i = 1,2,...,k, as
128: /// returned by DGEQRF in the first 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 DGEQRF.
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: // * DORGQR generates an M-by-N real matrix Q with orthonormal columns,
188: // * which is defined as the first N columns of a product of K elementary
189: // * reflectors of order M
190: // *
191: // * Q = H(1) H(2) . . . H(k)
192: // *
193: // * as returned by DGEQRF.
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 i-th column must contain the vector which
210: // * defines the elementary reflector H(i), for i = 1,2,...,k, as
211: // * returned by DGEQRF in the first 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 DGEQRF.
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: NB = this._ilaenv.Run(1, "DORGQR", " ", M, N, K, - 1);
264: LWKOPT = Math.Max(1, N) * 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 || N > M)
274: {
275: INFO = - 2;
276: }
277: else
278: {
279: if (K < 0 || K > N)
280: {
281: INFO = - 3;
282: }
283: else
284: {
285: if (LDA < Math.Max(1, M))
286: {
287: INFO = - 5;
288: }
289: else
290: {
291: if (LWORK < Math.Max(1, N) && !LQUERY)
292: {
293: INFO = - 8;
294: }
295: }
296: }
297: }
298: }
299: if (INFO != 0)
300: {
301: this._xerbla.Run("DORGQR", - INFO);
302: return;
303: }
304: else
305: {
306: if (LQUERY)
307: {
308: return;
309: }
310: }
311: // *
312: // * Quick return if possible
313: // *
314: if (N <= 0)
315: {
316: WORK[1 + o_work] = 1;
317: return;
318: }
319: // *
320: NBMIN = 2;
321: NX = 0;
322: IWS = N;
323: if (NB > 1 && NB < K)
324: {
325: // *
326: // * Determine when to cross over from blocked to unblocked code.
327: // *
328: NX = Math.Max(0, this._ilaenv.Run(3, "DORGQR", " ", M, N, K, - 1));
329: if (NX < K)
330: {
331: // *
332: // * Determine if workspace is large enough for blocked code.
333: // *
334: LDWORK = N;
335: IWS = LDWORK * NB;
336: if (LWORK < IWS)
337: {
338: // *
339: // * Not enough workspace to use optimal NB: reduce NB and
340: // * determine the minimum value of NB.
341: // *
342: NB = LWORK / LDWORK;
343: NBMIN = Math.Max(2, this._ilaenv.Run(2, "DORGQR", " ", M, N, K, - 1));
344: }
345: }
346: }
347: // *
348: if (NB >= NBMIN && NB < K && NX < K)
349: {
350: // *
351: // * Use blocked code after the last block.
352: // * The first kk columns are handled by the block method.
353: // *
354: KI = ((K - NX - 1) / NB) * NB;
355: KK = Math.Min(K, KI + NB);
356: // *
357: // * Set A(1:kk,kk+1:n) to zero.
358: // *
359: for (J = KK + 1; J <= N; J++)
360: {
361: for (I = 1; I <= KK; I++)
362: {
363: A[I+J * LDA + o_a] = ZERO;
364: }
365: }
366: }
367: else
368: {
369: KK = 0;
370: }
371: // *
372: // * Use unblocked code for the last or only block.
373: // *
374: if (KK < N)
375: {
376: this._dorg2r.Run(M - KK, N - KK, K - KK, ref A, KK + 1+(KK + 1) * LDA + o_a, LDA, TAU, KK + 1 + o_tau
377: , ref WORK, offset_work, ref IINFO);
378: }
379: // *
380: if (KK > 0)
381: {
382: // *
383: // * Use blocked code
384: // *
385: for (I = KI + 1; ( - NB >= 0) ? (I <= 1) : (I >= 1); I += - NB)
386: {
387: IB = Math.Min(NB, K - I + 1);
388: if (I + IB <= N)
389: {
390: // *
391: // * Form the triangular factor of the block reflector
392: // * H = H(i) H(i+1) . . . H(i+ib-1)
393: // *
394: this._dlarft.Run("Forward", "Columnwise", M - I + 1, IB, ref A, I+I * LDA + o_a, LDA
395: , TAU, I + o_tau, ref WORK, offset_work, LDWORK);
396: // *
397: // * Apply H to A(i:m,i+ib:n) from the left
398: // *
399: this._dlarfb.Run("Left", "No transpose", "Forward", "Columnwise", M - I + 1, N - I - IB + 1
400: , IB, A, I+I * LDA + o_a, LDA, WORK, offset_work, LDWORK, ref A, I+(I + IB) * LDA + o_a
401: , LDA, ref WORK, IB + 1 + o_work, LDWORK);
402: }
403: // *
404: // * Apply H to rows i:m of current block
405: // *
406: this._dorg2r.Run(M - I + 1, IB, IB, ref A, I+I * LDA + o_a, LDA, TAU, I + o_tau
407: , ref WORK, offset_work, ref IINFO);
408: // *
409: // * Set rows 1:i-1 of current block to zero
410: // *
411: for (J = I; J <= I + IB - 1; J++)
412: {
413: for (L = 1; L <= I - 1; L++)
414: {
415: A[L+J * LDA + o_a] = ZERO;
416: }
417: }
418: }
419: }
420: // *
421: WORK[1 + o_work] = IWS;
422: return;
423: // *
424: // * End of DORGQR
425: // *
426:
427: #endregion
428:
429: }
430: }
431: }