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