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