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