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: /// DGEQP3 computes a QR factorization with column pivoting of a
27: /// matrix A: A*P = Q*R using Level 3 BLAS.
28: ///
29: ///</summary>
30: public class DGEQP3
31: {
32:
33:
34: #region Dependencies
35:
36: DGEQRF _dgeqrf; DLAQP2 _dlaqp2; DLAQPS _dlaqps; DORMQR _dormqr; DSWAP _dswap; XERBLA _xerbla; ILAENV _ilaenv;
37: DNRM2 _dnrm2;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const int INB = 1; const int INBMIN = 2; const int IXOVER = 3; bool LQUERY = false; int FJB = 0; int IWS = 0; int J = 0;
45: int JB = 0;int LWKOPT = 0; int MINMN = 0; int MINWS = 0; int NA = 0; int NB = 0; int NBMIN = 0; int NFXD = 0; int NX = 0;
46: int SM = 0;int SMINMN = 0; int SN = 0; int TOPBMN = 0;
47:
48: #endregion
49:
50: public DGEQP3(DGEQRF dgeqrf, DLAQP2 dlaqp2, DLAQPS dlaqps, DORMQR dormqr, DSWAP dswap, XERBLA xerbla, ILAENV ilaenv, DNRM2 dnrm2)
51: {
52:
53:
54: #region Set Dependencies
55:
56: this._dgeqrf = dgeqrf; this._dlaqp2 = dlaqp2; this._dlaqps = dlaqps; this._dormqr = dormqr; this._dswap = dswap;
57: this._xerbla = xerbla;this._ilaenv = ilaenv; this._dnrm2 = dnrm2;
58:
59: #endregion
60:
61: }
62:
63: public DGEQP3()
64: {
65:
66:
67: #region Dependencies (Initialization)
68:
69: LSAME lsame = new LSAME();
70: XERBLA xerbla = new XERBLA();
71: DLAMC3 dlamc3 = new DLAMC3();
72: DLAPY2 dlapy2 = new DLAPY2();
73: DNRM2 dnrm2 = new DNRM2();
74: DSCAL dscal = new DSCAL();
75: DCOPY dcopy = new DCOPY();
76: IEEECK ieeeck = new IEEECK();
77: IPARMQ iparmq = new IPARMQ();
78: DSWAP dswap = new DSWAP();
79: IDAMAX idamax = new IDAMAX();
80: DGEMV dgemv = new DGEMV(lsame, xerbla);
81: DGER dger = new DGER(xerbla);
82: DLARF dlarf = new DLARF(dgemv, dger, lsame);
83: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
84: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
85: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
86: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
87: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
88: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
89: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
90: DGEMM dgemm = new DGEMM(lsame, xerbla);
91: DTRMM dtrmm = new DTRMM(lsame, xerbla);
92: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
93: DTRMV dtrmv = new DTRMV(lsame, xerbla);
94: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
95: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
96: DGEQRF dgeqrf = new DGEQRF(dgeqr2, dlarfb, dlarft, xerbla, ilaenv);
97: DLAQP2 dlaqp2 = new DLAQP2(dlarf, dlarfg, dswap, idamax, dlamch, dnrm2);
98: DLAQPS dlaqps = new DLAQPS(dgemm, dgemv, dlarfg, dswap, idamax, dlamch, dnrm2);
99: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
100: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
101:
102: #endregion
103:
104:
105: #region Set Dependencies
106:
107: this._dgeqrf = dgeqrf; this._dlaqp2 = dlaqp2; this._dlaqps = dlaqps; this._dormqr = dormqr; this._dswap = dswap;
108: this._xerbla = xerbla;this._ilaenv = ilaenv; this._dnrm2 = dnrm2;
109:
110: #endregion
111:
112: }
113: /// <summary>
114: /// Purpose
115: /// =======
116: ///
117: /// DGEQP3 computes a QR factorization with column pivoting of a
118: /// matrix A: A*P = Q*R using Level 3 BLAS.
119: ///
120: ///</summary>
121: /// <param name="M">
122: /// (input) INTEGER
123: /// The number of rows of the matrix A. M .GE. 0.
124: ///</param>
125: /// <param name="N">
126: /// (input) INTEGER
127: /// The number of columns of the matrix A. N .GE. 0.
128: ///</param>
129: /// <param name="A">
130: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
131: /// On entry, the M-by-N matrix A.
132: /// On exit, the upper triangle of the array contains the
133: /// min(M,N)-by-N upper trapezoidal matrix R; the elements below
134: /// the diagonal, together with the array TAU, represent the
135: /// orthogonal matrix Q as a product of min(M,N) elementary
136: /// reflectors.
137: ///</param>
138: /// <param name="LDA">
139: /// (input) INTEGER
140: /// The leading dimension of the array A. LDA .GE. max(1,M).
141: ///</param>
142: /// <param name="JPVT">
143: /// (input/output) INTEGER array, dimension (N)
144: /// On entry, if JPVT(J).ne.0, the J-th column of A is permuted
145: /// to the front of A*P (a leading column); if JPVT(J)=0,
146: /// the J-th column of A is a free column.
147: /// On exit, if JPVT(J)=K, then the J-th column of A*P was the
148: /// the K-th column of A.
149: ///</param>
150: /// <param name="TAU">
151: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
152: /// The scalar factors of the elementary reflectors.
153: ///</param>
154: /// <param name="WORK">
155: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
156: /// On exit, if INFO=0, WORK(1) returns the optimal LWORK.
157: ///</param>
158: /// <param name="LWORK">
159: /// (input) INTEGER
160: /// The dimension of the array WORK. LWORK .GE. 3*N+1.
161: /// For optimal performance LWORK .GE. 2*N+( N+1 )*NB, where NB
162: /// is the optimal blocksize.
163: ///
164: /// If LWORK = -1, then a workspace query is assumed; the routine
165: /// only calculates the optimal size of the WORK array, returns
166: /// this value as the first entry of the WORK array, and no error
167: /// message related to LWORK is issued by XERBLA.
168: ///</param>
169: /// <param name="INFO">
170: /// (output) INTEGER
171: /// = 0: successful exit.
172: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
173: ///</param>
174: public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref int[] JPVT, int offset_jpvt, ref double[] TAU, int offset_tau
175: , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
176: {
177:
178: #region Array Index Correction
179:
180: int o_a = -1 - LDA + offset_a; int o_jpvt = -1 + offset_jpvt; int o_tau = -1 + offset_tau;
181: int o_work = -1 + offset_work;
182:
183: #endregion
184:
185:
186: #region Prolog
187:
188: // *
189: // * -- LAPACK routine (version 3.1) --
190: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
191: // * November 2006
192: // *
193: // * .. Scalar Arguments ..
194: // * ..
195: // * .. Array Arguments ..
196: // * ..
197: // *
198: // * Purpose
199: // * =======
200: // *
201: // * DGEQP3 computes a QR factorization with column pivoting of a
202: // * matrix A: A*P = Q*R using Level 3 BLAS.
203: // *
204: // * Arguments
205: // * =========
206: // *
207: // * M (input) INTEGER
208: // * The number of rows of the matrix A. M >= 0.
209: // *
210: // * N (input) INTEGER
211: // * The number of columns of the matrix A. N >= 0.
212: // *
213: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
214: // * On entry, the M-by-N matrix A.
215: // * On exit, the upper triangle of the array contains the
216: // * min(M,N)-by-N upper trapezoidal matrix R; the elements below
217: // * the diagonal, together with the array TAU, represent the
218: // * orthogonal matrix Q as a product of min(M,N) elementary
219: // * reflectors.
220: // *
221: // * LDA (input) INTEGER
222: // * The leading dimension of the array A. LDA >= max(1,M).
223: // *
224: // * JPVT (input/output) INTEGER array, dimension (N)
225: // * On entry, if JPVT(J).ne.0, the J-th column of A is permuted
226: // * to the front of A*P (a leading column); if JPVT(J)=0,
227: // * the J-th column of A is a free column.
228: // * On exit, if JPVT(J)=K, then the J-th column of A*P was the
229: // * the K-th column of A.
230: // *
231: // * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
232: // * The scalar factors of the elementary reflectors.
233: // *
234: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
235: // * On exit, if INFO=0, WORK(1) returns the optimal LWORK.
236: // *
237: // * LWORK (input) INTEGER
238: // * The dimension of the array WORK. LWORK >= 3*N+1.
239: // * For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB
240: // * is the optimal blocksize.
241: // *
242: // * If LWORK = -1, then a workspace query is assumed; the routine
243: // * only calculates the optimal size of the WORK array, returns
244: // * this value as the first entry of the WORK array, and no error
245: // * message related to LWORK is issued by XERBLA.
246: // *
247: // * INFO (output) INTEGER
248: // * = 0: successful exit.
249: // * < 0: if INFO = -i, the i-th argument had an illegal value.
250: // *
251: // * Further Details
252: // * ===============
253: // *
254: // * The matrix Q is represented as a product of elementary reflectors
255: // *
256: // * Q = H(1) H(2) . . . H(k), where k = min(m,n).
257: // *
258: // * Each H(i) has the form
259: // *
260: // * H(i) = I - tau * v * v'
261: // *
262: // * where tau is a real/complex scalar, and v is a real/complex vector
263: // * with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
264: // * A(i+1:m,i), and tau in TAU(i).
265: // *
266: // * Based on contributions by
267: // * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
268: // * X. Sun, Computer Science Dept., Duke University, USA
269: // *
270: // * =====================================================================
271: // *
272: // * .. Parameters ..
273: // * ..
274: // * .. Local Scalars ..
275: // * ..
276: // * .. External Subroutines ..
277: // * ..
278: // * .. External Functions ..
279: // * ..
280: // * .. Intrinsic Functions ..
281: // INTRINSIC INT, MAX, MIN;
282: // * ..
283: // * .. Executable Statements ..
284: // *
285: // * Test input arguments
286: // * ====================
287: // *
288:
289: #endregion
290:
291:
292: #region Body
293:
294: INFO = 0;
295: LQUERY = (LWORK == - 1);
296: if (M < 0)
297: {
298: INFO = - 1;
299: }
300: else
301: {
302: if (N < 0)
303: {
304: INFO = - 2;
305: }
306: else
307: {
308: if (LDA < Math.Max(1, M))
309: {
310: INFO = - 4;
311: }
312: }
313: }
314: // *
315: if (INFO == 0)
316: {
317: MINMN = Math.Min(M, N);
318: if (MINMN == 0)
319: {
320: IWS = 1;
321: LWKOPT = 1;
322: }
323: else
324: {
325: IWS = 3 * N + 1;
326: NB = this._ilaenv.Run(INB, "DGEQRF", " ", M, N, - 1, - 1);
327: LWKOPT = 2 * N + (N + 1) * NB;
328: }
329: WORK[1 + o_work] = LWKOPT;
330: // *
331: if ((LWORK < IWS) && !LQUERY)
332: {
333: INFO = - 8;
334: }
335: }
336: // *
337: if (INFO != 0)
338: {
339: this._xerbla.Run("DGEQP3", - INFO);
340: return;
341: }
342: else
343: {
344: if (LQUERY)
345: {
346: return;
347: }
348: }
349: // *
350: // * Quick return if possible.
351: // *
352: if (MINMN == 0)
353: {
354: return;
355: }
356: // *
357: // * Move initial columns up front.
358: // *
359: NFXD = 1;
360: for (J = 1; J <= N; J++)
361: {
362: if (JPVT[J + o_jpvt] != 0)
363: {
364: if (J != NFXD)
365: {
366: this._dswap.Run(M, ref A, 1+J * LDA + o_a, 1, ref A, 1+NFXD * LDA + o_a, 1);
367: JPVT[J + o_jpvt] = JPVT[NFXD + o_jpvt];
368: JPVT[NFXD + o_jpvt] = J;
369: }
370: else
371: {
372: JPVT[J + o_jpvt] = J;
373: }
374: NFXD = NFXD + 1;
375: }
376: else
377: {
378: JPVT[J + o_jpvt] = J;
379: }
380: }
381: NFXD = NFXD - 1;
382: // *
383: // * Factorize fixed columns
384: // * =======================
385: // *
386: // * Compute the QR factorization of fixed columns and update
387: // * remaining columns.
388: // *
389: if (NFXD > 0)
390: {
391: NA = Math.Min(M, NFXD);
392: // *CC CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
393: this._dgeqrf.Run(M, NA, ref A, offset_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
394: , LWORK, ref INFO);
395: IWS = Math.Max(IWS, Convert.ToInt32(Math.Truncate(WORK[1 + o_work])));
396: if (NA < N)
397: {
398: // *CC CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
399: // *CC $ TAU, A( 1, NA+1 ), LDA, WORK, INFO )
400: this._dormqr.Run("Left", "Transpose", M, N - NA, NA, ref A, offset_a
401: , LDA, TAU, offset_tau, ref A, 1+(NA + 1) * LDA + o_a, LDA, ref WORK, offset_work, LWORK
402: , ref INFO);
403: IWS = Math.Max(IWS, Convert.ToInt32(Math.Truncate(WORK[1 + o_work])));
404: }
405: }
406: // *
407: // * Factorize free columns
408: // * ======================
409: // *
410: if (NFXD < MINMN)
411: {
412: // *
413: SM = M - NFXD;
414: SN = N - NFXD;
415: SMINMN = MINMN - NFXD;
416: // *
417: // * Determine the block size.
418: // *
419: NB = this._ilaenv.Run(INB, "DGEQRF", " ", SM, SN, - 1, - 1);
420: NBMIN = 2;
421: NX = 0;
422: // *
423: if ((NB > 1) && (NB < SMINMN))
424: {
425: // *
426: // * Determine when to cross over from blocked to unblocked code.
427: // *
428: NX = Math.Max(0, this._ilaenv.Run(IXOVER, "DGEQRF", " ", SM, SN, - 1, - 1));
429: // *
430: // *
431: if (NX < SMINMN)
432: {
433: // *
434: // * Determine if workspace is large enough for blocked code.
435: // *
436: MINWS = 2 * SN + (SN + 1) * NB;
437: IWS = Math.Max(IWS, MINWS);
438: if (LWORK < MINWS)
439: {
440: // *
441: // * Not enough workspace to use optimal NB: Reduce NB and
442: // * determine the minimum value of NB.
443: // *
444: NB = (LWORK - 2 * SN) / (SN + 1);
445: NBMIN = Math.Max(2, this._ilaenv.Run(INBMIN, "DGEQRF", " ", SM, SN, - 1, - 1));
446: // *
447: // *
448: }
449: }
450: }
451: // *
452: // * Initialize partial column norms. The first N elements of work
453: // * store the exact column norms.
454: // *
455: for (J = NFXD + 1; J <= N; J++)
456: {
457: WORK[J + o_work] = this._dnrm2.Run(SM, A, NFXD + 1+J * LDA + o_a, 1);
458: WORK[N + J + o_work] = WORK[J + o_work];
459: }
460: // *
461: if ((NB >= NBMIN) && (NB < SMINMN) && (NX < SMINMN))
462: {
463: // *
464: // * Use blocked code initially.
465: // *
466: J = NFXD + 1;
467: // *
468: // * Compute factorization: while loop.
469: // *
470: // *
471: TOPBMN = MINMN - NX;
472: LABEL30:;
473: if (J <= TOPBMN)
474: {
475: JB = Math.Min(NB, TOPBMN - J + 1);
476: // *
477: // * Factorize JB columns among columns J:N.
478: // *
479: this._dlaqps.Run(M, N - J + 1, J - 1, JB, ref FJB, ref A, 1+J * LDA + o_a
480: , LDA, ref JPVT, J + o_jpvt, ref TAU, J + o_tau, ref WORK, J + o_work, ref WORK, N + J + o_work, ref WORK, 2 * N + 1 + o_work
481: , ref WORK, 2 * N + JB + 1 + o_work, N - J + 1);
482: // *
483: J = J + FJB;
484: goto LABEL30;
485: }
486: }
487: else
488: {
489: J = NFXD + 1;
490: }
491: // *
492: // * Use unblocked code to factor the last or only block.
493: // *
494: // *
495: if (J <= MINMN)
496: {
497: this._dlaqp2.Run(M, N - J + 1, J - 1, ref A, 1+J * LDA + o_a, LDA, ref JPVT, J + o_jpvt
498: , ref TAU, J + o_tau, ref WORK, J + o_work, ref WORK, N + J + o_work, ref WORK, 2 * N + 1 + o_work);
499: }
500: // *
501: }
502: // *
503: WORK[1 + o_work] = IWS;
504: return;
505: // *
506: // * End of DGEQP3
507: // *
508:
509: #endregion
510:
511: }
512: }
513: }