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: /// DGEQRF computes a QR factorization of a real M-by-N matrix A:
27: /// A = Q * R.
28: ///
29: ///</summary>
30: public class DGEQRF
31: {
32:
33:
34: #region Dependencies
35:
36: DGEQR2 _dgeqr2; DLARFB _dlarfb; DLARFT _dlarft; XERBLA _xerbla; ILAENV _ilaenv;
37:
38: #endregion
39:
40:
41: #region Fields
42:
43: bool LQUERY = false; int I = 0; int IB = 0; int IINFO = 0; int IWS = 0; int K = 0; int LDWORK = 0; int LWKOPT = 0;
44: int NB = 0;int NBMIN = 0; int NX = 0;
45:
46: #endregion
47:
48: public DGEQRF(DGEQR2 dgeqr2, DLARFB dlarfb, DLARFT dlarft, XERBLA xerbla, ILAENV ilaenv)
49: {
50:
51:
52: #region Set Dependencies
53:
54: this._dgeqr2 = dgeqr2; this._dlarfb = dlarfb; this._dlarft = dlarft; this._xerbla = xerbla; this._ilaenv = ilaenv;
55:
56: #endregion
57:
58: }
59:
60: public DGEQRF()
61: {
62:
63:
64: #region Dependencies (Initialization)
65:
66: LSAME lsame = new LSAME();
67: XERBLA xerbla = new XERBLA();
68: DLAMC3 dlamc3 = new DLAMC3();
69: DLAPY2 dlapy2 = new DLAPY2();
70: DNRM2 dnrm2 = new DNRM2();
71: DSCAL dscal = new DSCAL();
72: DCOPY dcopy = new DCOPY();
73: IEEECK ieeeck = new IEEECK();
74: IPARMQ iparmq = new IPARMQ();
75: DGEMV dgemv = new DGEMV(lsame, xerbla);
76: DGER dger = new DGER(xerbla);
77: DLARF dlarf = new DLARF(dgemv, dger, lsame);
78: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
79: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
80: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
81: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
82: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
83: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
84: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
85: DGEMM dgemm = new DGEMM(lsame, xerbla);
86: DTRMM dtrmm = new DTRMM(lsame, xerbla);
87: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
88: DTRMV dtrmv = new DTRMV(lsame, xerbla);
89: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
90: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
91:
92: #endregion
93:
94:
95: #region Set Dependencies
96:
97: this._dgeqr2 = dgeqr2; this._dlarfb = dlarfb; this._dlarft = dlarft; this._xerbla = xerbla; this._ilaenv = ilaenv;
98:
99: #endregion
100:
101: }
102: /// <summary>
103: /// Purpose
104: /// =======
105: ///
106: /// DGEQRF computes a QR factorization of a real M-by-N matrix A:
107: /// A = Q * R.
108: ///
109: ///</summary>
110: /// <param name="M">
111: /// (input) INTEGER
112: /// The number of rows of the matrix A. M .GE. 0.
113: ///</param>
114: /// <param name="N">
115: /// (input) INTEGER
116: /// The number of columns of the matrix A. N .GE. 0.
117: ///</param>
118: /// <param name="A">
119: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
120: /// On entry, the M-by-N matrix A.
121: /// On exit, the elements on and above the diagonal of the array
122: /// contain the min(M,N)-by-N upper trapezoidal matrix R (R is
123: /// upper triangular if m .GE. n); the elements below the diagonal,
124: /// with the array TAU, represent the orthogonal matrix Q as a
125: /// product of min(m,n) elementary reflectors (see Further
126: /// Details).
127: ///</param>
128: /// <param name="LDA">
129: /// (input) INTEGER
130: /// The leading dimension of the array A. LDA .GE. max(1,M).
131: ///</param>
132: /// <param name="TAU">
133: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
134: /// The scalar factors of the elementary reflectors (see Further
135: /// Details).
136: ///</param>
137: /// <param name="WORK">
138: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
139: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
140: ///</param>
141: /// <param name="LWORK">
142: /// (input) INTEGER
143: /// The dimension of the array WORK. LWORK .GE. max(1,N).
144: /// For optimum performance LWORK .GE. N*NB, where NB is
145: /// the optimal blocksize.
146: ///
147: /// If LWORK = -1, then a workspace query is assumed; the routine
148: /// only calculates the optimal size of the WORK array, returns
149: /// this value as the first entry of the WORK array, and no error
150: /// message related to LWORK is issued by XERBLA.
151: ///</param>
152: /// <param name="INFO">
153: /// (output) INTEGER
154: /// = 0: successful exit
155: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
156: ///</param>
157: public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau, ref double[] WORK, int offset_work
158: , int LWORK, ref int INFO)
159: {
160:
161: #region Array Index Correction
162:
163: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
164:
165: #endregion
166:
167:
168: #region Prolog
169:
170: // *
171: // * -- LAPACK routine (version 3.1) --
172: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
173: // * November 2006
174: // *
175: // * .. Scalar Arguments ..
176: // * ..
177: // * .. Array Arguments ..
178: // * ..
179: // *
180: // * Purpose
181: // * =======
182: // *
183: // * DGEQRF computes a QR factorization of a real M-by-N matrix A:
184: // * A = Q * R.
185: // *
186: // * Arguments
187: // * =========
188: // *
189: // * M (input) INTEGER
190: // * The number of rows of the matrix A. M >= 0.
191: // *
192: // * N (input) INTEGER
193: // * The number of columns of the matrix A. N >= 0.
194: // *
195: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
196: // * On entry, the M-by-N matrix A.
197: // * On exit, the elements on and above the diagonal of the array
198: // * contain the min(M,N)-by-N upper trapezoidal matrix R (R is
199: // * upper triangular if m >= n); the elements below the diagonal,
200: // * with the array TAU, represent the orthogonal matrix Q as a
201: // * product of min(m,n) elementary reflectors (see Further
202: // * Details).
203: // *
204: // * LDA (input) INTEGER
205: // * The leading dimension of the array A. LDA >= max(1,M).
206: // *
207: // * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
208: // * The scalar factors of the elementary reflectors (see Further
209: // * Details).
210: // *
211: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
212: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
213: // *
214: // * LWORK (input) INTEGER
215: // * The dimension of the array WORK. LWORK >= max(1,N).
216: // * For optimum performance LWORK >= N*NB, where NB is
217: // * the optimal blocksize.
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: // *
224: // * INFO (output) INTEGER
225: // * = 0: successful exit
226: // * < 0: if INFO = -i, the i-th argument had an illegal value
227: // *
228: // * Further Details
229: // * ===============
230: // *
231: // * The matrix Q is represented as a product of elementary reflectors
232: // *
233: // * Q = H(1) H(2) . . . H(k), where k = min(m,n).
234: // *
235: // * Each H(i) has the form
236: // *
237: // * H(i) = I - tau * v * v'
238: // *
239: // * where tau is a real scalar, and v is a real vector with
240: // * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
241: // * and tau in TAU(i).
242: // *
243: // * =====================================================================
244: // *
245: // * .. Local Scalars ..
246: // * ..
247: // * .. External Subroutines ..
248: // * ..
249: // * .. Intrinsic Functions ..
250: // INTRINSIC MAX, MIN;
251: // * ..
252: // * .. External Functions ..
253: // * ..
254: // * .. Executable Statements ..
255: // *
256: // * Test the input arguments
257: // *
258:
259: #endregion
260:
261:
262: #region Body
263:
264: INFO = 0;
265: NB = this._ilaenv.Run(1, "DGEQRF", " ", M, N, - 1, - 1);
266: LWKOPT = N * NB;
267: WORK[1 + o_work] = LWKOPT;
268: LQUERY = (LWORK == - 1);
269: if (M < 0)
270: {
271: INFO = - 1;
272: }
273: else
274: {
275: if (N < 0)
276: {
277: INFO = - 2;
278: }
279: else
280: {
281: if (LDA < Math.Max(1, M))
282: {
283: INFO = - 4;
284: }
285: else
286: {
287: if (LWORK < Math.Max(1, N) && !LQUERY)
288: {
289: INFO = - 7;
290: }
291: }
292: }
293: }
294: if (INFO != 0)
295: {
296: this._xerbla.Run("DGEQRF", - INFO);
297: return;
298: }
299: else
300: {
301: if (LQUERY)
302: {
303: return;
304: }
305: }
306: // *
307: // * Quick return if possible
308: // *
309: K = Math.Min(M, N);
310: if (K == 0)
311: {
312: WORK[1 + o_work] = 1;
313: return;
314: }
315: // *
316: NBMIN = 2;
317: NX = 0;
318: IWS = N;
319: if (NB > 1 && NB < K)
320: {
321: // *
322: // * Determine when to cross over from blocked to unblocked code.
323: // *
324: NX = Math.Max(0, this._ilaenv.Run(3, "DGEQRF", " ", M, N, - 1, - 1));
325: if (NX < K)
326: {
327: // *
328: // * Determine if workspace is large enough for blocked code.
329: // *
330: LDWORK = N;
331: IWS = LDWORK * NB;
332: if (LWORK < IWS)
333: {
334: // *
335: // * Not enough workspace to use optimal NB: reduce NB and
336: // * determine the minimum value of NB.
337: // *
338: NB = LWORK / LDWORK;
339: NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGEQRF", " ", M, N, - 1, - 1));
340: }
341: }
342: }
343: // *
344: if (NB >= NBMIN && NB < K && NX < K)
345: {
346: // *
347: // * Use blocked code initially
348: // *
349: for (I = 1; (NB >= 0) ? (I <= K - NX) : (I >= K - NX); I += NB)
350: {
351: IB = Math.Min(K - I + 1, NB);
352: // *
353: // * Compute the QR factorization of the current block
354: // * A(i:m,i:i+ib-1)
355: // *
356: this._dgeqr2.Run(M - I + 1, IB, ref A, I+I * LDA + o_a, LDA, ref TAU, I + o_tau, ref WORK, offset_work
357: , ref IINFO);
358: if (I + IB <= N)
359: {
360: // *
361: // * Form the triangular factor of the block reflector
362: // * H = H(i) H(i+1) . . . H(i+ib-1)
363: // *
364: this._dlarft.Run("Forward", "Columnwise", M - I + 1, IB, ref A, I+I * LDA + o_a, LDA
365: , TAU, I + o_tau, ref WORK, offset_work, LDWORK);
366: // *
367: // * Apply H' to A(i:m,i+ib:n) from the left
368: // *
369: this._dlarfb.Run("Left", "Transpose", "Forward", "Columnwise", M - I + 1, N - I - IB + 1
370: , IB, A, I+I * LDA + o_a, LDA, WORK, offset_work, LDWORK, ref A, I+(I + IB) * LDA + o_a
371: , LDA, ref WORK, IB + 1 + o_work, LDWORK);
372: }
373: }
374: }
375: else
376: {
377: I = 1;
378: }
379: // *
380: // * Use unblocked code to factor the last or only block.
381: // *
382: if (I <= K)
383: {
384: this._dgeqr2.Run(M - I + 1, N - I + 1, ref A, I+I * LDA + o_a, LDA, ref TAU, I + o_tau, ref WORK, offset_work
385: , ref IINFO);
386: }
387: // *
388: WORK[1 + o_work] = IWS;
389: return;
390: // *
391: // * End of DGEQRF
392: // *
393:
394: #endregion
395:
396: }
397: }
398: }