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