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: /// DLAQPS computes a step of QR factorization with column pivoting
27: /// of a real M-by-N matrix A by using Blas-3. It tries to factorize
28: /// NB columns from A starting from the row OFFSET+1, and updates all
29: /// of the matrix with Blas-3 xGEMM.
30: ///
31: /// In some cases, due to catastrophic cancellations, it cannot
32: /// factorize NB columns. Hence, the actual number of factorized
33: /// columns is returned in KB.
34: ///
35: /// Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
36: ///
37: ///</summary>
38: public class DLAQPS
39: {
40:
41:
42: #region Dependencies
43:
44: DGEMM _dgemm; DGEMV _dgemv; DLARFG _dlarfg; DSWAP _dswap; IDAMAX _idamax; DLAMCH _dlamch; DNRM2 _dnrm2;
45:
46: #endregion
47:
48:
49: #region Fields
50:
51: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int ITEMP = 0; int J = 0; int K = 0; int LASTRK = 0;
52: int LSTICC = 0;int PVT = 0; int RK = 0; double AKK = 0; double TEMP = 0; double TEMP2 = 0; double TOL3Z = 0;
53:
54: #endregion
55:
56: public DLAQPS(DGEMM dgemm, DGEMV dgemv, DLARFG dlarfg, DSWAP dswap, IDAMAX idamax, DLAMCH dlamch, DNRM2 dnrm2)
57: {
58:
59:
60: #region Set Dependencies
61:
62: this._dgemm = dgemm; this._dgemv = dgemv; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax;
63: this._dlamch = dlamch;this._dnrm2 = dnrm2;
64:
65: #endregion
66:
67: }
68:
69: public DLAQPS()
70: {
71:
72:
73: #region Dependencies (Initialization)
74:
75: LSAME lsame = new LSAME();
76: XERBLA xerbla = new XERBLA();
77: DLAMC3 dlamc3 = new DLAMC3();
78: DLAPY2 dlapy2 = new DLAPY2();
79: DNRM2 dnrm2 = new DNRM2();
80: DSCAL dscal = new DSCAL();
81: DSWAP dswap = new DSWAP();
82: IDAMAX idamax = new IDAMAX();
83: DGEMM dgemm = new DGEMM(lsame, xerbla);
84: DGEMV dgemv = new DGEMV(lsame, xerbla);
85: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
86: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
87: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
88: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
89: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
90: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
91:
92: #endregion
93:
94:
95: #region Set Dependencies
96:
97: this._dgemm = dgemm; this._dgemv = dgemv; this._dlarfg = dlarfg; this._dswap = dswap; this._idamax = idamax;
98: this._dlamch = dlamch;this._dnrm2 = dnrm2;
99:
100: #endregion
101:
102: }
103: /// <summary>
104: /// Purpose
105: /// =======
106: ///
107: /// DLAQPS computes a step of QR factorization with column pivoting
108: /// of a real M-by-N matrix A by using Blas-3. It tries to factorize
109: /// NB columns from A starting from the row OFFSET+1, and updates all
110: /// of the matrix with Blas-3 xGEMM.
111: ///
112: /// In some cases, due to catastrophic cancellations, it cannot
113: /// factorize NB columns. Hence, the actual number of factorized
114: /// columns is returned in KB.
115: ///
116: /// Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
117: ///
118: ///</summary>
119: /// <param name="M">
120: /// (input) INTEGER
121: /// The number of rows of the matrix A. M .GE. 0.
122: ///</param>
123: /// <param name="N">
124: /// (input) INTEGER
125: /// The number of columns of the matrix A. N .GE. 0
126: ///</param>
127: /// <param name="OFFSET">
128: /// (input) INTEGER
129: /// The number of rows of A that have been factorized in
130: /// previous steps.
131: ///</param>
132: /// <param name="NB">
133: /// (input) INTEGER
134: /// The number of columns to factorize.
135: ///</param>
136: /// <param name="KB">
137: /// (output) INTEGER
138: /// The number of columns actually factorized.
139: ///</param>
140: /// <param name="A">
141: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
142: /// On entry, the M-by-N matrix A.
143: /// On exit, block A(OFFSET+1:M,1:KB) is the triangular
144: /// factor obtained and block A(1:OFFSET,1:N) has been
145: /// accordingly pivoted, but no factorized.
146: /// The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
147: /// been updated.
148: ///</param>
149: /// <param name="LDA">
150: /// (input) INTEGER
151: /// The leading dimension of the array A. LDA .GE. max(1,M).
152: ///</param>
153: /// <param name="JPVT">
154: /// (input/output) INTEGER array, dimension (N)
155: /// JPVT(I) = K .LE.=.GT. Column K of the full matrix A has been
156: /// permuted into position I in AP.
157: ///</param>
158: /// <param name="TAU">
159: /// (output) DOUBLE PRECISION array, dimension (KB)
160: /// The scalar factors of the elementary reflectors.
161: ///</param>
162: /// <param name="VN1">
163: /// (input/output) DOUBLE PRECISION array, dimension (N)
164: /// The vector with the partial column norms.
165: ///</param>
166: /// <param name="VN2">
167: /// (input/output) DOUBLE PRECISION array, dimension (N)
168: /// The vector with the exact column norms.
169: ///</param>
170: /// <param name="AUXV">
171: /// (input/output) DOUBLE PRECISION array, dimension (NB)
172: /// Auxiliar vector.
173: ///</param>
174: /// <param name="F">
175: /// (input/output) DOUBLE PRECISION array, dimension (LDF,NB)
176: /// Matrix F' = L*Y'*A.
177: ///</param>
178: /// <param name="LDF">
179: /// (input) INTEGER
180: /// The leading dimension of the array F. LDF .GE. max(1,N).
181: ///</param>
182: public void Run(int M, int N, int OFFSET, int NB, ref int KB, ref double[] A, int offset_a
183: , int LDA, ref int[] JPVT, int offset_jpvt, ref double[] TAU, int offset_tau, ref double[] VN1, int offset_vn1, ref double[] VN2, int offset_vn2, ref double[] AUXV, int offset_auxv
184: , ref double[] F, int offset_f, int LDF)
185: {
186:
187: #region Array Index Correction
188:
189: int o_a = -1 - LDA + offset_a; int o_jpvt = -1 + offset_jpvt; int o_tau = -1 + offset_tau;
190: int o_vn1 = -1 + offset_vn1; int o_vn2 = -1 + offset_vn2; int o_auxv = -1 + offset_auxv;
191: int o_f = -1 - LDF + offset_f;
192:
193: #endregion
194:
195:
196: #region Prolog
197:
198: // *
199: // * -- LAPACK auxiliary routine (version 3.1) --
200: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
201: // * November 2006
202: // *
203: // * .. Scalar Arguments ..
204: // * ..
205: // * .. Array Arguments ..
206: // * ..
207: // *
208: // * Purpose
209: // * =======
210: // *
211: // * DLAQPS computes a step of QR factorization with column pivoting
212: // * of a real M-by-N matrix A by using Blas-3. It tries to factorize
213: // * NB columns from A starting from the row OFFSET+1, and updates all
214: // * of the matrix with Blas-3 xGEMM.
215: // *
216: // * In some cases, due to catastrophic cancellations, it cannot
217: // * factorize NB columns. Hence, the actual number of factorized
218: // * columns is returned in KB.
219: // *
220: // * Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
221: // *
222: // * Arguments
223: // * =========
224: // *
225: // * M (input) INTEGER
226: // * The number of rows of the matrix A. M >= 0.
227: // *
228: // * N (input) INTEGER
229: // * The number of columns of the matrix A. N >= 0
230: // *
231: // * OFFSET (input) INTEGER
232: // * The number of rows of A that have been factorized in
233: // * previous steps.
234: // *
235: // * NB (input) INTEGER
236: // * The number of columns to factorize.
237: // *
238: // * KB (output) INTEGER
239: // * The number of columns actually factorized.
240: // *
241: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
242: // * On entry, the M-by-N matrix A.
243: // * On exit, block A(OFFSET+1:M,1:KB) is the triangular
244: // * factor obtained and block A(1:OFFSET,1:N) has been
245: // * accordingly pivoted, but no factorized.
246: // * The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has
247: // * been updated.
248: // *
249: // * LDA (input) INTEGER
250: // * The leading dimension of the array A. LDA >= max(1,M).
251: // *
252: // * JPVT (input/output) INTEGER array, dimension (N)
253: // * JPVT(I) = K <==> Column K of the full matrix A has been
254: // * permuted into position I in AP.
255: // *
256: // * TAU (output) DOUBLE PRECISION array, dimension (KB)
257: // * The scalar factors of the elementary reflectors.
258: // *
259: // * VN1 (input/output) DOUBLE PRECISION array, dimension (N)
260: // * The vector with the partial column norms.
261: // *
262: // * VN2 (input/output) DOUBLE PRECISION array, dimension (N)
263: // * The vector with the exact column norms.
264: // *
265: // * AUXV (input/output) DOUBLE PRECISION array, dimension (NB)
266: // * Auxiliar vector.
267: // *
268: // * F (input/output) DOUBLE PRECISION array, dimension (LDF,NB)
269: // * Matrix F' = L*Y'*A.
270: // *
271: // * LDF (input) INTEGER
272: // * The leading dimension of the array F. LDF >= max(1,N).
273: // *
274: // * Further Details
275: // * ===============
276: // *
277: // * Based on contributions by
278: // * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
279: // * X. Sun, Computer Science Dept., Duke University, USA
280: // *
281: // * Partial column norm updating strategy modified by
282: // * Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
283: // * University of Zagreb, Croatia.
284: // * June 2006.
285: // * For more details see LAPACK Working Note 176.
286: // * =====================================================================
287: // *
288: // * .. Parameters ..
289: // * ..
290: // * .. Local Scalars ..
291: // * ..
292: // * .. External Subroutines ..
293: // * ..
294: // * .. Intrinsic Functions ..
295: // INTRINSIC ABS, DBLE, MAX, MIN, NINT, SQRT;
296: // * ..
297: // * .. External Functions ..
298: // * ..
299: // * .. Executable Statements ..
300: // *
301:
302: #endregion
303:
304:
305: #region Body
306:
307: LASTRK = Math.Min(M, N + OFFSET);
308: LSTICC = 0;
309: K = 0;
310: TOL3Z = Math.Sqrt(this._dlamch.Run("Epsilon"));
311: // *
312: // * Beginning of while loop.
313: // *
314: LABEL10:;
315: if ((K < NB) && (LSTICC == 0))
316: {
317: K = K + 1;
318: RK = OFFSET + K;
319: // *
320: // * Determine ith pivot column and swap if necessary
321: // *
322: PVT = (K - 1) + this._idamax.Run(N - K + 1, VN1, K + o_vn1, 1);
323: if (PVT != K)
324: {
325: this._dswap.Run(M, ref A, 1+PVT * LDA + o_a, 1, ref A, 1+K * LDA + o_a, 1);
326: this._dswap.Run(K - 1, ref F, PVT+1 * LDF + o_f, LDF, ref F, K+1 * LDF + o_f, LDF);
327: ITEMP = JPVT[PVT + o_jpvt];
328: JPVT[PVT + o_jpvt] = JPVT[K + o_jpvt];
329: JPVT[K + o_jpvt] = ITEMP;
330: VN1[PVT + o_vn1] = VN1[K + o_vn1];
331: VN2[PVT + o_vn2] = VN2[K + o_vn2];
332: }
333: // *
334: // * Apply previous Householder reflectors to column K:
335: // * A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'.
336: // *
337: if (K > 1)
338: {
339: this._dgemv.Run("No transpose", M - RK + 1, K - 1, - ONE, A, RK+1 * LDA + o_a, LDA
340: , F, K+1 * LDF + o_f, LDF, ONE, ref A, RK+K * LDA + o_a, 1);
341: }
342: // *
343: // * Generate elementary reflector H(k).
344: // *
345: if (RK < M)
346: {
347: this._dlarfg.Run(M - RK + 1, ref A[RK+K * LDA + o_a], ref A, RK + 1+K * LDA + o_a, 1, ref TAU[K + o_tau]);
348: }
349: else
350: {
351: this._dlarfg.Run(1, ref A[RK+K * LDA + o_a], ref A, RK+K * LDA + o_a, 1, ref TAU[K + o_tau]);
352: }
353: // *
354: AKK = A[RK+K * LDA + o_a];
355: A[RK+K * LDA + o_a] = ONE;
356: // *
357: // * Compute Kth column of F:
358: // *
359: // * Compute F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K).
360: // *
361: if (K < N)
362: {
363: this._dgemv.Run("Transpose", M - RK + 1, N - K, TAU[K + o_tau], A, RK+(K + 1) * LDA + o_a, LDA
364: , A, RK+K * LDA + o_a, 1, ZERO, ref F, K + 1+K * LDF + o_f, 1);
365: }
366: // *
367: // * Padding F(1:K,K) with zeros.
368: // *
369: for (J = 1; J <= K; J++)
370: {
371: F[J+K * LDF + o_f] = ZERO;
372: }
373: // *
374: // * Incremental updating of F:
375: // * F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)'
376: // * *A(RK:M,K).
377: // *
378: if (K > 1)
379: {
380: this._dgemv.Run("Transpose", M - RK + 1, K - 1, - TAU[K + o_tau], A, RK+1 * LDA + o_a, LDA
381: , A, RK+K * LDA + o_a, 1, ZERO, ref AUXV, 1 + o_auxv, 1);
382: // *
383: this._dgemv.Run("No transpose", N, K - 1, ONE, F, 1+1 * LDF + o_f, LDF
384: , AUXV, 1 + o_auxv, 1, ONE, ref F, 1+K * LDF + o_f, 1);
385: }
386: // *
387: // * Update the current row of A:
388: // * A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'.
389: // *
390: if (K < N)
391: {
392: this._dgemv.Run("No transpose", N - K, K, - ONE, F, K + 1+1 * LDF + o_f, LDF
393: , A, RK+1 * LDA + o_a, LDA, ONE, ref A, RK+(K + 1) * LDA + o_a, LDA);
394: }
395: // *
396: // * Update partial column norms.
397: // *
398: if (RK < LASTRK)
399: {
400: for (J = K + 1; J <= N; J++)
401: {
402: if (VN1[J + o_vn1] != ZERO)
403: {
404: // *
405: // * NOTE: The following 4 lines follow from the analysis in
406: // * Lapack Working Note 176.
407: // *
408: TEMP = Math.Abs(A[RK+J * LDA + o_a]) / VN1[J + o_vn1];
409: TEMP = Math.Max(ZERO, (ONE + TEMP) * (ONE - TEMP));
410: TEMP2 = TEMP * Math.Pow(VN1[J + o_vn1] / VN2[J + o_vn2],2);
411: if (TEMP2 <= TOL3Z)
412: {
413: VN2[J + o_vn2] = Convert.ToDouble(LSTICC);
414: LSTICC = J;
415: }
416: else
417: {
418: VN1[J + o_vn1] = VN1[J + o_vn1] * Math.Sqrt(TEMP);
419: }
420: }
421: }
422: }
423: // *
424: A[RK+K * LDA + o_a] = AKK;
425: // *
426: // * End of while loop.
427: // *
428: goto LABEL10;
429: }
430: KB = K;
431: RK = OFFSET + KB;
432: // *
433: // * Apply the block reflector to the rest of the matrix:
434: // * A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
435: // * A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'.
436: // *
437: if (KB < Math.Min(N, M - OFFSET))
438: {
439: this._dgemm.Run("No transpose", "Transpose", M - RK, N - KB, KB, - ONE
440: , A, RK + 1+1 * LDA + o_a, LDA, F, KB + 1+1 * LDF + o_f, LDF, ONE, ref A, RK + 1+(KB + 1) * LDA + o_a
441: , LDA);
442: }
443: // *
444: // * Recomputation of difficult columns.
445: // *
446: LABEL40:;
447: if (LSTICC > 0)
448: {
449: ITEMP = (int)Math.Round(VN2[LSTICC + o_vn2]);
450: VN1[LSTICC + o_vn1] = this._dnrm2.Run(M - RK, A, RK + 1+LSTICC * LDA + o_a, 1);
451: // *
452: // * NOTE: The computation of VN1( LSTICC ) relies on the fact that
453: // * SNRM2 does not fail on vectors with norm below the value of
454: // * SQRT(DLAMCH('S'))
455: // *
456: VN2[LSTICC + o_vn2] = VN1[LSTICC + o_vn1];
457: LSTICC = ITEMP;
458: goto LABEL40;
459: }
460: // *
461: return;
462: // *
463: // * End of DLAQPS
464: // *
465:
466: #endregion
467:
468: }
469: }
470: }