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: /// DGGRQF computes a generalized RQ factorization of an M-by-N matrix A
27: /// and a P-by-N matrix B:
28: ///
29: /// A = R*Q, B = Z*T*Q,
30: ///
31: /// where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
32: /// matrix, and R and T assume one of the forms:
33: ///
34: /// if M .LE. N, R = ( 0 R12 ) M, or if M .GT. N, R = ( R11 ) M-N,
35: /// N-M M ( R21 ) N
36: /// N
37: ///
38: /// where R12 or R21 is upper triangular, and
39: ///
40: /// if P .GE. N, T = ( T11 ) N , or if P .LT. N, T = ( T11 T12 ) P,
41: /// ( 0 ) P-N P N-P
42: /// N
43: ///
44: /// where T11 is upper triangular.
45: ///
46: /// In particular, if B is square and nonsingular, the GRQ factorization
47: /// of A and B implicitly gives the RQ factorization of A*inv(B):
48: ///
49: /// A*inv(B) = (R*inv(T))*Z'
50: ///
51: /// where inv(B) denotes the inverse of the matrix B, and Z' denotes the
52: /// transpose of the matrix Z.
53: ///
54: ///</summary>
55: public class DGGRQF
56: {
57:
58:
59: #region Dependencies
60:
61: DGEQRF _dgeqrf; DGERQF _dgerqf; DORMRQ _dormrq; XERBLA _xerbla; ILAENV _ilaenv;
62:
63: #endregion
64:
65:
66: #region Fields
67:
68: bool LQUERY = false; int LOPT = 0; int LWKOPT = 0; int NB = 0; int NB1 = 0; int NB2 = 0; int NB3 = 0;
69:
70: #endregion
71:
72: public DGGRQF(DGEQRF dgeqrf, DGERQF dgerqf, DORMRQ dormrq, XERBLA xerbla, ILAENV ilaenv)
73: {
74:
75:
76: #region Set Dependencies
77:
78: this._dgeqrf = dgeqrf; this._dgerqf = dgerqf; this._dormrq = dormrq; this._xerbla = xerbla; this._ilaenv = ilaenv;
79:
80: #endregion
81:
82: }
83:
84: public DGGRQF()
85: {
86:
87:
88: #region Dependencies (Initialization)
89:
90: LSAME lsame = new LSAME();
91: XERBLA xerbla = new XERBLA();
92: DLAMC3 dlamc3 = new DLAMC3();
93: DLAPY2 dlapy2 = new DLAPY2();
94: DNRM2 dnrm2 = new DNRM2();
95: DSCAL dscal = new DSCAL();
96: DCOPY dcopy = new DCOPY();
97: IEEECK ieeeck = new IEEECK();
98: IPARMQ iparmq = new IPARMQ();
99: DGEMV dgemv = new DGEMV(lsame, xerbla);
100: DGER dger = new DGER(xerbla);
101: DLARF dlarf = new DLARF(dgemv, dger, lsame);
102: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
103: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
104: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
105: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
106: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
107: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
108: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
109: DGEMM dgemm = new DGEMM(lsame, xerbla);
110: DTRMM dtrmm = new DTRMM(lsame, xerbla);
111: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
112: DTRMV dtrmv = new DTRMV(lsame, xerbla);
113: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
114: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
115: DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
116: DGERQ2 dgerq2 = new DGERQ2(dlarf, dlarfg, xerbla);
117: DGERQF dgerqf = new DGERQF(dgerq2, dlarfb, dlarft, xerbla, ilaenv);
118: DORMR2 dormr2 = new DORMR2(lsame, dlarf, xerbla);
119: DORMRQ dormrq = new DORMRQ(lsame, ilaenv, dlarfb, dlarft, dormr2, xerbla);
120:
121: #endregion
122:
123:
124: #region Set Dependencies
125:
126: this._dgeqrf = dgeqrf; this._dgerqf = dgerqf; this._dormrq = dormrq; this._xerbla = xerbla; this._ilaenv = ilaenv;
127:
128: #endregion
129:
130: }
131: /// <summary>
132: /// Purpose
133: /// =======
134: ///
135: /// DGGRQF computes a generalized RQ factorization of an M-by-N matrix A
136: /// and a P-by-N matrix B:
137: ///
138: /// A = R*Q, B = Z*T*Q,
139: ///
140: /// where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
141: /// matrix, and R and T assume one of the forms:
142: ///
143: /// if M .LE. N, R = ( 0 R12 ) M, or if M .GT. N, R = ( R11 ) M-N,
144: /// N-M M ( R21 ) N
145: /// N
146: ///
147: /// where R12 or R21 is upper triangular, and
148: ///
149: /// if P .GE. N, T = ( T11 ) N , or if P .LT. N, T = ( T11 T12 ) P,
150: /// ( 0 ) P-N P N-P
151: /// N
152: ///
153: /// where T11 is upper triangular.
154: ///
155: /// In particular, if B is square and nonsingular, the GRQ factorization
156: /// of A and B implicitly gives the RQ factorization of A*inv(B):
157: ///
158: /// A*inv(B) = (R*inv(T))*Z'
159: ///
160: /// where inv(B) denotes the inverse of the matrix B, and Z' denotes the
161: /// transpose of the matrix Z.
162: ///
163: ///</summary>
164: /// <param name="M">
165: /// (input) INTEGER
166: /// The number of rows of the matrix A. M .GE. 0.
167: ///</param>
168: /// <param name="P">
169: /// (input) INTEGER
170: /// The number of rows of the matrix B. P .GE. 0.
171: ///</param>
172: /// <param name="N">
173: /// (input) INTEGER
174: /// The number of columns of the matrices A and B. N .GE. 0.
175: ///</param>
176: /// <param name="A">
177: /// = R*Q, B = Z*T*Q,
178: ///</param>
179: /// <param name="LDA">
180: /// (input) INTEGER
181: /// The leading dimension of the array A. LDA .GE. max(1,M).
182: ///</param>
183: /// <param name="TAUA">
184: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
185: /// The scalar factors of the elementary reflectors which
186: /// represent the orthogonal matrix Q (see Further Details).
187: ///</param>
188: /// <param name="B">
189: /// (input/output) DOUBLE PRECISION array, dimension (LDB,N)
190: /// On entry, the P-by-N matrix B.
191: /// On exit, the elements on and above the diagonal of the array
192: /// contain the min(P,N)-by-N upper trapezoidal matrix T (T is
193: /// upper triangular if P .GE. N); the elements below the diagonal,
194: /// with the array TAUB, represent the orthogonal matrix Z as a
195: /// product of elementary reflectors (see Further Details).
196: ///</param>
197: /// <param name="LDB">
198: /// (input) INTEGER
199: /// The leading dimension of the array B. LDB .GE. max(1,P).
200: ///</param>
201: /// <param name="TAUB">
202: /// (output) DOUBLE PRECISION array, dimension (min(P,N))
203: /// The scalar factors of the elementary reflectors which
204: /// represent the orthogonal matrix Z (see Further Details).
205: ///</param>
206: /// <param name="WORK">
207: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
208: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
209: ///</param>
210: /// <param name="LWORK">
211: /// (input) INTEGER
212: /// The dimension of the array WORK. LWORK .GE. max(1,N,M,P).
213: /// For optimum performance LWORK .GE. max(N,M,P)*max(NB1,NB2,NB3),
214: /// where NB1 is the optimal blocksize for the RQ factorization
215: /// of an M-by-N matrix, NB2 is the optimal blocksize for the
216: /// QR factorization of a P-by-N matrix, and NB3 is the optimal
217: /// blocksize for a call of DORMRQ.
218: ///
219: /// If LWORK = -1, then a workspace query is assumed; the routine
220: /// only calculates the optimal size of the WORK array, returns
221: /// this value as the first entry of the WORK array, and no error
222: /// message related to LWORK is issued by XERBLA.
223: ///</param>
224: /// <param name="INFO">
225: /// (output) INTEGER
226: /// = 0: successful exit
227: /// .LT. 0: if INF0= -i, the i-th argument had an illegal value.
228: ///</param>
229: public void Run(int M, int P, int N, ref double[] A, int offset_a, int LDA, ref double[] TAUA, int offset_taua
230: , ref double[] B, int offset_b, int LDB, ref double[] TAUB, int offset_taub, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
231: {
232:
233: #region Array Index Correction
234:
235: int o_a = -1 - LDA + offset_a; int o_taua = -1 + offset_taua; int o_b = -1 - LDB + offset_b;
236: int o_taub = -1 + offset_taub; int o_work = -1 + offset_work;
237:
238: #endregion
239:
240:
241: #region Prolog
242:
243: // *
244: // * -- LAPACK routine (version 3.1) --
245: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
246: // * November 2006
247: // *
248: // * .. Scalar Arguments ..
249: // * ..
250: // * .. Array Arguments ..
251: // * ..
252: // *
253: // * Purpose
254: // * =======
255: // *
256: // * DGGRQF computes a generalized RQ factorization of an M-by-N matrix A
257: // * and a P-by-N matrix B:
258: // *
259: // * A = R*Q, B = Z*T*Q,
260: // *
261: // * where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal
262: // * matrix, and R and T assume one of the forms:
263: // *
264: // * if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N,
265: // * N-M M ( R21 ) N
266: // * N
267: // *
268: // * where R12 or R21 is upper triangular, and
269: // *
270: // * if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P,
271: // * ( 0 ) P-N P N-P
272: // * N
273: // *
274: // * where T11 is upper triangular.
275: // *
276: // * In particular, if B is square and nonsingular, the GRQ factorization
277: // * of A and B implicitly gives the RQ factorization of A*inv(B):
278: // *
279: // * A*inv(B) = (R*inv(T))*Z'
280: // *
281: // * where inv(B) denotes the inverse of the matrix B, and Z' denotes the
282: // * transpose of the matrix Z.
283: // *
284: // * Arguments
285: // * =========
286: // *
287: // * M (input) INTEGER
288: // * The number of rows of the matrix A. M >= 0.
289: // *
290: // * P (input) INTEGER
291: // * The number of rows of the matrix B. P >= 0.
292: // *
293: // * N (input) INTEGER
294: // * The number of columns of the matrices A and B. N >= 0.
295: // *
296: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
297: // * On entry, the M-by-N matrix A.
298: // * On exit, if M <= N, the upper triangle of the subarray
299: // * A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R;
300: // * if M > N, the elements on and above the (M-N)-th subdiagonal
301: // * contain the M-by-N upper trapezoidal matrix R; the remaining
302: // * elements, with the array TAUA, represent the orthogonal
303: // * matrix Q as a product of elementary reflectors (see Further
304: // * Details).
305: // *
306: // * LDA (input) INTEGER
307: // * The leading dimension of the array A. LDA >= max(1,M).
308: // *
309: // * TAUA (output) DOUBLE PRECISION array, dimension (min(M,N))
310: // * The scalar factors of the elementary reflectors which
311: // * represent the orthogonal matrix Q (see Further Details).
312: // *
313: // * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
314: // * On entry, the P-by-N matrix B.
315: // * On exit, the elements on and above the diagonal of the array
316: // * contain the min(P,N)-by-N upper trapezoidal matrix T (T is
317: // * upper triangular if P >= N); the elements below the diagonal,
318: // * with the array TAUB, represent the orthogonal matrix Z as a
319: // * product of elementary reflectors (see Further Details).
320: // *
321: // * LDB (input) INTEGER
322: // * The leading dimension of the array B. LDB >= max(1,P).
323: // *
324: // * TAUB (output) DOUBLE PRECISION array, dimension (min(P,N))
325: // * The scalar factors of the elementary reflectors which
326: // * represent the orthogonal matrix Z (see Further Details).
327: // *
328: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
329: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
330: // *
331: // * LWORK (input) INTEGER
332: // * The dimension of the array WORK. LWORK >= max(1,N,M,P).
333: // * For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
334: // * where NB1 is the optimal blocksize for the RQ factorization
335: // * of an M-by-N matrix, NB2 is the optimal blocksize for the
336: // * QR factorization of a P-by-N matrix, and NB3 is the optimal
337: // * blocksize for a call of DORMRQ.
338: // *
339: // * If LWORK = -1, then a workspace query is assumed; the routine
340: // * only calculates the optimal size of the WORK array, returns
341: // * this value as the first entry of the WORK array, and no error
342: // * message related to LWORK is issued by XERBLA.
343: // *
344: // * INFO (output) INTEGER
345: // * = 0: successful exit
346: // * < 0: if INF0= -i, the i-th argument had an illegal value.
347: // *
348: // * Further Details
349: // * ===============
350: // *
351: // * The matrix Q is represented as a product of elementary reflectors
352: // *
353: // * Q = H(1) H(2) . . . H(k), where k = min(m,n).
354: // *
355: // * Each H(i) has the form
356: // *
357: // * H(i) = I - taua * v * v'
358: // *
359: // * where taua is a real scalar, and v is a real vector with
360: // * v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in
361: // * A(m-k+i,1:n-k+i-1), and taua in TAUA(i).
362: // * To form Q explicitly, use LAPACK subroutine DORGRQ.
363: // * To use Q to update another matrix, use LAPACK subroutine DORMRQ.
364: // *
365: // * The matrix Z is represented as a product of elementary reflectors
366: // *
367: // * Z = H(1) H(2) . . . H(k), where k = min(p,n).
368: // *
369: // * Each H(i) has the form
370: // *
371: // * H(i) = I - taub * v * v'
372: // *
373: // * where taub is a real scalar, and v is a real vector with
374: // * v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i),
375: // * and taub in TAUB(i).
376: // * To form Z explicitly, use LAPACK subroutine DORGQR.
377: // * To use Z to update another matrix, use LAPACK subroutine DORMQR.
378: // *
379: // * =====================================================================
380: // *
381: // * .. Local Scalars ..
382: // * ..
383: // * .. External Subroutines ..
384: // * ..
385: // * .. External Functions ..
386: // * ..
387: // * .. Intrinsic Functions ..
388: // INTRINSIC INT, MAX, MIN;
389: // * ..
390: // * .. Executable Statements ..
391: // *
392: // * Test the input parameters
393: // *
394:
395: #endregion
396:
397:
398: #region Body
399:
400: INFO = 0;
401: NB1 = this._ilaenv.Run(1, "DGERQF", " ", M, N, - 1, - 1);
402: NB2 = this._ilaenv.Run(1, "DGEQRF", " ", P, N, - 1, - 1);
403: NB3 = this._ilaenv.Run(1, "DORMRQ", " ", M, N, P, - 1);
404: NB = Math.Max(NB1, Math.Max(NB2, NB3));
405: LWKOPT = Math.Max(N, Math.Max(M, P)) * NB;
406: WORK[1 + o_work] = LWKOPT;
407: LQUERY = (LWORK == - 1);
408: if (M < 0)
409: {
410: INFO = - 1;
411: }
412: else
413: {
414: if (P < 0)
415: {
416: INFO = - 2;
417: }
418: else
419: {
420: if (N < 0)
421: {
422: INFO = - 3;
423: }
424: else
425: {
426: if (LDA < Math.Max(1, M))
427: {
428: INFO = - 5;
429: }
430: else
431: {
432: if (LDB < Math.Max(1, P))
433: {
434: INFO = - 8;
435: }
436: else
437: {
438: if (LWORK < Math.Max(1, Math.Max(M, Math.Max(P, N))) && !LQUERY)
439: {
440: INFO = - 11;
441: }
442: }
443: }
444: }
445: }
446: }
447: if (INFO != 0)
448: {
449: this._xerbla.Run("DGGRQF", - INFO);
450: return;
451: }
452: else
453: {
454: if (LQUERY)
455: {
456: return;
457: }
458: }
459: // *
460: // * RQ factorization of M-by-N matrix A: A = R*Q
461: // *
462: this._dgerqf.Run(M, N, ref A, offset_a, LDA, ref TAUA, offset_taua, ref WORK, offset_work
463: , LWORK, ref INFO);
464: LOPT = (int)WORK[1 + o_work];
465: // *
466: // * Update B := B*Q'
467: // *
468: this._dormrq.Run("Right", "Transpose", P, N, Math.Min(M, N), ref A, Math.Max(1, M - N + 1)+1 * LDA + o_a
469: , LDA, TAUA, offset_taua, ref B, offset_b, LDB, ref WORK, offset_work, LWORK
470: , ref INFO);
471: LOPT = Math.Max(LOPT, Convert.ToInt32(Math.Truncate(WORK[1 + o_work])));
472: // *
473: // * QR factorization of P-by-N matrix B: B = Z*T
474: // *
475: this._dgeqrf.Run(P, N, ref B, offset_b, LDB, ref TAUB, offset_taub, ref WORK, offset_work
476: , LWORK, ref INFO);
477: WORK[1 + o_work] = Math.Max(LOPT, Convert.ToInt32(Math.Truncate(WORK[1 + o_work])));
478: // *
479: return;
480: // *
481: // * End of DGGRQF
482: // *
483:
484: #endregion
485:
486: }
487: }
488: }