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 auxiliary routine (version 3.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// November 2006
23: /// Purpose
24: /// =======
25: ///
26: /// DLAQP2 computes a QR factorization with column pivoting of
27: /// the block A(OFFSET+1:M,1:N).
28: /// The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
29: ///
30: ///</summary>
31: public class DLAQP2
32: {
33:
34:
35: #region Dependencies
36:
37: DLARF _dlarf; DLARFG _dlarfg; DSWAP _dswap; IDAMAX _idamax; DLAMCH _dlamch; DNRM2 _dnrm2;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0; int ITEMP = 0; int J = 0; int MN = 0; int OFFPI = 0;
45: int PVT = 0;double AII = 0; double TEMP = 0; double TEMP2 = 0; double TOL3Z = 0;
46:
47: #endregion
48:
49: public DLAQP2(DLARF dlarf, DLARFG dlarfg, DSWAP dswap, IDAMAX idamax, DLAMCH dlamch, DNRM2 dnrm2)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._dlarf = dlarf; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax; this._dlamch = dlamch;
56: this._dnrm2 = dnrm2;
57:
58: #endregion
59:
60: }
61:
62: public DLAQP2()
63: {
64:
65:
66: #region Dependencies (Initialization)
67:
68: LSAME lsame = new LSAME();
69: XERBLA xerbla = new XERBLA();
70: DLAMC3 dlamc3 = new DLAMC3();
71: DLAPY2 dlapy2 = new DLAPY2();
72: DNRM2 dnrm2 = new DNRM2();
73: DSCAL dscal = new DSCAL();
74: DSWAP dswap = new DSWAP();
75: IDAMAX idamax = new IDAMAX();
76: DGEMV dgemv = new DGEMV(lsame, xerbla);
77: DGER dger = new DGER(xerbla);
78: DLARF dlarf = new DLARF(dgemv, dger, lsame);
79: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
80: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
81: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
82: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
83: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
84: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
85:
86: #endregion
87:
88:
89: #region Set Dependencies
90:
91: this._dlarf = dlarf; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax; this._dlamch = dlamch;
92: this._dnrm2 = dnrm2;
93:
94: #endregion
95:
96: }
97: /// <summary>
98: /// Purpose
99: /// =======
100: ///
101: /// DLAQP2 computes a QR factorization with column pivoting of
102: /// the block A(OFFSET+1:M,1:N).
103: /// The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
104: ///
105: ///</summary>
106: /// <param name="M">
107: /// (input) INTEGER
108: /// The number of rows of the matrix A. M .GE. 0.
109: ///</param>
110: /// <param name="N">
111: /// (input) INTEGER
112: /// The number of columns of the matrix A. N .GE. 0.
113: ///</param>
114: /// <param name="OFFSET">
115: /// (input) INTEGER
116: /// The number of rows of the matrix A that must be pivoted
117: /// but no factorized. OFFSET .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 block A(OFFSET+1:M,1:N) is
123: /// the triangular factor obtained; the elements in block
124: /// A(OFFSET+1:M,1:N) below the diagonal, together with the
125: /// array TAU, represent the orthogonal matrix Q as a product of
126: /// elementary reflectors. Block A(1:OFFSET,1:N) has been
127: /// accordingly pivoted, but no factorized.
128: ///</param>
129: /// <param name="LDA">
130: /// (input) INTEGER
131: /// The leading dimension of the array A. LDA .GE. max(1,M).
132: ///</param>
133: /// <param name="JPVT">
134: /// (input/output) INTEGER array, dimension (N)
135: /// On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
136: /// to the front of A*P (a leading column); if JPVT(i) = 0,
137: /// the i-th column of A is a free column.
138: /// On exit, if JPVT(i) = k, then the i-th column of A*P
139: /// was the k-th column of A.
140: ///</param>
141: /// <param name="TAU">
142: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
143: /// The scalar factors of the elementary reflectors.
144: ///</param>
145: /// <param name="VN1">
146: /// (input/output) DOUBLE PRECISION array, dimension (N)
147: /// The vector with the partial column norms.
148: ///</param>
149: /// <param name="VN2">
150: /// (input/output) DOUBLE PRECISION array, dimension (N)
151: /// The vector with the exact column norms.
152: ///</param>
153: /// <param name="WORK">
154: /// (workspace) DOUBLE PRECISION array, dimension (N)
155: ///</param>
156: public void Run(int M, int N, int OFFSET, ref double[] A, int offset_a, int LDA, ref int[] JPVT, int offset_jpvt
157: , ref double[] TAU, int offset_tau, ref double[] VN1, int offset_vn1, ref double[] VN2, int offset_vn2, ref double[] WORK, int offset_work)
158: {
159:
160: #region Array Index Correction
161:
162: int o_a = -1 - LDA + offset_a; int o_jpvt = -1 + offset_jpvt; int o_tau = -1 + offset_tau;
163: int o_vn1 = -1 + offset_vn1; int o_vn2 = -1 + offset_vn2; int o_work = -1 + offset_work;
164:
165: #endregion
166:
167:
168: #region Prolog
169:
170: // *
171: // * -- LAPACK auxiliary 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: // * DLAQP2 computes a QR factorization with column pivoting of
184: // * the block A(OFFSET+1:M,1:N).
185: // * The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
186: // *
187: // * Arguments
188: // * =========
189: // *
190: // * M (input) INTEGER
191: // * The number of rows of the matrix A. M >= 0.
192: // *
193: // * N (input) INTEGER
194: // * The number of columns of the matrix A. N >= 0.
195: // *
196: // * OFFSET (input) INTEGER
197: // * The number of rows of the matrix A that must be pivoted
198: // * but no factorized. OFFSET >= 0.
199: // *
200: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
201: // * On entry, the M-by-N matrix A.
202: // * On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
203: // * the triangular factor obtained; the elements in block
204: // * A(OFFSET+1:M,1:N) below the diagonal, together with the
205: // * array TAU, represent the orthogonal matrix Q as a product of
206: // * elementary reflectors. Block A(1:OFFSET,1:N) has been
207: // * accordingly pivoted, but no factorized.
208: // *
209: // * LDA (input) INTEGER
210: // * The leading dimension of the array A. LDA >= max(1,M).
211: // *
212: // * JPVT (input/output) INTEGER array, dimension (N)
213: // * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
214: // * to the front of A*P (a leading column); if JPVT(i) = 0,
215: // * the i-th column of A is a free column.
216: // * On exit, if JPVT(i) = k, then the i-th column of A*P
217: // * was the k-th column of A.
218: // *
219: // * TAU (output) DOUBLE PRECISION array, dimension (min(M,N))
220: // * The scalar factors of the elementary reflectors.
221: // *
222: // * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
223: // * The vector with the partial column norms.
224: // *
225: // * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
226: // * The vector with the exact column norms.
227: // *
228: // * WORK (workspace) DOUBLE PRECISION array, dimension (N)
229: // *
230: // * Further Details
231: // * ===============
232: // *
233: // * Based on contributions by
234: // * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
235: // * X. Sun, Computer Science Dept., Duke University, USA
236: // *
237: // * Partial column norm updating strategy modified by
238: // * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
239: // * University of Zagreb, Croatia.
240: // * June 2006.
241: // * For more details see LAPACK Working Note 176.
242: // * =====================================================================
243: // *
244: // * .. Parameters ..
245: // * ..
246: // * .. Local Scalars ..
247: // * ..
248: // * .. External Subroutines ..
249: // * ..
250: // * .. Intrinsic Functions ..
251: // INTRINSIC ABS, MAX, MIN, SQRT;
252: // * ..
253: // * .. External Functions ..
254: // * ..
255: // * .. Executable Statements ..
256: // *
257:
258: #endregion
259:
260:
261: #region Body
262:
263: MN = Math.Min(M - OFFSET, N);
264: TOL3Z = Math.Sqrt(this._dlamch.Run("Epsilon"));
265: // *
266: // * Compute factorization.
267: // *
268: for (I = 1; I <= MN; I++)
269: {
270: // *
271: OFFPI = OFFSET + I;
272: // *
273: // * Determine ith pivot column and swap if necessary.
274: // *
275: PVT = (I - 1) + this._idamax.Run(N - I + 1, VN1, I + o_vn1, 1);
276: // *
277: if (PVT != I)
278: {
279: this._dswap.Run(M, ref A, 1+PVT * LDA + o_a, 1, ref A, 1+I * LDA + o_a, 1);
280: ITEMP = JPVT[PVT + o_jpvt];
281: JPVT[PVT + o_jpvt] = JPVT[I + o_jpvt];
282: JPVT[I + o_jpvt] = ITEMP;
283: VN1[PVT + o_vn1] = VN1[I + o_vn1];
284: VN2[PVT + o_vn2] = VN2[I + o_vn2];
285: }
286: // *
287: // * Generate elementary reflector H(i).
288: // *
289: if (OFFPI < M)
290: {
291: this._dlarfg.Run(M - OFFPI + 1, ref A[OFFPI+I * LDA + o_a], ref A, OFFPI + 1+I * LDA + o_a, 1, ref TAU[I + o_tau]);
292: }
293: else
294: {
295: this._dlarfg.Run(1, ref A[M+I * LDA + o_a], ref A, M+I * LDA + o_a, 1, ref TAU[I + o_tau]);
296: }
297: // *
298: if (I < N)
299: {
300: // *
301: // * Apply H(i)' to A(offset+i:m,i+1:n) from the left.
302: // *
303: AII = A[OFFPI+I * LDA + o_a];
304: A[OFFPI+I * LDA + o_a] = ONE;
305: this._dlarf.Run("Left", M - OFFPI + 1, N - I, A, OFFPI+I * LDA + o_a, 1, TAU[I + o_tau]
306: , ref A, OFFPI+(I + 1) * LDA + o_a, LDA, ref WORK, 1 + o_work);
307: A[OFFPI+I * LDA + o_a] = AII;
308: }
309: // *
310: // * Update partial column norms.
311: // *
312: for (J = I + 1; J <= N; J++)
313: {
314: if (VN1[J + o_vn1] != ZERO)
315: {
316: // *
317: // * NOTE: The following 4 lines follow from the analysis in
318: // * Lapack Working Note 176.
319: // *
320: TEMP = ONE - Math.Pow(Math.Abs(A[OFFPI+J * LDA + o_a]) / VN1[J + o_vn1],2);
321: TEMP = Math.Max(TEMP, ZERO);
322: TEMP2 = TEMP * Math.Pow(VN1[J + o_vn1] / VN2[J + o_vn2],2);
323: if (TEMP2 <= TOL3Z)
324: {
325: if (OFFPI < M)
326: {
327: VN1[J + o_vn1] = this._dnrm2.Run(M - OFFPI, A, OFFPI + 1+J * LDA + o_a, 1);
328: VN2[J + o_vn2] = VN1[J + o_vn1];
329: }
330: else
331: {
332: VN1[J + o_vn1] = ZERO;
333: VN2[J + o_vn2] = ZERO;
334: }
335: }
336: else
337: {
338: VN1[J + o_vn1] = VN1[J + o_vn1] * Math.Sqrt(TEMP);
339: }
340: }
341: }
342: // *
343: }
344: // *
345: return;
346: // *
347: // * End of DLAQP2
348: // *
349:
350: #endregion
351:
352: }
353: }
354: }