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: /// DORMTR 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 of order nq, with nq = m if
33: /// SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
34: /// nq-1 elementary reflectors, as returned by DSYTRD:
35: ///
36: /// if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
37: ///
38: /// if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
39: ///
40: ///</summary>
41: public class DORMTR
42: {
43:
44:
45: #region Dependencies
46:
47: LSAME _lsame; ILAENV _ilaenv; DORMQL _dormql; DORMQR _dormqr; XERBLA _xerbla;
48:
49: #endregion
50:
51:
52: #region Fields
53:
54: bool LEFT = false; bool LQUERY = false; bool UPPER = false; int I1 = 0; int I2 = 0; int IINFO = 0; int LWKOPT = 0;
55: int MI = 0;int NB = 0; int NI = 0; int NQ = 0; int NW = 0;
56:
57: #endregion
58:
59: public DORMTR(LSAME lsame, ILAENV ilaenv, DORMQL dormql, DORMQR dormqr, XERBLA xerbla)
60: {
61:
62:
63: #region Set Dependencies
64:
65: this._lsame = lsame; this._ilaenv = ilaenv; this._dormql = dormql; this._dormqr = dormqr; this._xerbla = xerbla;
66:
67: #endregion
68:
69: }
70:
71: public DORMTR()
72: {
73:
74:
75: #region Dependencies (Initialization)
76:
77: LSAME lsame = new LSAME();
78: IEEECK ieeeck = new IEEECK();
79: IPARMQ iparmq = new IPARMQ();
80: DCOPY dcopy = new DCOPY();
81: XERBLA xerbla = new XERBLA();
82: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
83: DGEMM dgemm = new DGEMM(lsame, xerbla);
84: DTRMM dtrmm = new DTRMM(lsame, xerbla);
85: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
86: DGEMV dgemv = new DGEMV(lsame, xerbla);
87: DTRMV dtrmv = new DTRMV(lsame, xerbla);
88: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
89: DGER dger = new DGER(xerbla);
90: DLARF dlarf = new DLARF(dgemv, dger, lsame);
91: DORM2L dorm2l = new DORM2L(lsame, dlarf, xerbla);
92: DORMQL dormql = new DORMQL(lsame, ilaenv, dlarfb, dlarft, dorm2l, xerbla);
93: DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
94: DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
95:
96: #endregion
97:
98:
99: #region Set Dependencies
100:
101: this._lsame = lsame; this._ilaenv = ilaenv; this._dormql = dormql; this._dormqr = dormqr; this._xerbla = xerbla;
102:
103: #endregion
104:
105: }
106: /// <summary>
107: /// Purpose
108: /// =======
109: ///
110: /// DORMTR 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 of order nq, with nq = m if
117: /// SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
118: /// nq-1 elementary reflectors, as returned by DSYTRD:
119: ///
120: /// if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
121: ///
122: /// if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
123: ///
124: ///</summary>
125: /// <param name="SIDE">
126: /// = 'L' SIDE = 'R'
127: ///</param>
128: /// <param name="UPLO">
129: /// (input) CHARACTER*1
130: /// = 'U': Upper triangle of A contains elementary reflectors
131: /// from DSYTRD;
132: /// = 'L': Lower triangle of A contains elementary reflectors
133: /// from DSYTRD.
134: ///</param>
135: /// <param name="TRANS">
136: /// (input) CHARACTER*1
137: /// = 'N': No transpose, apply Q;
138: /// = 'T': Transpose, apply Q**T.
139: ///</param>
140: /// <param name="M">
141: /// (input) INTEGER
142: /// The number of rows of the matrix C. M .GE. 0.
143: ///</param>
144: /// <param name="N">
145: /// (input) INTEGER
146: /// The number of columns of the matrix C. N .GE. 0.
147: ///</param>
148: /// <param name="A">
149: /// (input) DOUBLE PRECISION array, dimension
150: /// (LDA,M) if SIDE = 'L'
151: /// (LDA,N) if SIDE = 'R'
152: /// The vectors which define the elementary reflectors, as
153: /// returned by DSYTRD.
154: ///</param>
155: /// <param name="LDA">
156: /// (input) INTEGER
157: /// The leading dimension of the array A.
158: /// LDA .GE. max(1,M) if SIDE = 'L'; LDA .GE. max(1,N) if SIDE = 'R'.
159: ///</param>
160: /// <param name="TAU">
161: /// (input) DOUBLE PRECISION array, dimension
162: /// (M-1) if SIDE = 'L'
163: /// (N-1) if SIDE = 'R'
164: /// TAU(i) must contain the scalar factor of the elementary
165: /// reflector H(i), as returned by DSYTRD.
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 UPLO, string TRANS, int M, int N, 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); UPLO = UPLO.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: // * DORMTR 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 of order nq, with nq = m if
241: // * SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
242: // * nq-1 elementary reflectors, as returned by DSYTRD:
243: // *
244: // * if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
245: // *
246: // * if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
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: // * UPLO (input) CHARACTER*1
256: // * = 'U': Upper triangle of A contains elementary reflectors
257: // * from DSYTRD;
258: // * = 'L': Lower triangle of A contains elementary reflectors
259: // * from DSYTRD.
260: // *
261: // * TRANS (input) CHARACTER*1
262: // * = 'N': No transpose, apply Q;
263: // * = 'T': Transpose, apply Q**T.
264: // *
265: // * M (input) INTEGER
266: // * The number of rows of the matrix C. M >= 0.
267: // *
268: // * N (input) INTEGER
269: // * The number of columns of the matrix C. N >= 0.
270: // *
271: // * A (input) DOUBLE PRECISION array, dimension
272: // * (LDA,M) if SIDE = 'L'
273: // * (LDA,N) if SIDE = 'R'
274: // * The vectors which define the elementary reflectors, as
275: // * returned by DSYTRD.
276: // *
277: // * LDA (input) INTEGER
278: // * The leading dimension of the array A.
279: // * LDA >= max(1,M) if SIDE = 'L'; LDA >= max(1,N) if SIDE = 'R'.
280: // *
281: // * TAU (input) DOUBLE PRECISION array, dimension
282: // * (M-1) if SIDE = 'L'
283: // * (N-1) if SIDE = 'R'
284: // * TAU(i) must contain the scalar factor of the elementary
285: // * reflector H(i), as returned by DSYTRD.
286: // *
287: // * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
288: // * On entry, the M-by-N matrix C.
289: // * On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
290: // *
291: // * LDC (input) INTEGER
292: // * The leading dimension of the array C. LDC >= max(1,M).
293: // *
294: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
295: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
296: // *
297: // * LWORK (input) INTEGER
298: // * The dimension of the array WORK.
299: // * If SIDE = 'L', LWORK >= max(1,N);
300: // * if SIDE = 'R', LWORK >= max(1,M).
301: // * For optimum performance LWORK >= N*NB if SIDE = 'L', and
302: // * LWORK >= M*NB if SIDE = 'R', where NB is the optimal
303: // * blocksize.
304: // *
305: // * If LWORK = -1, then a workspace query is assumed; the routine
306: // * only calculates the optimal size of the WORK array, returns
307: // * this value as the first entry of the WORK array, and no error
308: // * message related to LWORK is issued by XERBLA.
309: // *
310: // * INFO (output) INTEGER
311: // * = 0: successful exit
312: // * < 0: if INFO = -i, the i-th argument had an illegal value
313: // *
314: // * =====================================================================
315: // *
316: // * .. Local Scalars ..
317: // * ..
318: // * .. External Functions ..
319: // * ..
320: // * .. External Subroutines ..
321: // * ..
322: // * .. Intrinsic Functions ..
323: // INTRINSIC MAX;
324: // * ..
325: // * .. Executable Statements ..
326: // *
327: // * Test the input arguments
328: // *
329:
330: #endregion
331:
332:
333: #region Body
334:
335: INFO = 0;
336: LEFT = this._lsame.Run(SIDE, "L");
337: UPPER = this._lsame.Run(UPLO, "U");
338: LQUERY = (LWORK == - 1);
339: // *
340: // * NQ is the order of Q and NW is the minimum dimension of WORK
341: // *
342: if (LEFT)
343: {
344: NQ = M;
345: NW = N;
346: }
347: else
348: {
349: NQ = N;
350: NW = M;
351: }
352: if (!LEFT && !this._lsame.Run(SIDE, "R"))
353: {
354: INFO = - 1;
355: }
356: else
357: {
358: if (!UPPER && !this._lsame.Run(UPLO, "L"))
359: {
360: INFO = - 2;
361: }
362: else
363: {
364: if (!this._lsame.Run(TRANS, "N") && !this._lsame.Run(TRANS, "T"))
365: {
366: INFO = - 3;
367: }
368: else
369: {
370: if (M < 0)
371: {
372: INFO = - 4;
373: }
374: else
375: {
376: if (N < 0)
377: {
378: INFO = - 5;
379: }
380: else
381: {
382: if (LDA < Math.Max(1, NQ))
383: {
384: INFO = - 7;
385: }
386: else
387: {
388: if (LDC < Math.Max(1, M))
389: {
390: INFO = - 10;
391: }
392: else
393: {
394: if (LWORK < Math.Max(1, NW) && !LQUERY)
395: {
396: INFO = - 12;
397: }
398: }
399: }
400: }
401: }
402: }
403: }
404: }
405: // *
406: if (INFO == 0)
407: {
408: if (UPPER)
409: {
410: if (LEFT)
411: {
412: NB = this._ilaenv.Run(1, "DORMQL", SIDE + TRANS, M - 1, N, M - 1, - 1);
413: }
414: else
415: {
416: NB = this._ilaenv.Run(1, "DORMQL", SIDE + TRANS, M, N - 1, N - 1, - 1);
417: }
418: }
419: else
420: {
421: if (LEFT)
422: {
423: NB = this._ilaenv.Run(1, "DORMQR", SIDE + TRANS, M - 1, N, M - 1, - 1);
424: }
425: else
426: {
427: NB = this._ilaenv.Run(1, "DORMQR", SIDE + TRANS, M, N - 1, N - 1, - 1);
428: }
429: }
430: LWKOPT = Math.Max(1, NW) * NB;
431: WORK[1 + o_work] = LWKOPT;
432: }
433: // *
434: if (INFO != 0)
435: {
436: this._xerbla.Run("DORMTR", - INFO);
437: return;
438: }
439: else
440: {
441: if (LQUERY)
442: {
443: return;
444: }
445: }
446: // *
447: // * Quick return if possible
448: // *
449: if (M == 0 || N == 0 || NQ == 1)
450: {
451: WORK[1 + o_work] = 1;
452: return;
453: }
454: // *
455: if (LEFT)
456: {
457: MI = M - 1;
458: NI = N;
459: }
460: else
461: {
462: MI = M;
463: NI = N - 1;
464: }
465: // *
466: if (UPPER)
467: {
468: // *
469: // * Q was determined by a call to DSYTRD with UPLO = 'U'
470: // *
471: this._dormql.Run(SIDE, TRANS, MI, NI, NQ - 1, ref A, 1+2 * LDA + o_a
472: , LDA, TAU, offset_tau, ref C, offset_c, LDC, ref WORK, offset_work, LWORK
473: , ref IINFO);
474: }
475: else
476: {
477: // *
478: // * Q was determined by a call to DSYTRD with UPLO = 'L'
479: // *
480: if (LEFT)
481: {
482: I1 = 2;
483: I2 = 1;
484: }
485: else
486: {
487: I1 = 1;
488: I2 = 2;
489: }
490: this._dormqr.Run(SIDE, TRANS, MI, NI, NQ - 1, ref A, 2+1 * LDA + o_a
491: , LDA, TAU, offset_tau, ref C, I1+I2 * LDC + o_c, LDC, ref WORK, offset_work, LWORK
492: , ref IINFO);
493: }
494: WORK[1 + o_work] = LWKOPT;
495: return;
496: // *
497: // * End of DORMTR
498: // *
499:
500: #endregion
501:
502: }
503: }
504: }