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: /// DGEBD2 reduces a real general m by n matrix A to upper or lower
27: /// bidiagonal form B by an orthogonal transformation: Q' * 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 DGEBD2
33: {
34:
35:
36: #region Dependencies
37:
38: DLARF _dlarf; DLARFG _dlarfg; XERBLA _xerbla;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; int I = 0;
46:
47: #endregion
48:
49: public DGEBD2(DLARF dlarf, DLARFG dlarfg, XERBLA xerbla)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla;
56:
57: #endregion
58:
59: }
60:
61: public DGEBD2()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: LSAME lsame = new LSAME();
68: XERBLA xerbla = new XERBLA();
69: DLAMC3 dlamc3 = new DLAMC3();
70: DLAPY2 dlapy2 = new DLAPY2();
71: DNRM2 dnrm2 = new DNRM2();
72: DSCAL dscal = new DSCAL();
73: DGEMV dgemv = new DGEMV(lsame, xerbla);
74: DGER dger = new DGER(xerbla);
75: DLARF dlarf = new DLARF(dgemv, dger, lsame);
76: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
77: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
78: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
79: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
80: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
81: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
82:
83: #endregion
84:
85:
86: #region Set Dependencies
87:
88: this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla;
89:
90: #endregion
91:
92: }
93: /// <summary>
94: /// Purpose
95: /// =======
96: ///
97: /// DGEBD2 reduces a real general m by n matrix A to upper or lower
98: /// bidiagonal form B by an orthogonal transformation: Q' * A * P = B.
99: ///
100: /// If m .GE. n, B is upper bidiagonal; if m .LT. n, B is lower bidiagonal.
101: ///
102: ///</summary>
103: /// <param name="M">
104: /// (input) INTEGER
105: /// The number of rows in the matrix A. M .GE. 0.
106: ///</param>
107: /// <param name="N">
108: /// (input) INTEGER
109: /// The number of columns in the matrix A. N .GE. 0.
110: ///</param>
111: /// <param name="A">
112: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
113: /// On entry, the m by n general matrix to be reduced.
114: /// On exit,
115: /// if m .GE. n, the diagonal and the first superdiagonal are
116: /// overwritten with the upper bidiagonal matrix B; the
117: /// elements below the diagonal, with the array TAUQ, represent
118: /// the orthogonal matrix Q as a product of elementary
119: /// reflectors, and the elements above the first superdiagonal,
120: /// with the array TAUP, represent the orthogonal matrix P as
121: /// a product of elementary reflectors;
122: /// if m .LT. n, the diagonal and the first subdiagonal are
123: /// overwritten with the lower bidiagonal matrix B; the
124: /// elements below the first subdiagonal, with the array TAUQ,
125: /// represent the orthogonal matrix Q as a product of
126: /// elementary reflectors, and the elements above the diagonal,
127: /// with the array TAUP, represent the orthogonal matrix P as
128: /// a product of elementary reflectors.
129: /// See Further Details.
130: ///</param>
131: /// <param name="LDA">
132: /// (input) INTEGER
133: /// The leading dimension of the array A. LDA .GE. max(1,M).
134: ///</param>
135: /// <param name="D">
136: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
137: /// The diagonal elements of the bidiagonal matrix B:
138: /// D(i) = A(i,i).
139: ///</param>
140: /// <param name="E">
141: /// (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
142: /// The off-diagonal elements of the bidiagonal matrix B:
143: /// if m .GE. n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
144: /// if m .LT. n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
145: ///</param>
146: /// <param name="TAUQ">
147: /// (output) DOUBLE PRECISION array dimension (min(M,N))
148: /// The scalar factors of the elementary reflectors which
149: /// represent the orthogonal matrix Q. See Further Details.
150: ///</param>
151: /// <param name="TAUP">
152: /// (output) DOUBLE PRECISION array, dimension (min(M,N))
153: /// The scalar factors of the elementary reflectors which
154: /// represent the orthogonal matrix P. See Further Details.
155: ///</param>
156: /// <param name="WORK">
157: /// (workspace) DOUBLE PRECISION array, dimension (max(M,N))
158: ///</param>
159: /// <param name="INFO">
160: /// (output) INTEGER
161: /// = 0: successful exit.
162: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
163: ///</param>
164: 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
165: , ref double[] TAUQ, int offset_tauq, ref double[] TAUP, int offset_taup, ref double[] WORK, int offset_work, ref int INFO)
166: {
167:
168: #region Array Index Correction
169:
170: 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;
171: int o_taup = -1 + offset_taup; int o_work = -1 + offset_work;
172:
173: #endregion
174:
175:
176: #region Prolog
177:
178: // *
179: // * -- LAPACK routine (version 3.1) --
180: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
181: // * November 2006
182: // *
183: // * .. Scalar Arguments ..
184: // * ..
185: // * .. Array Arguments ..
186: // * ..
187: // *
188: // * Purpose
189: // * =======
190: // *
191: // * DGEBD2 reduces a real general m by n matrix A to upper or lower
192: // * bidiagonal form B by an orthogonal transformation: Q' * A * P = B.
193: // *
194: // * If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
195: // *
196: // * Arguments
197: // * =========
198: // *
199: // * M (input) INTEGER
200: // * The number of rows in the matrix A. M >= 0.
201: // *
202: // * N (input) INTEGER
203: // * The number of columns in the matrix A. N >= 0.
204: // *
205: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
206: // * On entry, the m by n general matrix to be reduced.
207: // * On exit,
208: // * if m >= n, the diagonal and the first superdiagonal are
209: // * overwritten with the upper bidiagonal matrix B; the
210: // * elements below the diagonal, with the array TAUQ, represent
211: // * the orthogonal matrix Q as a product of elementary
212: // * reflectors, and the elements above the first superdiagonal,
213: // * with the array TAUP, represent the orthogonal matrix P as
214: // * a product of elementary reflectors;
215: // * if m < n, the diagonal and the first subdiagonal are
216: // * overwritten with the lower bidiagonal matrix B; the
217: // * elements below the first subdiagonal, with the array TAUQ,
218: // * represent the orthogonal matrix Q as a product of
219: // * elementary reflectors, and the elements above the diagonal,
220: // * with the array TAUP, represent the orthogonal matrix P as
221: // * a product of elementary reflectors.
222: // * See Further Details.
223: // *
224: // * LDA (input) INTEGER
225: // * The leading dimension of the array A. LDA >= max(1,M).
226: // *
227: // * D (output) DOUBLE PRECISION array, dimension (min(M,N))
228: // * The diagonal elements of the bidiagonal matrix B:
229: // * D(i) = A(i,i).
230: // *
231: // * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
232: // * The off-diagonal elements of the bidiagonal matrix B:
233: // * if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
234: // * if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
235: // *
236: // * TAUQ (output) DOUBLE PRECISION array dimension (min(M,N))
237: // * The scalar factors of the elementary reflectors which
238: // * represent the orthogonal matrix Q. See Further Details.
239: // *
240: // * TAUP (output) DOUBLE PRECISION array, dimension (min(M,N))
241: // * The scalar factors of the elementary reflectors which
242: // * represent the orthogonal matrix P. See Further Details.
243: // *
244: // * WORK (workspace) DOUBLE PRECISION array, dimension (max(M,N))
245: // *
246: // * INFO (output) INTEGER
247: // * = 0: successful exit.
248: // * < 0: if INFO = -i, the i-th argument had an illegal value.
249: // *
250: // * Further Details
251: // * ===============
252: // *
253: // * The matrices Q and P are represented as products of elementary
254: // * reflectors:
255: // *
256: // * If m >= n,
257: // *
258: // * Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1)
259: // *
260: // * Each H(i) and G(i) has the form:
261: // *
262: // * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
263: // *
264: // * where tauq and taup are real scalars, and v and u are real vectors;
265: // * v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
266: // * u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
267: // * tauq is stored in TAUQ(i) and taup in TAUP(i).
268: // *
269: // * If m < n,
270: // *
271: // * Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m)
272: // *
273: // * Each H(i) and G(i) has the form:
274: // *
275: // * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
276: // *
277: // * where tauq and taup are real scalars, and v and u are real vectors;
278: // * v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
279: // * u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
280: // * tauq is stored in TAUQ(i) and taup in TAUP(i).
281: // *
282: // * The contents of A on exit are illustrated by the following examples:
283: // *
284: // * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
285: // *
286: // * ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
287: // * ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
288: // * ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
289: // * ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
290: // * ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
291: // * ( v1 v2 v3 v4 v5 )
292: // *
293: // * where d and e denote diagonal and off-diagonal elements of B, vi
294: // * denotes an element of the vector defining H(i), and ui an element of
295: // * the vector defining G(i).
296: // *
297: // * =====================================================================
298: // *
299: // * .. Parameters ..
300: // * ..
301: // * .. Local Scalars ..
302: // * ..
303: // * .. External Subroutines ..
304: // * ..
305: // * .. Intrinsic Functions ..
306: // INTRINSIC MAX, MIN;
307: // * ..
308: // * .. Executable Statements ..
309: // *
310: // * Test the input parameters
311: // *
312:
313: #endregion
314:
315:
316: #region Body
317:
318: INFO = 0;
319: if (M < 0)
320: {
321: INFO = - 1;
322: }
323: else
324: {
325: if (N < 0)
326: {
327: INFO = - 2;
328: }
329: else
330: {
331: if (LDA < Math.Max(1, M))
332: {
333: INFO = - 4;
334: }
335: }
336: }
337: if (INFO < 0)
338: {
339: this._xerbla.Run("DGEBD2", - INFO);
340: return;
341: }
342: // *
343: if (M >= N)
344: {
345: // *
346: // * Reduce to upper bidiagonal form
347: // *
348: for (I = 1; I <= N; I++)
349: {
350: // *
351: // * Generate elementary reflector H(i) to annihilate A(i+1:m,i)
352: // *
353: this._dlarfg.Run(M - I + 1, ref A[I+I * LDA + o_a], ref A, Math.Min(I + 1, M)+I * LDA + o_a, 1, ref TAUQ[I + o_tauq]);
354: D[I + o_d] = A[I+I * LDA + o_a];
355: A[I+I * LDA + o_a] = ONE;
356: // *
357: // * Apply H(i) to A(i:m,i+1:n) from the left
358: // *
359: if (I < N)
360: {
361: this._dlarf.Run("Left", M - I + 1, N - I, A, I+I * LDA + o_a, 1, TAUQ[I + o_tauq]
362: , ref A, I+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
363: }
364: A[I+I * LDA + o_a] = D[I + o_d];
365: // *
366: if (I < N)
367: {
368: // *
369: // * Generate elementary reflector G(i) to annihilate
370: // * A(i,i+2:n)
371: // *
372: this._dlarfg.Run(N - I, ref A[I+(I + 1) * LDA + o_a], ref A, I+Math.Min(I + 2, N) * LDA + o_a, LDA, ref TAUP[I + o_taup]);
373: E[I + o_e] = A[I+(I + 1) * LDA + o_a];
374: A[I+(I + 1) * LDA + o_a] = ONE;
375: // *
376: // * Apply G(i) to A(i+1:m,i+1:n) from the right
377: // *
378: this._dlarf.Run("Right", M - I, N - I, A, I+(I + 1) * LDA + o_a, LDA, TAUP[I + o_taup]
379: , ref A, I + 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
380: A[I+(I + 1) * LDA + o_a] = E[I + o_e];
381: }
382: else
383: {
384: TAUP[I + o_taup] = ZERO;
385: }
386: }
387: }
388: else
389: {
390: // *
391: // * Reduce to lower bidiagonal form
392: // *
393: for (I = 1; I <= M; I++)
394: {
395: // *
396: // * Generate elementary reflector G(i) to annihilate A(i,i+1:n)
397: // *
398: this._dlarfg.Run(N - I + 1, ref A[I+I * LDA + o_a], ref A, I+Math.Min(I + 1, N) * LDA + o_a, LDA, ref TAUP[I + o_taup]);
399: D[I + o_d] = A[I+I * LDA + o_a];
400: A[I+I * LDA + o_a] = ONE;
401: // *
402: // * Apply G(i) to A(i+1:m,i:n) from the right
403: // *
404: if (I < M)
405: {
406: this._dlarf.Run("Right", M - I, N - I + 1, A, I+I * LDA + o_a, LDA, TAUP[I + o_taup]
407: , ref A, I + 1+I * LDA + o_a, LDA, ref WORK, offset_work);
408: }
409: A[I+I * LDA + o_a] = D[I + o_d];
410: // *
411: if (I < M)
412: {
413: // *
414: // * Generate elementary reflector H(i) to annihilate
415: // * A(i+2:m,i)
416: // *
417: this._dlarfg.Run(M - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, M)+I * LDA + o_a, 1, ref TAUQ[I + o_tauq]);
418: E[I + o_e] = A[I + 1+I * LDA + o_a];
419: A[I + 1+I * LDA + o_a] = ONE;
420: // *
421: // * Apply H(i) to A(i+1:m,i+1:n) from the left
422: // *
423: this._dlarf.Run("Left", M - I, N - I, A, I + 1+I * LDA + o_a, 1, TAUQ[I + o_tauq]
424: , ref A, I + 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
425: A[I + 1+I * LDA + o_a] = E[I + o_e];
426: }
427: else
428: {
429: TAUQ[I + o_tauq] = ZERO;
430: }
431: }
432: }
433: return;
434: // *
435: // * End of DGEBD2
436: // *
437:
438: #endregion
439:
440: }
441: }
442: }