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