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.1) --
21: /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
22: /// January 2007
23: /// Purpose
24: /// =======
25: ///
26: /// DORMRZ overwrites the general real M-by-N matrix C with
27: ///
28: /// SIDE = 'L' SIDE = 'R'
29: /// TRANS = 'N': Q * C C * Q
30: /// TRANS = 'T': Q**T * C C * Q**T
31: ///
32: /// where Q is a real orthogonal matrix defined as the product of k
33: /// elementary reflectors
34: ///
35: /// Q = H(1) H(2) . . . H(k)
36: ///
37: /// as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N
38: /// if SIDE = 'R'.
39: ///
40: ///</summary>
41: public class DORMRZ
42: {
43:
44:
45: #region Dependencies
46:
47: LSAME _lsame; ILAENV _ilaenv; DLARZB _dlarzb; DLARZT _dlarzt; DORMR3 _dormr3; XERBLA _xerbla;
48:
49: #endregion
50:
51:
52: #region Fields
53:
54: const int NBMAX = 64; const int LDT = NBMAX + 1; bool LEFT = false; bool LQUERY = false; bool NOTRAN = false;
55: string TRANST = new string(' ', 1);int I = 0; int I1 = 0; int I2 = 0; int I3 = 0; int IB = 0; int IC = 0; int IINFO = 0;
56: int IWS = 0;int JA = 0; int JC = 0; int LDWORK = 0; int LWKOPT = 0; int MI = 0; int NB = 0; int NBMIN = 0; int NI = 0;
57: int NQ = 0;int NW = 0; double[] T = new double[LDT * NBMAX]; int offset_t = 0;
58:
59: #endregion
60:
61: public DORMRZ(LSAME lsame, ILAENV ilaenv, DLARZB dlarzb, DLARZT dlarzt, DORMR3 dormr3, XERBLA xerbla)
62: {
63:
64:
65: #region Set Dependencies
66:
67: this._lsame = lsame; this._ilaenv = ilaenv; this._dlarzb = dlarzb; this._dlarzt = dlarzt; this._dormr3 = dormr3;
68: this._xerbla = xerbla;
69:
70: #endregion
71:
72: }
73:
74: public DORMRZ()
75: {
76:
77:
78: #region Dependencies (Initialization)
79:
80: LSAME lsame = new LSAME();
81: IEEECK ieeeck = new IEEECK();
82: IPARMQ iparmq = new IPARMQ();
83: DCOPY dcopy = new DCOPY();
84: XERBLA xerbla = new XERBLA();
85: DAXPY daxpy = new DAXPY();
86: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
87: DGEMM dgemm = new DGEMM(lsame, xerbla);
88: DTRMM dtrmm = new DTRMM(lsame, xerbla);
89: DLARZB dlarzb = new DLARZB(lsame, dcopy, dgemm, dtrmm, xerbla);
90: DGEMV dgemv = new DGEMV(lsame, xerbla);
91: DTRMV dtrmv = new DTRMV(lsame, xerbla);
92: DLARZT dlarzt = new DLARZT(dgemv, dtrmv, xerbla, lsame);
93: DGER dger = new DGER(xerbla);
94: DLARZ dlarz = new DLARZ(daxpy, dcopy, dgemv, dger, lsame);
95: DORMR3 dormr3 = new DORMR3(lsame, dlarz, xerbla);
96:
97: #endregion
98:
99:
100: #region Set Dependencies
101:
102: this._lsame = lsame; this._ilaenv = ilaenv; this._dlarzb = dlarzb; this._dlarzt = dlarzt; this._dormr3 = dormr3;
103: this._xerbla = xerbla;
104:
105: #endregion
106:
107: }
108: /// <summary>
109: /// Purpose
110: /// =======
111: ///
112: /// DORMRZ overwrites the general real M-by-N matrix C with
113: ///
114: /// SIDE = 'L' SIDE = 'R'
115: /// TRANS = 'N': Q * C C * Q
116: /// TRANS = 'T': Q**T * C C * Q**T
117: ///
118: /// where Q is a real orthogonal matrix defined as the product of k
119: /// elementary reflectors
120: ///
121: /// Q = H(1) H(2) . . . H(k)
122: ///
123: /// as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N
124: /// if SIDE = 'R'.
125: ///
126: ///</summary>
127: /// <param name="SIDE">
128: /// = 'L' SIDE = 'R'
129: ///</param>
130: /// <param name="TRANS">
131: /// (input) CHARACTER*1
132: /// = 'N': No transpose, apply Q;
133: /// = 'T': Transpose, apply Q**T.
134: ///</param>
135: /// <param name="M">
136: /// (input) INTEGER
137: /// The number of rows of the matrix C. M .GE. 0.
138: ///</param>
139: /// <param name="N">
140: /// (input) INTEGER
141: /// The number of columns of the matrix C. N .GE. 0.
142: ///</param>
143: /// <param name="K">
144: /// (input) INTEGER
145: /// The number of elementary reflectors whose product defines
146: /// the matrix Q.
147: /// If SIDE = 'L', M .GE. K .GE. 0;
148: /// if SIDE = 'R', N .GE. K .GE. 0.
149: ///</param>
150: /// <param name="L">
151: /// (input) INTEGER
152: /// The number of columns of the matrix A containing
153: /// the meaningful part of the Householder reflectors.
154: /// If SIDE = 'L', M .GE. L .GE. 0, if SIDE = 'R', N .GE. L .GE. 0.
155: ///</param>
156: /// <param name="A">
157: /// (input) DOUBLE PRECISION array, dimension
158: /// (LDA,M) if SIDE = 'L',
159: /// (LDA,N) if SIDE = 'R'
160: /// The i-th row must contain the vector which defines the
161: /// elementary reflector H(i), for i = 1,2,...,k, as returned by
162: /// DTZRZF in the last k rows of its array argument A.
163: /// A is modified by the routine but restored on exit.
164: ///</param>
165: /// <param name="LDA">
166: /// (input) INTEGER
167: /// The leading dimension of the array A. LDA .GE. max(1,K).
168: ///</param>
169: /// <param name="TAU">
170: /// (input) DOUBLE PRECISION array, dimension (K)
171: /// TAU(i) must contain the scalar factor of the elementary
172: /// reflector H(i), as returned by DTZRZF.
173: ///</param>
174: /// <param name="C">
175: /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
176: /// On entry, the M-by-N matrix C.
177: /// On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
178: ///</param>
179: /// <param name="LDC">
180: /// (input) INTEGER
181: /// The leading dimension of the array C. LDC .GE. max(1,M).
182: ///</param>
183: /// <param name="WORK">
184: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
185: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
186: ///</param>
187: /// <param name="LWORK">
188: /// (input) INTEGER
189: /// The dimension of the array WORK.
190: /// If SIDE = 'L', LWORK .GE. max(1,N);
191: /// if SIDE = 'R', LWORK .GE. max(1,M).
192: /// For optimum performance LWORK .GE. N*NB if SIDE = 'L', and
193: /// LWORK .GE. M*NB if SIDE = 'R', where NB is the optimal
194: /// blocksize.
195: ///
196: /// If LWORK = -1, then a workspace query is assumed; the routine
197: /// only calculates the optimal size of the WORK array, returns
198: /// this value as the first entry of the WORK array, and no error
199: /// message related to LWORK is issued by XERBLA.
200: ///</param>
201: /// <param name="INFO">
202: /// (output) INTEGER
203: /// = 0: successful exit
204: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
205: ///</param>
206: public void Run(string SIDE, string TRANS, int M, int N, int K, int L
207: , double[] A, int offset_a, int LDA, double[] TAU, int offset_tau, ref double[] C, int offset_c, int LDC, ref double[] WORK, int offset_work
208: , int LWORK, ref int INFO)
209: {
210:
211: #region Array Index Correction
212:
213: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_c = -1 - LDC + offset_c;
214: int o_work = -1 + offset_work;
215:
216: #endregion
217:
218:
219: #region Strings
220:
221: SIDE = SIDE.Substring(0, 1); TRANS = TRANS.Substring(0, 1);
222:
223: #endregion
224:
225:
226: #region Prolog
227:
228: // *
229: // * -- LAPACK routine (version 3.1.1) --
230: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
231: // * January 2007
232: // *
233: // * .. Scalar Arguments ..
234: // * ..
235: // * .. Array Arguments ..
236: // * ..
237: // *
238: // * Purpose
239: // * =======
240: // *
241: // * DORMRZ overwrites the general real M-by-N matrix C with
242: // *
243: // * SIDE = 'L' SIDE = 'R'
244: // * TRANS = 'N': Q * C C * Q
245: // * TRANS = 'T': Q**T * C C * Q**T
246: // *
247: // * where Q is a real orthogonal matrix defined as the product of k
248: // * elementary reflectors
249: // *
250: // * Q = H(1) H(2) . . . H(k)
251: // *
252: // * as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N
253: // * if SIDE = 'R'.
254: // *
255: // * Arguments
256: // * =========
257: // *
258: // * SIDE (input) CHARACTER*1
259: // * = 'L': apply Q or Q**T from the Left;
260: // * = 'R': apply Q or Q**T from the Right.
261: // *
262: // * TRANS (input) CHARACTER*1
263: // * = 'N': No transpose, apply Q;
264: // * = 'T': Transpose, apply Q**T.
265: // *
266: // * M (input) INTEGER
267: // * The number of rows of the matrix C. M >= 0.
268: // *
269: // * N (input) INTEGER
270: // * The number of columns of the matrix C. N >= 0.
271: // *
272: // * K (input) INTEGER
273: // * The number of elementary reflectors whose product defines
274: // * the matrix Q.
275: // * If SIDE = 'L', M >= K >= 0;
276: // * if SIDE = 'R', N >= K >= 0.
277: // *
278: // * L (input) INTEGER
279: // * The number of columns of the matrix A containing
280: // * the meaningful part of the Householder reflectors.
281: // * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
282: // *
283: // * A (input) DOUBLE PRECISION array, dimension
284: // * (LDA,M) if SIDE = 'L',
285: // * (LDA,N) if SIDE = 'R'
286: // * The i-th row must contain the vector which defines the
287: // * elementary reflector H(i), for i = 1,2,...,k, as returned by
288: // * DTZRZF in the last k rows of its array argument A.
289: // * A is modified by the routine but restored on exit.
290: // *
291: // * LDA (input) INTEGER
292: // * The leading dimension of the array A. LDA >= max(1,K).
293: // *
294: // * TAU (input) DOUBLE PRECISION array, dimension (K)
295: // * TAU(i) must contain the scalar factor of the elementary
296: // * reflector H(i), as returned by DTZRZF.
297: // *
298: // * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
299: // * On entry, the M-by-N matrix C.
300: // * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
301: // *
302: // * LDC (input) INTEGER
303: // * The leading dimension of the array C. LDC >= max(1,M).
304: // *
305: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
306: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
307: // *
308: // * LWORK (input) INTEGER
309: // * The dimension of the array WORK.
310: // * If SIDE = 'L', LWORK >= max(1,N);
311: // * if SIDE = 'R', LWORK >= max(1,M).
312: // * For optimum performance LWORK >= N*NB if SIDE = 'L', and
313: // * LWORK >= M*NB if SIDE = 'R', where NB is the optimal
314: // * blocksize.
315: // *
316: // * If LWORK = -1, then a workspace query is assumed; the routine
317: // * only calculates the optimal size of the WORK array, returns
318: // * this value as the first entry of the WORK array, and no error
319: // * message related to LWORK is issued by XERBLA.
320: // *
321: // * INFO (output) INTEGER
322: // * = 0: successful exit
323: // * < 0: if INFO = -i, the i-th argument had an illegal value
324: // *
325: // * Further Details
326: // * ===============
327: // *
328: // * Based on contributions by
329: // * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
330: // *
331: // * =====================================================================
332: // *
333: // * .. Parameters ..
334: // * ..
335: // * .. Local Scalars ..
336: // * ..
337: // * .. Local Arrays ..
338: // * ..
339: // * .. External Functions ..
340: // * ..
341: // * .. External Subroutines ..
342: // * ..
343: // * .. Intrinsic Functions ..
344: // INTRINSIC MAX, MIN;
345: // * ..
346: // * .. Executable Statements ..
347: // *
348: // * Test the input arguments
349: // *
350:
351: #endregion
352:
353:
354: #region Body
355:
356: INFO = 0;
357: LEFT = this._lsame.Run(SIDE, "L");
358: NOTRAN = this._lsame.Run(TRANS, "N");
359: LQUERY = (LWORK == - 1);
360: // *
361: // * NQ is the order of Q and NW is the minimum dimension of WORK
362: // *
363: if (LEFT)
364: {
365: NQ = M;
366: NW = Math.Max(1, N);
367: }
368: else
369: {
370: NQ = N;
371: NW = Math.Max(1, M);
372: }
373: if (!LEFT && !this._lsame.Run(SIDE, "R"))
374: {
375: INFO = - 1;
376: }
377: else
378: {
379: if (!NOTRAN && !this._lsame.Run(TRANS, "T"))
380: {
381: INFO = - 2;
382: }
383: else
384: {
385: if (M < 0)
386: {
387: INFO = - 3;
388: }
389: else
390: {
391: if (N < 0)
392: {
393: INFO = - 4;
394: }
395: else
396: {
397: if (K < 0 || K > NQ)
398: {
399: INFO = - 5;
400: }
401: else
402: {
403: if (L < 0 || (LEFT && (L > M)) || (!LEFT && (L > N)))
404: {
405: INFO = - 6;
406: }
407: else
408: {
409: if (LDA < Math.Max(1, K))
410: {
411: INFO = - 8;
412: }
413: else
414: {
415: if (LDC < Math.Max(1, M))
416: {
417: INFO = - 11;
418: }
419: }
420: }
421: }
422: }
423: }
424: }
425: }
426: // *
427: if (INFO == 0)
428: {
429: if (M == 0 || N == 0)
430: {
431: LWKOPT = 1;
432: }
433: else
434: {
435: // *
436: // * Determine the block size. NB may be at most NBMAX, where
437: // * NBMAX is used to define the local array T.
438: // *
439: NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DORMRQ", SIDE + TRANS, M, N, K, - 1));
440: LWKOPT = NW * NB;
441: }
442: WORK[1 + o_work] = LWKOPT;
443: // *
444: if (LWORK < Math.Max(1, NW) && !LQUERY)
445: {
446: INFO = - 13;
447: }
448: }
449: // *
450: if (INFO != 0)
451: {
452: this._xerbla.Run("DORMRZ", - INFO);
453: return;
454: }
455: else
456: {
457: if (LQUERY)
458: {
459: return;
460: }
461: }
462: // *
463: // * Quick return if possible
464: // *
465: if (M == 0 || N == 0)
466: {
467: WORK[1 + o_work] = 1;
468: return;
469: }
470: // *
471: NBMIN = 2;
472: LDWORK = NW;
473: if (NB > 1 && NB < K)
474: {
475: IWS = NW * NB;
476: if (LWORK < IWS)
477: {
478: NB = LWORK / LDWORK;
479: NBMIN = Math.Max(2, this._ilaenv.Run(2, "DORMRQ", SIDE + TRANS, M, N, K, - 1));
480: }
481: }
482: else
483: {
484: IWS = NW;
485: }
486: // *
487: if (NB < NBMIN || NB >= K)
488: {
489: // *
490: // * Use unblocked code
491: // *
492: this._dormr3.Run(SIDE, TRANS, M, N, K, L
493: , A, offset_a, LDA, TAU, offset_tau, ref C, offset_c, LDC, ref WORK, offset_work
494: , ref IINFO);
495: }
496: else
497: {
498: // *
499: // * Use blocked code
500: // *
501: if ((LEFT && !NOTRAN) || (!LEFT && NOTRAN))
502: {
503: I1 = 1;
504: I2 = K;
505: I3 = NB;
506: }
507: else
508: {
509: I1 = ((K - 1) / NB) * NB + 1;
510: I2 = 1;
511: I3 = - NB;
512: }
513: // *
514: if (LEFT)
515: {
516: NI = N;
517: JC = 1;
518: JA = M - L + 1;
519: }
520: else
521: {
522: MI = M;
523: IC = 1;
524: JA = N - L + 1;
525: }
526: // *
527: if (NOTRAN)
528: {
529: FortranLib.Copy(ref TRANST , "T");
530: }
531: else
532: {
533: FortranLib.Copy(ref TRANST , "N");
534: }
535: // *
536: for (I = I1; (I3 >= 0) ? (I <= I2) : (I >= I2); I += I3)
537: {
538: IB = Math.Min(NB, K - I + 1);
539: // *
540: // * Form the triangular factor of the block reflector
541: // * H = H(i+ib-1) . . . H(i+1) H(i)
542: // *
543: this._dlarzt.Run("Backward", "Rowwise", L, IB, A, I+JA * LDA + o_a, LDA
544: , TAU, I + o_tau, ref T, offset_t, LDT);
545: // *
546: if (LEFT)
547: {
548: // *
549: // * H or H' is applied to C(i:m,1:n)
550: // *
551: MI = M - I + 1;
552: IC = I;
553: }
554: else
555: {
556: // *
557: // * H or H' is applied to C(1:m,i:n)
558: // *
559: NI = N - I + 1;
560: JC = I;
561: }
562: // *
563: // * Apply H or H'
564: // *
565: this._dlarzb.Run(SIDE, TRANST, "Backward", "Rowwise", MI, NI
566: , IB, L, A, I+JA * LDA + o_a, LDA, T, offset_t, LDT
567: , ref C, IC+JC * LDC + o_c, LDC, ref WORK, offset_work, LDWORK);
568: }
569: // *
570: }
571: // *
572: WORK[1 + o_work] = LWKOPT;
573: // *
574: return;
575: // *
576: // * End of DORMRZ
577: // *
578:
579: #endregion
580:
581: }
582: }
583: }