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: /// DGEBRD reduces a general real M-by-N matrix A to upper or lower
27: /// bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
28: ///
29: /// If m .GE. n, B is upper bidiagonal; if m .LT. n, B is lower bidiagonal.
30: ///
31: ///</summary>
32: public class DGEBRD
33: {
34:
35:
36: #region Dependencies
37:
38: DGEBD2 _dgebd2; DGEMM _dgemm; DLABRD _dlabrd; XERBLA _xerbla; ILAENV _ilaenv;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ONE = 1.0E+0; bool LQUERY = false; int I = 0; int IINFO = 0; int J = 0; int LDWRKX = 0; int LDWRKY = 0;
46: int LWKOPT = 0;int MINMN = 0; int NB = 0; int NBMIN = 0; int NX = 0; double WS = 0;
47:
48: #endregion
49:
50: public DGEBRD(DGEBD2 dgebd2, DGEMM dgemm, DLABRD dlabrd, XERBLA xerbla, ILAENV ilaenv)
51: {
52:
53:
54: #region Set Dependencies
55:
56: this._dgebd2 = dgebd2; this._dgemm = dgemm; this._dlabrd = dlabrd; this._xerbla = xerbla; this._ilaenv = ilaenv;
57:
58: #endregion
59:
60: }
61:
62: public DGEBRD()
63: {
64:
65:
66: #region Dependencies (Initialization)
67:
68: LSAME lsame = new LSAME();
69: XERBLA xerbla = new XERBLA();
70: DLAMC3 dlamc3 = new DLAMC3();
71: DLAPY2 dlapy2 = new DLAPY2();
72: DNRM2 dnrm2 = new DNRM2();
73: DSCAL dscal = new DSCAL();
74: IEEECK ieeeck = new IEEECK();
75: IPARMQ iparmq = new IPARMQ();
76: DGEMV dgemv = new DGEMV(lsame, xerbla);
77: DGER dger = new DGER(xerbla);
78: DLARF dlarf = new DLARF(dgemv, dger, lsame);
79: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
80: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
81: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
82: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
83: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
84: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
85: DGEBD2 dgebd2 = new DGEBD2(dlarf, dlarfg, xerbla);
86: DGEMM dgemm = new DGEMM(lsame, xerbla);
87: DLABRD dlabrd = new DLABRD(dgemv, dlarfg, dscal);
88: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
89:
90: #endregion
91:
92:
93: #region Set Dependencies
94:
95: this._dgebd2 = dgebd2; this._dgemm = dgemm; this._dlabrd = dlabrd; this._xerbla = xerbla; this._ilaenv = ilaenv;
96:
97: #endregion
98:
99: }
100: /// <summary>
101: /// Purpose
102: /// =======
103: ///
104: /// DGEBRD reduces a general real M-by-N matrix A to upper or lower
105: /// bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
106: ///
107: /// If m .GE. n, B is upper bidiagonal; if m .LT. n, B is lower bidiagonal.
108: ///
109: ///</summary>
110: /// <param name="M">
111: /// (input) INTEGER
112: /// The number of rows in the matrix A. M .GE. 0.
113: ///</param>
114: /// <param name="N">
115: /// (input) INTEGER
116: /// The number of columns in the matrix A. N .GE. 0.
117: ///</param>
118: /// <param name="A">
119: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
120: /// On entry, the M-by-N general matrix to be reduced.
121: /// On exit,
122: /// if m .GE. n, the diagonal and the first superdiagonal are
123: /// overwritten with the upper bidiagonal matrix B; the
124: /// elements below the diagonal, with the array TAUQ, represent
125: /// the orthogonal matrix Q as a product of elementary
126: /// reflectors, and the elements above the first superdiagonal,
127: /// with the array TAUP, represent the orthogonal matrix P as
128: /// a product of elementary reflectors;
129: /// if m .LT. n, the diagonal and the first subdiagonal are
130: /// overwritten with the lower bidiagonal matrix B; the
131: /// elements below the first subdiagonal, with the array TAUQ,
132: /// represent the orthogonal matrix Q as a product of
133: /// elementary reflectors, and the elements above the diagonal,
134: /// with the array TAUP, represent the orthogonal matrix P as
135: /// a product of elementary reflectors.
136: /// See Further Details.
137: ///</param>
138: /// <param name="LDA">
139: /// (input) INTEGER
140: /// The leading dimension of the array A. LDA .GE. max(1,M).
141: ///</param>
142: /// <param name="D">
143: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
144: /// The diagonal elements of the bidiagonal matrix B:
145: /// D(i) = A(i,i).
146: ///</param>
147: /// <param name="E">
148: /// (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
149: /// The off-diagonal elements of the bidiagonal matrix B:
150: /// if m .GE. n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
151: /// if m .LT. n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
152: ///</param>
153: /// <param name="TAUQ">
154: /// (output) DOUBLE PRECISION array dimension (min(M,N))
155: /// The scalar factors of the elementary reflectors which
156: /// represent the orthogonal matrix Q. See Further Details.
157: ///</param>
158: /// <param name="TAUP">
159: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
160: /// The scalar factors of the elementary reflectors which
161: /// represent the orthogonal matrix P. See Further Details.
162: ///</param>
163: /// <param name="WORK">
164: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
165: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
166: ///</param>
167: /// <param name="LWORK">
168: /// (input) INTEGER
169: /// The length of the array WORK. LWORK .GE. max(1,M,N).
170: /// For optimum performance LWORK .GE. (M+N)*NB, where NB
171: /// is the optimal blocksize.
172: ///
173: /// If LWORK = -1, then a workspace query is assumed; the routine
174: /// only calculates the optimal size of the WORK array, returns
175: /// this value as the first entry of the WORK array, and no error
176: /// message related to LWORK is issued by XERBLA.
177: ///</param>
178: /// <param name="INFO">
179: /// (output) INTEGER
180: /// = 0: successful exit
181: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
182: ///</param>
183: public void Run(int M, int N, ref double[] A, int offset_a, int LDA, ref double[] D, int offset_d, ref double[] E, int offset_e
184: , ref double[] TAUQ, int offset_tauq, ref double[] TAUP, int offset_taup, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
185: {
186:
187: #region Array Index Correction
188:
189: int o_a = -1 - LDA + offset_a; int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_tauq = -1 + offset_tauq;
190: int o_taup = -1 + offset_taup; int o_work = -1 + offset_work;
191:
192: #endregion
193:
194:
195: #region Prolog
196:
197: // *
198: // * -- LAPACK routine (version 3.1) --
199: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
200: // * November 2006
201: // *
202: // * .. Scalar Arguments ..
203: // * ..
204: // * .. Array Arguments ..
205: // * ..
206: // *
207: // * Purpose
208: // * =======
209: // *
210: // * DGEBRD reduces a general real M-by-N matrix A to upper or lower
211: // * bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
212: // *
213: // * If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
214: // *
215: // * Arguments
216: // * =========
217: // *
218: // * M (input) INTEGER
219: // * The number of rows in the matrix A. M >= 0.
220: // *
221: // * N (input) INTEGER
222: // * The number of columns in the matrix A. N >= 0.
223: // *
224: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
225: // * On entry, the M-by-N general matrix to be reduced.
226: // * On exit,
227: // * if m >= n, the diagonal and the first superdiagonal are
228: // * overwritten with the upper bidiagonal matrix B; the
229: // * elements below the diagonal, with the array TAUQ, represent
230: // * the orthogonal matrix Q as a product of elementary
231: // * reflectors, and the elements above the first superdiagonal,
232: // * with the array TAUP, represent the orthogonal matrix P as
233: // * a product of elementary reflectors;
234: // * if m < n, the diagonal and the first subdiagonal are
235: // * overwritten with the lower bidiagonal matrix B; the
236: // * elements below the first subdiagonal, with the array TAUQ,
237: // * represent the orthogonal matrix Q as a product of
238: // * elementary reflectors, and the elements above the diagonal,
239: // * with the array TAUP, represent the orthogonal matrix P as
240: // * a product of elementary reflectors.
241: // * See Further Details.
242: // *
243: // * LDA (input) INTEGER
244: // * The leading dimension of the array A. LDA >= max(1,M).
245: // *
246: // * D (output) DOUBLE PRECISION array, dimension (min(M,N))
247: // * The diagonal elements of the bidiagonal matrix B:
248: // * D(i) = A(i,i).
249: // *
250: // * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
251: // * The off-diagonal elements of the bidiagonal matrix B:
252: // * if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
253: // * if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
254: // *
255: // * TAUQ (output) DOUBLE PRECISION array dimension (min(M,N))
256: // * The scalar factors of the elementary reflectors which
257: // * represent the orthogonal matrix Q. See Further Details.
258: // *
259: // * TAUP (output) DOUBLE PRECISION array, dimension (min(M,N))
260: // * The scalar factors of the elementary reflectors which
261: // * represent the orthogonal matrix P. See Further Details.
262: // *
263: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
264: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
265: // *
266: // * LWORK (input) INTEGER
267: // * The length of the array WORK. LWORK >= max(1,M,N).
268: // * For optimum performance LWORK >= (M+N)*NB, where NB
269: // * is the optimal blocksize.
270: // *
271: // * If LWORK = -1, then a workspace query is assumed; the routine
272: // * only calculates the optimal size of the WORK array, returns
273: // * this value as the first entry of the WORK array, and no error
274: // * message related to LWORK is issued by XERBLA.
275: // *
276: // * INFO (output) INTEGER
277: // * = 0: successful exit
278: // * < 0: if INFO = -i, the i-th argument had an illegal value.
279: // *
280: // * Further Details
281: // * ===============
282: // *
283: // * The matrices Q and P are represented as products of elementary
284: // * reflectors:
285: // *
286: // * If m >= n,
287: // *
288: // * Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1)
289: // *
290: // * Each H(i) and G(i) has the form:
291: // *
292: // * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
293: // *
294: // * where tauq and taup are real scalars, and v and u are real vectors;
295: // * v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
296: // * u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
297: // * tauq is stored in TAUQ(i) and taup in TAUP(i).
298: // *
299: // * If m < n,
300: // *
301: // * Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m)
302: // *
303: // * Each H(i) and G(i) has the form:
304: // *
305: // * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
306: // *
307: // * where tauq and taup are real scalars, and v and u are real vectors;
308: // * v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
309: // * u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
310: // * tauq is stored in TAUQ(i) and taup in TAUP(i).
311: // *
312: // * The contents of A on exit are illustrated by the following examples:
313: // *
314: // * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
315: // *
316: // * ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
317: // * ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
318: // * ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
319: // * ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
320: // * ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
321: // * ( v1 v2 v3 v4 v5 )
322: // *
323: // * where d and e denote diagonal and off-diagonal elements of B, vi
324: // * denotes an element of the vector defining H(i), and ui an element of
325: // * the vector defining G(i).
326: // *
327: // * =====================================================================
328: // *
329: // * .. Parameters ..
330: // * ..
331: // * .. Local Scalars ..
332: // * ..
333: // * .. External Subroutines ..
334: // * ..
335: // * .. Intrinsic Functions ..
336: // INTRINSIC DBLE, MAX, MIN;
337: // * ..
338: // * .. External Functions ..
339: // * ..
340: // * .. Executable Statements ..
341: // *
342: // * Test the input parameters
343: // *
344:
345: #endregion
346:
347:
348: #region Body
349:
350: INFO = 0;
351: NB = Math.Max(1, this._ilaenv.Run(1, "DGEBRD", " ", M, N, - 1, - 1));
352: LWKOPT = (M + N) * NB;
353: WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
354: LQUERY = (LWORK == - 1);
355: if (M < 0)
356: {
357: INFO = - 1;
358: }
359: else
360: {
361: if (N < 0)
362: {
363: INFO = - 2;
364: }
365: else
366: {
367: if (LDA < Math.Max(1, M))
368: {
369: INFO = - 4;
370: }
371: else
372: {
373: if (LWORK < Math.Max(1, Math.Max(M, N)) && !LQUERY)
374: {
375: INFO = - 10;
376: }
377: }
378: }
379: }
380: if (INFO < 0)
381: {
382: this._xerbla.Run("DGEBRD", - INFO);
383: return;
384: }
385: else
386: {
387: if (LQUERY)
388: {
389: return;
390: }
391: }
392: // *
393: // * Quick return if possible
394: // *
395: MINMN = Math.Min(M, N);
396: if (MINMN == 0)
397: {
398: WORK[1 + o_work] = 1;
399: return;
400: }
401: // *
402: WS = Math.Max(M, N);
403: LDWRKX = M;
404: LDWRKY = N;
405: // *
406: if (NB > 1 && NB < MINMN)
407: {
408: // *
409: // * Set the crossover point NX.
410: // *
411: NX = Math.Max(NB, this._ilaenv.Run(3, "DGEBRD", " ", M, N, - 1, - 1));
412: // *
413: // * Determine when to switch from blocked to unblocked code.
414: // *
415: if (NX < MINMN)
416: {
417: WS = (M + N) * NB;
418: if (LWORK < WS)
419: {
420: // *
421: // * Not enough work space for the optimal NB, consider using
422: // * a smaller block size.
423: // *
424: NBMIN = this._ilaenv.Run(2, "DGEBRD", " ", M, N, - 1, - 1);
425: if (LWORK >= (M + N) * NBMIN)
426: {
427: NB = LWORK / (M + N);
428: }
429: else
430: {
431: NB = 1;
432: NX = MINMN;
433: }
434: }
435: }
436: }
437: else
438: {
439: NX = MINMN;
440: }
441: // *
442: for (I = 1; (NB >= 0) ? (I <= MINMN - NX) : (I >= MINMN - NX); I += NB)
443: {
444: // *
445: // * Reduce rows and columns i:i+nb-1 to bidiagonal form and return
446: // * the matrices X and Y which are needed to update the unreduced
447: // * part of the matrix
448: // *
449: this._dlabrd.Run(M - I + 1, N - I + 1, NB, ref A, I+I * LDA + o_a, LDA, ref D, I + o_d
450: , ref E, I + o_e, ref TAUQ, I + o_tauq, ref TAUP, I + o_taup, ref WORK, offset_work, LDWRKX, ref WORK, LDWRKX * NB + 1 + o_work
451: , LDWRKY);
452: // *
453: // * Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
454: // * of the form A := A - V*Y' - X*U'
455: // *
456: this._dgemm.Run("No transpose", "Transpose", M - I - NB + 1, N - I - NB + 1, NB, - ONE
457: , A, I + NB+I * LDA + o_a, LDA, WORK, LDWRKX * NB + NB + 1 + o_work, LDWRKY, ONE, ref A, I + NB+(I + NB) * LDA + o_a
458: , LDA);
459: this._dgemm.Run("No transpose", "No transpose", M - I - NB + 1, N - I - NB + 1, NB, - ONE
460: , WORK, NB + 1 + o_work, LDWRKX, A, I+(I + NB) * LDA + o_a, LDA, ONE, ref A, I + NB+(I + NB) * LDA + o_a
461: , LDA);
462: // *
463: // * Copy diagonal and off-diagonal elements of B back into A
464: // *
465: if (M >= N)
466: {
467: for (J = I; J <= I + NB - 1; J++)
468: {
469: A[J+J * LDA + o_a] = D[J + o_d];
470: A[J+(J + 1) * LDA + o_a] = E[J + o_e];
471: }
472: }
473: else
474: {
475: for (J = I; J <= I + NB - 1; J++)
476: {
477: A[J+J * LDA + o_a] = D[J + o_d];
478: A[J + 1+J * LDA + o_a] = E[J + o_e];
479: }
480: }
481: }
482: // *
483: // * Use unblocked code to reduce the remainder of the matrix
484: // *
485: this._dgebd2.Run(M - I + 1, N - I + 1, ref A, I+I * LDA + o_a, LDA, ref D, I + o_d, ref E, I + o_e
486: , ref TAUQ, I + o_tauq, ref TAUP, I + o_taup, ref WORK, offset_work, ref IINFO);
487: WORK[1 + o_work] = WS;
488: return;
489: // *
490: // * End of DGEBRD
491: // *
492:
493: #endregion
494:
495: }
496: }
497: }