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 deprecated driver routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// This routine is deprecated and has been replaced by routine DGEQP3.
27: ///
28: /// DGEQPF computes a QR factorization with column pivoting of a
29: /// real M-by-N matrix A: A*P = Q*R.
30: ///
31: ///</summary>
32: public class DGEQPF
33: {
34:
35:
36: #region Dependencies
37:
38: DGEQR2 _dgeqr2; DLARF _dlarf; DLARFG _dlarfg; DORM2R _dorm2r; DSWAP _dswap; XERBLA _xerbla; IDAMAX _idamax;
39: DLAMCH _dlamch;DNRM2 _dnrm2;
40:
41: #endregion
42:
43:
44: #region Fields
45:
46: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; int ITEMP = 0; int J = 0; int MA = 0; int MN = 0;
47: int PVT = 0;double AII = 0; double TEMP = 0; double TEMP2 = 0; double TOL3Z = 0;
48:
49: #endregion
50:
51: public DGEQPF(DGEQR2 dgeqr2, DLARF dlarf, DLARFG dlarfg, DORM2R dorm2r, DSWAP dswap, XERBLA xerbla, IDAMAX idamax, DLAMCH dlamch, DNRM2 dnrm2)
52: {
53:
54:
55: #region Set Dependencies
56:
57: this._dgeqr2 = dgeqr2; this._dlarf = dlarf; this._dlarfg = dlarfg; this._dorm2r = dorm2r; this._dswap = dswap;
58: this._xerbla = xerbla;this._idamax = idamax; this._dlamch = dlamch; this._dnrm2 = dnrm2;
59:
60: #endregion
61:
62: }
63:
64: public DGEQPF()
65: {
66:
67:
68: #region Dependencies (Initialization)
69:
70: LSAME lsame = new LSAME();
71: XERBLA xerbla = new XERBLA();
72: DLAMC3 dlamc3 = new DLAMC3();
73: DLAPY2 dlapy2 = new DLAPY2();
74: DNRM2 dnrm2 = new DNRM2();
75: DSCAL dscal = new DSCAL();
76: DSWAP dswap = new DSWAP();
77: IDAMAX idamax = new IDAMAX();
78: DGEMV dgemv = new DGEMV(lsame, xerbla);
79: DGER dger = new DGER(xerbla);
80: DLARF dlarf = new DLARF(dgemv, dger, lsame);
81: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
82: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
83: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
84: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
85: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
86: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
87: DGEQR2 dgeqr2 = new DGEQR2(dlarf, dlarfg, xerbla);
88: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
89:
90: #endregion
91:
92:
93: #region Set Dependencies
94:
95: this._dgeqr2 = dgeqr2; this._dlarf = dlarf; this._dlarfg = dlarfg; this._dorm2r = dorm2r; this._dswap = dswap;
96: this._xerbla = xerbla;this._idamax = idamax; this._dlamch = dlamch; this._dnrm2 = dnrm2;
97:
98: #endregion
99:
100: }
101: /// <summary>
102: /// Purpose
103: /// =======
104: ///
105: /// This routine is deprecated and has been replaced by routine DGEQP3.
106: ///
107: /// DGEQPF computes a QR factorization with column pivoting of a
108: /// real M-by-N matrix A: A*P = Q*R.
109: ///
110: ///</summary>
111: /// <param name="M">
112: /// (input) INTEGER
113: /// The number of rows of the matrix A. M .GE. 0.
114: ///</param>
115: /// <param name="N">
116: /// (input) INTEGER
117: /// The number of columns of the matrix A. N .GE. 0
118: ///</param>
119: /// <param name="A">
120: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
121: /// On entry, the M-by-N matrix A.
122: /// On exit, the upper triangle of the array contains the
123: /// min(M,N)-by-N upper triangular matrix R; the elements
124: /// below the diagonal, together with the array TAU,
125: /// represent the orthogonal matrix Q as a product of
126: /// min(m,n) elementary reflectors.
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="JPVT">
133: /// (input/output) INTEGER array, dimension (N)
134: /// On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
135: /// to the front of A*P (a leading column); if JPVT(i) = 0,
136: /// the i-th column of A is a free column.
137: /// On exit, if JPVT(i) = k, then the i-th column of A*P
138: /// was the k-th column of A.
139: ///</param>
140: /// <param name="TAU">
141: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
142: /// The scalar factors of the elementary reflectors.
143: ///</param>
144: /// <param name="WORK">
145: /// (workspace) DOUBLE PRECISION array, dimension (3*N)
146: ///</param>
147: /// <param name="INFO">
148: /// (output) INTEGER
149: /// = 0: successful exit
150: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
151: ///</param>
152: 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
153: , ref double[] WORK, int offset_work, ref int INFO)
154: {
155:
156: #region Array Index Correction
157:
158: int o_a = -1 - LDA + offset_a; int o_jpvt = -1 + offset_jpvt; int o_tau = -1 + offset_tau;
159: int o_work = -1 + offset_work;
160:
161: #endregion
162:
163:
164: #region Prolog
165:
166: // *
167: // * -- LAPACK deprecated driver routine (version 3.1) --
168: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
169: // * November 2006
170: // *
171: // * .. Scalar Arguments ..
172: // * ..
173: // * .. Array Arguments ..
174: // * ..
175: // *
176: // * Purpose
177: // * =======
178: // *
179: // * This routine is deprecated and has been replaced by routine DGEQP3.
180: // *
181: // * DGEQPF computes a QR factorization with column pivoting of a
182: // * real M-by-N matrix A: A*P = Q*R.
183: // *
184: // * Arguments
185: // * =========
186: // *
187: // * M (input) INTEGER
188: // * The number of rows of the matrix A. M >= 0.
189: // *
190: // * N (input) INTEGER
191: // * The number of columns of the matrix A. N >= 0
192: // *
193: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
194: // * On entry, the M-by-N matrix A.
195: // * On exit, the upper triangle of the array contains the
196: // * min(M,N)-by-N upper triangular matrix R; the elements
197: // * below the diagonal, together with the array TAU,
198: // * represent the orthogonal matrix Q as a product of
199: // * min(m,n) elementary reflectors.
200: // *
201: // * LDA (input) INTEGER
202: // * The leading dimension of the array A. LDA >= max(1,M).
203: // *
204: // * JPVT (input/output) INTEGER array, dimension (N)
205: // * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
206: // * to the front of A*P (a leading column); if JPVT(i) = 0,
207: // * the i-th column of A is a free column.
208: // * On exit, if JPVT(i) = k, then the i-th column of A*P
209: // * was the k-th column of A.
210: // *
211: // * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
212: // * The scalar factors of the elementary reflectors.
213: // *
214: // * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
215: // *
216: // * INFO (output) INTEGER
217: // * = 0: successful exit
218: // * < 0: if INFO = -i, the i-th argument had an illegal value
219: // *
220: // * Further Details
221: // * ===============
222: // *
223: // * The matrix Q is represented as a product of elementary reflectors
224: // *
225: // * Q = H(1) H(2) . . . H(n)
226: // *
227: // * Each H(i) has the form
228: // *
229: // * H = I - tau * v * v'
230: // *
231: // * where tau is a real scalar, and v is a real vector with
232: // * v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
233: // *
234: // * The matrix P is represented in jpvt as follows: If
235: // * jpvt(j) = i
236: // * then the jth column of P is the ith canonical unit vector.
237: // *
238: // * Partial column norm updating strategy modified by
239: // * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
240: // * University of Zagreb, Croatia.
241: // * June 2006.
242: // * For more details see LAPACK Working Note 176.
243: // *
244: // * =====================================================================
245: // *
246: // * .. Parameters ..
247: // * ..
248: // * .. Local Scalars ..
249: // * ..
250: // * .. External Subroutines ..
251: // * ..
252: // * .. Intrinsic Functions ..
253: // INTRINSIC ABS, MAX, MIN, SQRT;
254: // * ..
255: // * .. External Functions ..
256: // * ..
257: // * .. Executable Statements ..
258: // *
259: // * Test the input arguments
260: // *
261:
262: #endregion
263:
264:
265: #region Body
266:
267: INFO = 0;
268: if (M < 0)
269: {
270: INFO = - 1;
271: }
272: else
273: {
274: if (N < 0)
275: {
276: INFO = - 2;
277: }
278: else
279: {
280: if (LDA < Math.Max(1, M))
281: {
282: INFO = - 4;
283: }
284: }
285: }
286: if (INFO != 0)
287: {
288: this._xerbla.Run("DGEQPF", - INFO);
289: return;
290: }
291: // *
292: MN = Math.Min(M, N);
293: TOL3Z = Math.Sqrt(this._dlamch.Run("Epsilon"));
294: // *
295: // * Move initial columns up front
296: // *
297: ITEMP = 1;
298: for (I = 1; I <= N; I++)
299: {
300: if (JPVT[I + o_jpvt] != 0)
301: {
302: if (I != ITEMP)
303: {
304: this._dswap.Run(M, ref A, 1+I * LDA + o_a, 1, ref A, 1+ITEMP * LDA + o_a, 1);
305: JPVT[I + o_jpvt] = JPVT[ITEMP + o_jpvt];
306: JPVT[ITEMP + o_jpvt] = I;
307: }
308: else
309: {
310: JPVT[I + o_jpvt] = I;
311: }
312: ITEMP = ITEMP + 1;
313: }
314: else
315: {
316: JPVT[I + o_jpvt] = I;
317: }
318: }
319: ITEMP = ITEMP - 1;
320: // *
321: // * Compute the QR factorization and update remaining columns
322: // *
323: if (ITEMP > 0)
324: {
325: MA = Math.Min(ITEMP, M);
326: this._dgeqr2.Run(M, MA, ref A, offset_a, LDA, ref TAU, offset_tau, ref WORK, offset_work
327: , ref INFO);
328: if (MA < N)
329: {
330: this._dorm2r.Run("Left", "Transpose", M, N - MA, MA, ref A, offset_a
331: , LDA, TAU, offset_tau, ref A, 1+(MA + 1) * LDA + o_a, LDA, ref WORK, offset_work, ref INFO);
332: }
333: }
334: // *
335: if (ITEMP < MN)
336: {
337: // *
338: // * Initialize partial column norms. The first n elements of
339: // * work store the exact column norms.
340: // *
341: for (I = ITEMP + 1; I <= N; I++)
342: {
343: WORK[I + o_work] = this._dnrm2.Run(M - ITEMP, A, ITEMP + 1+I * LDA + o_a, 1);
344: WORK[N + I + o_work] = WORK[I + o_work];
345: }
346: // *
347: // * Compute factorization
348: // *
349: for (I = ITEMP + 1; I <= MN; I++)
350: {
351: // *
352: // * Determine ith pivot column and swap if necessary
353: // *
354: PVT = (I - 1) + this._idamax.Run(N - I + 1, WORK, I + o_work, 1);
355: // *
356: if (PVT != I)
357: {
358: this._dswap.Run(M, ref A, 1+PVT * LDA + o_a, 1, ref A, 1+I * LDA + o_a, 1);
359: ITEMP = JPVT[PVT + o_jpvt];
360: JPVT[PVT + o_jpvt] = JPVT[I + o_jpvt];
361: JPVT[I + o_jpvt] = ITEMP;
362: WORK[PVT + o_work] = WORK[I + o_work];
363: WORK[N + PVT + o_work] = WORK[N + I + o_work];
364: }
365: // *
366: // * Generate elementary reflector H(i)
367: // *
368: if (I < M)
369: {
370: this._dlarfg.Run(M - I + 1, ref A[I+I * LDA + o_a], ref A, I + 1+I * LDA + o_a, 1, ref TAU[I + o_tau]);
371: }
372: else
373: {
374: this._dlarfg.Run(1, ref A[M+M * LDA + o_a], ref A, M+M * LDA + o_a, 1, ref TAU[M + o_tau]);
375: }
376: // *
377: if (I < N)
378: {
379: // *
380: // * Apply H(i) to A(i:m,i+1:n) from the left
381: // *
382: AII = A[I+I * LDA + o_a];
383: A[I+I * LDA + o_a] = ONE;
384: this._dlarf.Run("LEFT", M - I + 1, N - I, A, I+I * LDA + o_a, 1, TAU[I + o_tau]
385: , ref A, I+(I + 1) * LDA + o_a, LDA, ref WORK, 2 * N + 1 + o_work);
386: A[I+I * LDA + o_a] = AII;
387: }
388: // *
389: // * Update partial column norms
390: // *
391: for (J = I + 1; J <= N; J++)
392: {
393: if (WORK[J + o_work] != ZERO)
394: {
395: // *
396: // * NOTE: The following 4 lines follow from the analysis in
397: // * Lapack Working Note 176.
398: // *
399: TEMP = Math.Abs(A[I+J * LDA + o_a]) / WORK[J + o_work];
400: TEMP = Math.Max(ZERO, (ONE + TEMP) * (ONE - TEMP));
401: TEMP2 = TEMP * Math.Pow(WORK[J + o_work] / WORK[N + J + o_work],2);
402: if (TEMP2 <= TOL3Z)
403: {
404: if (M - I > 0)
405: {
406: WORK[J + o_work] = this._dnrm2.Run(M - I, A, I + 1+J * LDA + o_a, 1);
407: WORK[N + J + o_work] = WORK[J + o_work];
408: }
409: else
410: {
411: WORK[J + o_work] = ZERO;
412: WORK[N + J + o_work] = ZERO;
413: }
414: }
415: else
416: {
417: WORK[J + o_work] = WORK[J + o_work] * Math.Sqrt(TEMP);
418: }
419: }
420: }
421: // *
422: }
423: }
424: return;
425: // *
426: // * End of DGEQPF
427: // *
428:
429: #endregion
430:
431: }
432: }
433: }