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