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: /// DORGBR generates one of the real orthogonal matrices Q or P**T
27: /// determined by DGEBRD when reducing a real matrix A to bidiagonal
28: /// form: A = Q * B * P**T. Q and P**T are defined as products of
29: /// elementary reflectors H(i) or G(i) respectively.
30: ///
31: /// If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
32: /// is of order M:
33: /// if m .GE. k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
34: /// columns of Q, where m .GE. n .GE. k;
35: /// if m .LT. k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
36: /// M-by-M matrix.
37: ///
38: /// If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
39: /// is of order N:
40: /// if k .LT. n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
41: /// rows of P**T, where n .GE. m .GE. k;
42: /// if k .GE. n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
43: /// an N-by-N matrix.
44: ///
45: ///</summary>
46: public class DORGBR
47: {
48:
49:
50: #region Dependencies
51:
52: LSAME _lsame; ILAENV _ilaenv; DORGLQ _dorglq; DORGQR _dorgqr; XERBLA _xerbla;
53:
54: #endregion
55:
56:
57: #region Fields
58:
59: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; bool WANTQ = false; int I = 0; int IINFO = 0;
60: int J = 0;int LWKOPT = 0; int MN = 0; int NB = 0;
61:
62: #endregion
63:
64: public DORGBR(LSAME lsame, ILAENV ilaenv, DORGLQ dorglq, DORGQR dorgqr, XERBLA xerbla)
65: {
66:
67:
68: #region Set Dependencies
69:
70: this._lsame = lsame; this._ilaenv = ilaenv; this._dorglq = dorglq; this._dorgqr = dorgqr; this._xerbla = xerbla;
71:
72: #endregion
73:
74: }
75:
76: public DORGBR()
77: {
78:
79:
80: #region Dependencies (Initialization)
81:
82: LSAME lsame = new LSAME();
83: IEEECK ieeeck = new IEEECK();
84: IPARMQ iparmq = new IPARMQ();
85: DCOPY dcopy = new DCOPY();
86: XERBLA xerbla = new XERBLA();
87: DSCAL dscal = new DSCAL();
88: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
89: DGEMM dgemm = new DGEMM(lsame, xerbla);
90: DTRMM dtrmm = new DTRMM(lsame, xerbla);
91: DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
92: DGEMV dgemv = new DGEMV(lsame, xerbla);
93: DTRMV dtrmv = new DTRMV(lsame, xerbla);
94: DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
95: DGER dger = new DGER(xerbla);
96: DLARF dlarf = new DLARF(dgemv, dger, lsame);
97: DORGL2 dorgl2 = new DORGL2(dlarf, dscal, xerbla);
98: DORGLQ dorglq = new DORGLQ(dlarfb, dlarft, dorgl2, xerbla, ilaenv);
99: DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
100: DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
101:
102: #endregion
103:
104:
105: #region Set Dependencies
106:
107: this._lsame = lsame; this._ilaenv = ilaenv; this._dorglq = dorglq; this._dorgqr = dorgqr; this._xerbla = xerbla;
108:
109: #endregion
110:
111: }
112: /// <summary>
113: /// Purpose
114: /// =======
115: ///
116: /// DORGBR generates one of the real orthogonal matrices Q or P**T
117: /// determined by DGEBRD when reducing a real matrix A to bidiagonal
118: /// form: A = Q * B * P**T. Q and P**T are defined as products of
119: /// elementary reflectors H(i) or G(i) respectively.
120: ///
121: /// If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
122: /// is of order M:
123: /// if m .GE. k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
124: /// columns of Q, where m .GE. n .GE. k;
125: /// if m .LT. k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
126: /// M-by-M matrix.
127: ///
128: /// If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
129: /// is of order N:
130: /// if k .LT. n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
131: /// rows of P**T, where n .GE. m .GE. k;
132: /// if k .GE. n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
133: /// an N-by-N matrix.
134: ///
135: ///</summary>
136: /// <param name="VECT">
137: /// (input) CHARACTER*1
138: /// Specifies whether the matrix Q or the matrix P**T is
139: /// required, as defined in the transformation applied by DGEBRD:
140: /// = 'Q': generate Q;
141: /// = 'P': generate P**T.
142: ///</param>
143: /// <param name="M">
144: /// (input) INTEGER
145: /// The number of rows of the matrix Q or P**T to be returned.
146: /// M .GE. 0.
147: ///</param>
148: /// <param name="N">
149: /// (input) INTEGER
150: /// The number of columns of the matrix Q or P**T to be returned.
151: /// N .GE. 0.
152: /// If VECT = 'Q', M .GE. N .GE. min(M,K);
153: /// if VECT = 'P', N .GE. M .GE. min(N,K).
154: ///</param>
155: /// <param name="K">
156: /// (input) INTEGER
157: /// If VECT = 'Q', the number of columns in the original M-by-K
158: /// matrix reduced by DGEBRD.
159: /// If VECT = 'P', the number of rows in the original K-by-N
160: /// matrix reduced by DGEBRD.
161: /// K .GE. 0.
162: ///</param>
163: /// <param name="A">
164: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
165: /// On entry, the vectors which define the elementary reflectors,
166: /// as returned by DGEBRD.
167: /// On exit, the M-by-N matrix Q or P**T.
168: ///</param>
169: /// <param name="LDA">
170: /// (input) INTEGER
171: /// The leading dimension of the array A. LDA .GE. max(1,M).
172: ///</param>
173: /// <param name="TAU">
174: /// (input) DOUBLE PRECISION array, dimension
175: /// (min(M,K)) if VECT = 'Q'
176: /// (min(N,K)) if VECT = 'P'
177: /// TAU(i) must contain the scalar factor of the elementary
178: /// reflector H(i) or G(i), which determines Q or P**T, as
179: /// returned by DGEBRD in its array argument TAUQ or TAUP.
180: ///</param>
181: /// <param name="WORK">
182: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
183: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
184: ///</param>
185: /// <param name="LWORK">
186: /// (input) INTEGER
187: /// The dimension of the array WORK. LWORK .GE. max(1,min(M,N)).
188: /// For optimum performance LWORK .GE. min(M,N)*NB, where NB
189: /// is the optimal blocksize.
190: ///
191: /// If LWORK = -1, then a workspace query is assumed; the routine
192: /// only calculates the optimal size of the WORK array, returns
193: /// this value as the first entry of the WORK array, and no error
194: /// message related to LWORK is issued by XERBLA.
195: ///</param>
196: /// <param name="INFO">
197: /// (output) INTEGER
198: /// = 0: successful exit
199: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
200: ///</param>
201: public void Run(string VECT, int M, int N, int K, ref double[] A, int offset_a, int LDA
202: , double[] TAU, int offset_tau, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
203: {
204:
205: #region Array Index Correction
206:
207: int o_a = -1 - LDA + offset_a; int o_tau = -1 + offset_tau; int o_work = -1 + offset_work;
208:
209: #endregion
210:
211:
212: #region Strings
213:
214: VECT = VECT.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: // * DORGBR generates one of the real orthogonal matrices Q or P**T
235: // * determined by DGEBRD when reducing a real matrix A to bidiagonal
236: // * form: A = Q * B * P**T. Q and P**T are defined as products of
237: // * elementary reflectors H(i) or G(i) respectively.
238: // *
239: // * If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
240: // * is of order M:
241: // * if m >= k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
242: // * columns of Q, where m >= n >= k;
243: // * if m < k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
244: // * M-by-M matrix.
245: // *
246: // * If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
247: // * is of order N:
248: // * if k < n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
249: // * rows of P**T, where n >= m >= k;
250: // * if k >= n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
251: // * an N-by-N matrix.
252: // *
253: // * Arguments
254: // * =========
255: // *
256: // * VECT (input) CHARACTER*1
257: // * Specifies whether the matrix Q or the matrix P**T is
258: // * required, as defined in the transformation applied by DGEBRD:
259: // * = 'Q': generate Q;
260: // * = 'P': generate P**T.
261: // *
262: // * M (input) INTEGER
263: // * The number of rows of the matrix Q or P**T to be returned.
264: // * M >= 0.
265: // *
266: // * N (input) INTEGER
267: // * The number of columns of the matrix Q or P**T to be returned.
268: // * N >= 0.
269: // * If VECT = 'Q', M >= N >= min(M,K);
270: // * if VECT = 'P', N >= M >= min(N,K).
271: // *
272: // * K (input) INTEGER
273: // * If VECT = 'Q', the number of columns in the original M-by-K
274: // * matrix reduced by DGEBRD.
275: // * If VECT = 'P', the number of rows in the original K-by-N
276: // * matrix reduced by DGEBRD.
277: // * K >= 0.
278: // *
279: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
280: // * On entry, the vectors which define the elementary reflectors,
281: // * as returned by DGEBRD.
282: // * On exit, the M-by-N matrix Q or P**T.
283: // *
284: // * LDA (input) INTEGER
285: // * The leading dimension of the array A. LDA >= max(1,M).
286: // *
287: // * TAU (input) DOUBLE PRECISION array, dimension
288: // * (min(M,K)) if VECT = 'Q'
289: // * (min(N,K)) if VECT = 'P'
290: // * TAU(i) must contain the scalar factor of the elementary
291: // * reflector H(i) or G(i), which determines Q or P**T, as
292: // * returned by DGEBRD in its array argument TAUQ or TAUP.
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. LWORK >= max(1,min(M,N)).
299: // * For optimum performance LWORK >= min(M,N)*NB, where NB
300: // * is the optimal blocksize.
301: // *
302: // * If LWORK = -1, then a workspace query is assumed; the routine
303: // * only calculates the optimal size of the WORK array, returns
304: // * this value as the first entry of the WORK array, and no error
305: // * message related to LWORK is issued by XERBLA.
306: // *
307: // * INFO (output) INTEGER
308: // * = 0: successful exit
309: // * < 0: if INFO = -i, the i-th argument had an illegal value
310: // *
311: // * =====================================================================
312: // *
313: // * .. Parameters ..
314: // * ..
315: // * .. Local Scalars ..
316: // * ..
317: // * .. External Functions ..
318: // * ..
319: // * .. External Subroutines ..
320: // * ..
321: // * .. Intrinsic Functions ..
322: // INTRINSIC MAX, MIN;
323: // * ..
324: // * .. Executable Statements ..
325: // *
326: // * Test the input arguments
327: // *
328:
329: #endregion
330:
331:
332: #region Body
333:
334: INFO = 0;
335: WANTQ = this._lsame.Run(VECT, "Q");
336: MN = Math.Min(M, N);
337: LQUERY = (LWORK == - 1);
338: if (!WANTQ && !this._lsame.Run(VECT, "P"))
339: {
340: INFO = - 1;
341: }
342: else
343: {
344: if (M < 0)
345: {
346: INFO = - 2;
347: }
348: else
349: {
350: if (N < 0 || (WANTQ && (N > M || N < Math.Min(M, K))) || (!WANTQ && (M > N || M < Math.Min(N, K))))
351: {
352: INFO = - 3;
353: }
354: else
355: {
356: if (K < 0)
357: {
358: INFO = - 4;
359: }
360: else
361: {
362: if (LDA < Math.Max(1, M))
363: {
364: INFO = - 6;
365: }
366: else
367: {
368: if (LWORK < Math.Max(1, MN) && !LQUERY)
369: {
370: INFO = - 9;
371: }
372: }
373: }
374: }
375: }
376: }
377: // *
378: if (INFO == 0)
379: {
380: if (WANTQ)
381: {
382: NB = this._ilaenv.Run(1, "DORGQR", " ", M, N, K, - 1);
383: }
384: else
385: {
386: NB = this._ilaenv.Run(1, "DORGLQ", " ", M, N, K, - 1);
387: }
388: LWKOPT = Math.Max(1, MN) * NB;
389: WORK[1 + o_work] = LWKOPT;
390: }
391: // *
392: if (INFO != 0)
393: {
394: this._xerbla.Run("DORGBR", - INFO);
395: return;
396: }
397: else
398: {
399: if (LQUERY)
400: {
401: return;
402: }
403: }
404: // *
405: // * Quick return if possible
406: // *
407: if (M == 0 || N == 0)
408: {
409: WORK[1 + o_work] = 1;
410: return;
411: }
412: // *
413: if (WANTQ)
414: {
415: // *
416: // * Form Q, determined by a call to DGEBRD to reduce an m-by-k
417: // * matrix
418: // *
419: if (M >= K)
420: {
421: // *
422: // * If m >= k, assume m >= n >= k
423: // *
424: this._dorgqr.Run(M, N, K, ref A, offset_a, LDA, TAU, offset_tau
425: , ref WORK, offset_work, LWORK, ref IINFO);
426: // *
427: }
428: else
429: {
430: // *
431: // * If m < k, assume m = n
432: // *
433: // * Shift the vectors which define the elementary reflectors one
434: // * column to the right, and set the first row and column of Q
435: // * to those of the unit matrix
436: // *
437: for (J = M; J >= 2; J += - 1)
438: {
439: A[1+J * LDA + o_a] = ZERO;
440: for (I = J + 1; I <= M; I++)
441: {
442: A[I+J * LDA + o_a] = A[I+(J - 1) * LDA + o_a];
443: }
444: }
445: A[1+1 * LDA + o_a] = ONE;
446: for (I = 2; I <= M; I++)
447: {
448: A[I+1 * LDA + o_a] = ZERO;
449: }
450: if (M > 1)
451: {
452: // *
453: // * Form Q(2:m,2:m)
454: // *
455: this._dorgqr.Run(M - 1, M - 1, M - 1, ref A, 2+2 * LDA + o_a, LDA, TAU, offset_tau
456: , ref WORK, offset_work, LWORK, ref IINFO);
457: }
458: }
459: }
460: else
461: {
462: // *
463: // * Form P', determined by a call to DGEBRD to reduce a k-by-n
464: // * matrix
465: // *
466: if (K < N)
467: {
468: // *
469: // * If k < n, assume k <= m <= n
470: // *
471: this._dorglq.Run(M, N, K, ref A, offset_a, LDA, TAU, offset_tau
472: , ref WORK, offset_work, LWORK, ref IINFO);
473: // *
474: }
475: else
476: {
477: // *
478: // * If k >= n, assume m = n
479: // *
480: // * Shift the vectors which define the elementary reflectors one
481: // * row downward, and set the first row and column of P' to
482: // * those of the unit matrix
483: // *
484: A[1+1 * LDA + o_a] = ONE;
485: for (I = 2; I <= N; I++)
486: {
487: A[I+1 * LDA + o_a] = ZERO;
488: }
489: for (J = 2; J <= N; J++)
490: {
491: for (I = J - 1; I >= 2; I += - 1)
492: {
493: A[I+J * LDA + o_a] = A[I - 1+J * LDA + o_a];
494: }
495: A[1+J * LDA + o_a] = ZERO;
496: }
497: if (N > 1)
498: {
499: // *
500: // * Form P'(2:n,2:n)
501: // *
502: this._dorglq.Run(N - 1, N - 1, N - 1, ref A, 2+2 * LDA + o_a, LDA, TAU, offset_tau
503: , ref WORK, offset_work, LWORK, ref IINFO);
504: }
505: }
506: }
507: WORK[1 + o_work] = LWKOPT;
508: return;
509: // *
510: // * End of DORGBR
511: // *
512:
513: #endregion
514:
515: }
516: }
517: }