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: /// DTZRZF reduces the M-by-N ( M.LE.N ) real upper trapezoidal matrix A
27: /// to upper triangular form by means of orthogonal transformations.
28: ///
29: /// The upper trapezoidal matrix A is factored as
30: ///
31: /// A = ( R 0 ) * Z,
32: ///
33: /// where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
34: /// triangular matrix.
35: ///
36: ///</summary>
37: public class DTZRZF
38: {
39:
40:
41: #region Dependencies
42:
43: DLARZB _dlarzb; DLARZT _dlarzt; DLATRZ _dlatrz; XERBLA _xerbla; ILAENV _ilaenv;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const double ZERO = 0.0E+0; bool LQUERY = false; int I = 0; int IB = 0; int IWS = 0; int KI = 0; int KK = 0;
51: int LDWORK = 0;int LWKOPT = 0; int M1 = 0; int MU = 0; int NB = 0; int NBMIN = 0; int NX = 0;
52:
53: #endregion
54:
55: public DTZRZF(DLARZB dlarzb, DLARZT dlarzt, DLATRZ dlatrz, XERBLA xerbla, ILAENV ilaenv)
56: {
57:
58:
59: #region Set Dependencies
60:
61: this._dlarzb = dlarzb; this._dlarzt = dlarzt; this._dlatrz = dlatrz; this._xerbla = xerbla; this._ilaenv = ilaenv;
62:
63: #endregion
64:
65: }
66:
67: public DTZRZF()
68: {
69:
70:
71: #region Dependencies (Initialization)
72:
73: LSAME lsame = new LSAME();
74: DCOPY dcopy = new DCOPY();
75: XERBLA xerbla = new XERBLA();
76: DLAMC3 dlamc3 = new DLAMC3();
77: DLAPY2 dlapy2 = new DLAPY2();
78: DNRM2 dnrm2 = new DNRM2();
79: DSCAL dscal = new DSCAL();
80: DAXPY daxpy = new DAXPY();
81: IEEECK ieeeck = new IEEECK();
82: IPARMQ iparmq = new IPARMQ();
83: DGEMM dgemm = new DGEMM(lsame, xerbla);
84: DTRMM dtrmm = new DTRMM(lsame, xerbla);
85: DLARZB dlarzb = new DLARZB(lsame, dcopy, dgemm, dtrmm, xerbla);
86: DGEMV dgemv = new DGEMV(lsame, xerbla);
87: DTRMV dtrmv = new DTRMV(lsame, xerbla);
88: DLARZT dlarzt = new DLARZT(dgemv, dtrmv, xerbla, lsame);
89: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
90: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
91: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
92: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
93: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
94: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
95: DGER dger = new DGER(xerbla);
96: DLARZ dlarz = new DLARZ(daxpy, dcopy, dgemv, dger, lsame);
97: DLATRZ dlatrz = new DLATRZ(dlarfg, dlarz);
98: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
99:
100: #endregion
101:
102:
103: #region Set Dependencies
104:
105: this._dlarzb = dlarzb; this._dlarzt = dlarzt; this._dlatrz = dlatrz; this._xerbla = xerbla; this._ilaenv = ilaenv;
106:
107: #endregion
108:
109: }
110: /// <summary>
111: /// Purpose
112: /// =======
113: ///
114: /// DTZRZF reduces the M-by-N ( M.LE.N ) real upper trapezoidal matrix A
115: /// to upper triangular form by means of orthogonal transformations.
116: ///
117: /// The upper trapezoidal matrix A is factored as
118: ///
119: /// A = ( R 0 ) * Z,
120: ///
121: /// where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
122: /// triangular matrix.
123: ///
124: ///</summary>
125: /// <param name="M">
126: /// (input) INTEGER
127: /// The number of rows of the matrix A. M .GE. 0.
128: ///</param>
129: /// <param name="N">
130: /// (input) INTEGER
131: /// The number of columns of the matrix A. N .GE. M.
132: ///</param>
133: /// <param name="A">
134: /// = ( R 0 ) * Z,
135: ///</param>
136: /// <param name="LDA">
137: /// (input) INTEGER
138: /// The leading dimension of the array A. LDA .GE. max(1,M).
139: ///</param>
140: /// <param name="TAU">
141: /// (output) DOUBLE PRECISION array, dimension (M)
142: /// The scalar factors of the elementary reflectors.
143: ///</param>
144: /// <param name="WORK">
145: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
146: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
147: ///</param>
148: /// <param name="LWORK">
149: /// (input) INTEGER
150: /// The dimension of the array WORK. LWORK .GE. max(1,M).
151: /// For optimum performance LWORK .GE. M*NB, where NB is
152: /// the optimal blocksize.
153: ///
154: /// If LWORK = -1, then a workspace query is assumed; the routine
155: /// only calculates the optimal size of the WORK array, returns
156: /// this value as the first entry of the WORK array, and no error
157: /// message related to LWORK is issued by XERBLA.
158: ///</param>
159: /// <param name="INFO">
160: /// (output) INTEGER
161: /// = 0: successful exit
162: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
163: ///</param>
164: 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
165: , int LWORK, ref int INFO)
166: {
167:
168: #region Array Index Correction
169:
170: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
171:
172: #endregion
173:
174:
175: #region Prolog
176:
177: // *
178: // * -- LAPACK routine (version 3.1) --
179: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
180: // * November 2006
181: // *
182: // * .. Scalar Arguments ..
183: // * ..
184: // * .. Array Arguments ..
185: // * ..
186: // *
187: // * Purpose
188: // * =======
189: // *
190: // * DTZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A
191: // * to upper triangular form by means of orthogonal transformations.
192: // *
193: // * The upper trapezoidal matrix A is factored as
194: // *
195: // * A = ( R 0 ) * Z,
196: // *
197: // * where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
198: // * triangular matrix.
199: // *
200: // * Arguments
201: // * =========
202: // *
203: // * M (input) INTEGER
204: // * The number of rows of the matrix A. M >= 0.
205: // *
206: // * N (input) INTEGER
207: // * The number of columns of the matrix A. N >= M.
208: // *
209: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
210: // * On entry, the leading M-by-N upper trapezoidal part of the
211: // * array A must contain the matrix to be factorized.
212: // * On exit, the leading M-by-M upper triangular part of A
213: // * contains the upper triangular matrix R, and elements M+1 to
214: // * N of the first M rows of A, with the array TAU, represent the
215: // * orthogonal matrix Z as a product of M elementary reflectors.
216: // *
217: // * LDA (input) INTEGER
218: // * The leading dimension of the array A. LDA >= max(1,M).
219: // *
220: // * TAU (output) DOUBLE PRECISION array, dimension (M)
221: // * The scalar factors of the elementary reflectors.
222: // *
223: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
224: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
225: // *
226: // * LWORK (input) INTEGER
227: // * The dimension of the array WORK. LWORK >= max(1,M).
228: // * For optimum performance LWORK >= M*NB, where NB is
229: // * the optimal blocksize.
230: // *
231: // * If LWORK = -1, then a workspace query is assumed; the routine
232: // * only calculates the optimal size of the WORK array, returns
233: // * this value as the first entry of the WORK array, and no error
234: // * message related to LWORK is issued by XERBLA.
235: // *
236: // * INFO (output) INTEGER
237: // * = 0: successful exit
238: // * < 0: if INFO = -i, the i-th argument had an illegal value
239: // *
240: // * Further Details
241: // * ===============
242: // *
243: // * Based on contributions by
244: // * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
245: // *
246: // * The factorization is obtained by Householder's method. The kth
247: // * transformation matrix, Z( k ), which is used to introduce zeros into
248: // * the ( m - k + 1 )th row of A, is given in the form
249: // *
250: // * Z( k ) = ( I 0 ),
251: // * ( 0 T( k ) )
252: // *
253: // * where
254: // *
255: // * T( k ) = I - tau*u( k )*u( k )', u( k ) = ( 1 ),
256: // * ( 0 )
257: // * ( z( k ) )
258: // *
259: // * tau is a scalar and z( k ) is an ( n - m ) element vector.
260: // * tau and z( k ) are chosen to annihilate the elements of the kth row
261: // * of X.
262: // *
263: // * The scalar tau is returned in the kth element of TAU and the vector
264: // * u( k ) in the kth row of A, such that the elements of z( k ) are
265: // * in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
266: // * the upper triangular part of A.
267: // *
268: // * Z is given by
269: // *
270: // * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ).
271: // *
272: // * =====================================================================
273: // *
274: // * .. Parameters ..
275: // * ..
276: // * .. Local Scalars ..
277: // * ..
278: // * .. External Subroutines ..
279: // * ..
280: // * .. Intrinsic Functions ..
281: // INTRINSIC MAX, MIN;
282: // * ..
283: // * .. External Functions ..
284: // * ..
285: // * .. Executable Statements ..
286: // *
287: // * Test the input arguments
288: // *
289:
290: #endregion
291:
292:
293: #region Body
294:
295: INFO = 0;
296: LQUERY = (LWORK == - 1);
297: if (M < 0)
298: {
299: INFO = - 1;
300: }
301: else
302: {
303: if (N < M)
304: {
305: INFO = - 2;
306: }
307: else
308: {
309: if (LDA < Math.Max(1, M))
310: {
311: INFO = - 4;
312: }
313: }
314: }
315: // *
316: if (INFO == 0)
317: {
318: if (M == 0 || M == N)
319: {
320: LWKOPT = 1;
321: }
322: else
323: {
324: // *
325: // * Determine the block size.
326: // *
327: NB = this._ilaenv.Run(1, "DGERQF", " ", M, N, - 1, - 1);
328: LWKOPT = M * NB;
329: }
330: WORK[1 + o_work] = LWKOPT;
331: // *
332: if (LWORK < Math.Max(1, M) && !LQUERY)
333: {
334: INFO = - 7;
335: }
336: }
337: // *
338: if (INFO != 0)
339: {
340: this._xerbla.Run("DTZRZF", - INFO);
341: return;
342: }
343: else
344: {
345: if (LQUERY)
346: {
347: return;
348: }
349: }
350: // *
351: // * Quick return if possible
352: // *
353: if (M == 0)
354: {
355: return;
356: }
357: else
358: {
359: if (M == N)
360: {
361: for (I = 1; I <= N; I++)
362: {
363: TAU[I + o_tau] = ZERO;
364: }
365: return;
366: }
367: }
368: // *
369: NBMIN = 2;
370: NX = 1;
371: IWS = M;
372: if (NB > 1 && NB < M)
373: {
374: // *
375: // * Determine when to cross over from blocked to unblocked code.
376: // *
377: NX = Math.Max(0, this._ilaenv.Run(3, "DGERQF", " ", M, N, - 1, - 1));
378: if (NX < M)
379: {
380: // *
381: // * Determine if workspace is large enough for blocked code.
382: // *
383: LDWORK = M;
384: IWS = LDWORK * NB;
385: if (LWORK < IWS)
386: {
387: // *
388: // * Not enough workspace to use optimal NB: reduce NB and
389: // * determine the minimum value of NB.
390: // *
391: NB = LWORK / LDWORK;
392: NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGERQF", " ", M, N, - 1, - 1));
393: }
394: }
395: }
396: // *
397: if (NB >= NBMIN && NB < M && NX < M)
398: {
399: // *
400: // * Use blocked code initially.
401: // * The last kk rows are handled by the block method.
402: // *
403: M1 = Math.Min(M + 1, N);
404: KI = ((M - NX - 1) / NB) * NB;
405: KK = Math.Min(M, KI + NB);
406: // *
407: for (I = M - KK + KI + 1; ( - NB >= 0) ? (I <= M - KK + 1) : (I >= M - KK + 1); I += - NB)
408: {
409: IB = Math.Min(M - I + 1, NB);
410: // *
411: // * Compute the TZ factorization of the current block
412: // * A(i:i+ib-1,i:n)
413: // *
414: this._dlatrz.Run(IB, N - I + 1, N - M, ref A, I+I * LDA + o_a, LDA, ref TAU, I + o_tau
415: , ref WORK, offset_work);
416: if (I > 1)
417: {
418: // *
419: // * Form the triangular factor of the block reflector
420: // * H = H(i+ib-1) . . . H(i+1) H(i)
421: // *
422: this._dlarzt.Run("Backward", "Rowwise", N - M, IB, A, I+M1 * LDA + o_a, LDA
423: , TAU, I + o_tau, ref WORK, offset_work, LDWORK);
424: // *
425: // * Apply H to A(1:i-1,i:n) from the right
426: // *
427: this._dlarzb.Run("Right", "No transpose", "Backward", "Rowwise", I - 1, N - I + 1
428: , IB, N - M, A, I+M1 * LDA + o_a, LDA, WORK, offset_work, LDWORK
429: , ref A, 1+I * LDA + o_a, LDA, ref WORK, IB + 1 + o_work, LDWORK);
430: }
431: }
432: MU = I + NB - 1;
433: }
434: else
435: {
436: MU = M;
437: }
438: // *
439: // * Use unblocked code to factor the last or only block
440: // *
441: if (MU > 0)
442: {
443: this._dlatrz.Run(MU, N, N - M, ref A, offset_a, LDA, ref TAU, offset_tau
444: , ref WORK, offset_work);
445: }
446: // *
447: WORK[1 + o_work] = LWKOPT;
448: // *
449: return;
450: // *
451: // * End of DTZRZF
452: // *
453:
454: #endregion
455:
456: }
457: }
458: }