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