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: /// DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
27: /// form T by an orthogonal similarity transformation: Q' * A * Q = T.
28: ///
29: ///</summary>
30: public class DSYTD2
31: {
32:
33:
34: #region Dependencies
35:
36: DAXPY _daxpy; DLARFG _dlarfg; DSYMV _dsymv; DSYR2 _dsyr2; XERBLA _xerbla; LSAME _lsame; DDOT _ddot;
37:
38: #endregion
39:
40:
41: #region Fields
42:
43: const double ONE = 1.0E0; const double ZERO = 0.0E0; const double HALF = 1.0E0 / 2.0E0; bool UPPER = false; int I = 0;
44: double ALPHA = 0;double TAUI = 0;
45:
46: #endregion
47:
48: public DSYTD2(DAXPY daxpy, DLARFG dlarfg, DSYMV dsymv, DSYR2 dsyr2, XERBLA xerbla, LSAME lsame, DDOT ddot)
49: {
50:
51:
52: #region Set Dependencies
53:
54: this._daxpy = daxpy; this._dlarfg = dlarfg; this._dsymv = dsymv; this._dsyr2 = dsyr2; this._xerbla = xerbla;
55: this._lsame = lsame;this._ddot = ddot;
56:
57: #endregion
58:
59: }
60:
61: public DSYTD2()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: DAXPY daxpy = new DAXPY();
68: LSAME lsame = new LSAME();
69: DLAMC3 dlamc3 = new DLAMC3();
70: DLAPY2 dlapy2 = new DLAPY2();
71: DNRM2 dnrm2 = new DNRM2();
72: DSCAL dscal = new DSCAL();
73: XERBLA xerbla = new XERBLA();
74: DDOT ddot = new DDOT();
75: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
76: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
77: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
78: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
79: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
80: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
81: DSYMV dsymv = new DSYMV(lsame, xerbla);
82: DSYR2 dsyr2 = new DSYR2(lsame, xerbla);
83:
84: #endregion
85:
86:
87: #region Set Dependencies
88:
89: this._daxpy = daxpy; this._dlarfg = dlarfg; this._dsymv = dsymv; this._dsyr2 = dsyr2; this._xerbla = xerbla;
90: this._lsame = lsame;this._ddot = ddot;
91:
92: #endregion
93:
94: }
95: /// <summary>
96: /// Purpose
97: /// =======
98: ///
99: /// DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
100: /// form T by an orthogonal similarity transformation: Q' * A * Q = T.
101: ///
102: ///</summary>
103: /// <param name="UPLO">
104: /// (input) CHARACTER*1
105: /// Specifies whether the upper or lower triangular part of the
106: /// symmetric matrix A is stored:
107: /// = 'U': Upper triangular
108: /// = 'L': Lower triangular
109: ///</param>
110: /// <param name="N">
111: /// (input) INTEGER
112: /// The order of the matrix A. N .GE. 0.
113: ///</param>
114: /// <param name="A">
115: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
116: /// On entry, the symmetric matrix A. If UPLO = 'U', the leading
117: /// n-by-n upper triangular part of A contains the upper
118: /// triangular part of the matrix A, and the strictly lower
119: /// triangular part of A is not referenced. If UPLO = 'L', the
120: /// leading n-by-n lower triangular part of A contains the lower
121: /// triangular part of the matrix A, and the strictly upper
122: /// triangular part of A is not referenced.
123: /// On exit, if UPLO = 'U', the diagonal and first superdiagonal
124: /// of A are overwritten by the corresponding elements of the
125: /// tridiagonal matrix T, and the elements above the first
126: /// superdiagonal, with the array TAU, represent the orthogonal
127: /// matrix Q as a product of elementary reflectors; if UPLO
128: /// = 'L', the diagonal and first subdiagonal of A are over-
129: /// written by the corresponding elements of the tridiagonal
130: /// matrix T, and the elements below the first subdiagonal, with
131: /// the array TAU, represent the orthogonal matrix Q as a product
132: /// of elementary reflectors. See Further Details.
133: ///</param>
134: /// <param name="LDA">
135: /// (input) INTEGER
136: /// The leading dimension of the array A. LDA .GE. max(1,N).
137: ///</param>
138: /// <param name="D">
139: /// (output) DOUBLE PRECISION array, dimension (N)
140: /// The diagonal elements of the tridiagonal matrix T:
141: /// D(i) = A(i,i).
142: ///</param>
143: /// <param name="E">
144: /// (output) DOUBLE PRECISION array, dimension (N-1)
145: /// The off-diagonal elements of the tridiagonal matrix T:
146: /// E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
147: ///</param>
148: /// <param name="TAU">
149: /// (output) DOUBLE PRECISION array, dimension (N-1)
150: /// The scalar factors of the elementary reflectors (see Further
151: /// Details).
152: ///</param>
153: /// <param name="INFO">
154: /// (output) INTEGER
155: /// = 0: successful exit
156: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
157: ///</param>
158: public void Run(string UPLO, int N, ref double[] A, int offset_a, int LDA, ref double[] D, int offset_d, ref double[] E, int offset_e
159: , ref double[] TAU, int offset_tau, ref int INFO)
160: {
161:
162: #region Array Index Correction
163:
164: int o_a = -1 - LDA + offset_a; int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_tau = -1 + offset_tau;
165:
166: #endregion
167:
168:
169: #region Strings
170:
171: UPLO = UPLO.Substring(0, 1);
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: // * DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
192: // * form T by an orthogonal similarity transformation: Q' * A * Q = T.
193: // *
194: // * Arguments
195: // * =========
196: // *
197: // * UPLO (input) CHARACTER*1
198: // * Specifies whether the upper or lower triangular part of the
199: // * symmetric matrix A is stored:
200: // * = 'U': Upper triangular
201: // * = 'L': Lower triangular
202: // *
203: // * N (input) INTEGER
204: // * The order of the matrix A. N >= 0.
205: // *
206: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
207: // * On entry, the symmetric matrix A. If UPLO = 'U', the leading
208: // * n-by-n upper triangular part of A contains the upper
209: // * triangular part of the matrix A, and the strictly lower
210: // * triangular part of A is not referenced. If UPLO = 'L', the
211: // * leading n-by-n lower triangular part of A contains the lower
212: // * triangular part of the matrix A, and the strictly upper
213: // * triangular part of A is not referenced.
214: // * On exit, if UPLO = 'U', the diagonal and first superdiagonal
215: // * of A are overwritten by the corresponding elements of the
216: // * tridiagonal matrix T, and the elements above the first
217: // * superdiagonal, with the array TAU, represent the orthogonal
218: // * matrix Q as a product of elementary reflectors; if UPLO
219: // * = 'L', the diagonal and first subdiagonal of A are over-
220: // * written by the corresponding elements of the tridiagonal
221: // * matrix T, and the elements below the first subdiagonal, with
222: // * the array TAU, represent the orthogonal matrix Q as a product
223: // * of elementary reflectors. See Further Details.
224: // *
225: // * LDA (input) INTEGER
226: // * The leading dimension of the array A. LDA >= max(1,N).
227: // *
228: // * D (output) DOUBLE PRECISION array, dimension (N)
229: // * The diagonal elements of the tridiagonal matrix T:
230: // * D(i) = A(i,i).
231: // *
232: // * E (output) DOUBLE PRECISION array, dimension (N-1)
233: // * The off-diagonal elements of the tridiagonal matrix T:
234: // * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
235: // *
236: // * TAU (output) DOUBLE PRECISION array, dimension (N-1)
237: // * The scalar factors of the elementary reflectors (see Further
238: // * Details).
239: // *
240: // * INFO (output) INTEGER
241: // * = 0: successful exit
242: // * < 0: if INFO = -i, the i-th argument had an illegal value.
243: // *
244: // * Further Details
245: // * ===============
246: // *
247: // * If UPLO = 'U', the matrix Q is represented as a product of elementary
248: // * reflectors
249: // *
250: // * Q = H(n-1) . . . H(2) H(1).
251: // *
252: // * Each H(i) has the form
253: // *
254: // * H(i) = I - tau * v * v'
255: // *
256: // * where tau is a real scalar, and v is a real vector with
257: // * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
258: // * A(1:i-1,i+1), and tau in TAU(i).
259: // *
260: // * If UPLO = 'L', the matrix Q is represented as a product of elementary
261: // * reflectors
262: // *
263: // * Q = H(1) H(2) . . . H(n-1).
264: // *
265: // * Each H(i) has the form
266: // *
267: // * H(i) = I - tau * v * v'
268: // *
269: // * where tau is a real scalar, and v is a real vector with
270: // * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
271: // * and tau in TAU(i).
272: // *
273: // * The contents of A on exit are illustrated by the following examples
274: // * with n = 5:
275: // *
276: // * if UPLO = 'U': if UPLO = 'L':
277: // *
278: // * ( d e v2 v3 v4 ) ( d )
279: // * ( d e v3 v4 ) ( e d )
280: // * ( d e v4 ) ( v1 e d )
281: // * ( d e ) ( v1 v2 e d )
282: // * ( d ) ( v1 v2 v3 e d )
283: // *
284: // * where d and e denote diagonal and off-diagonal elements of T, and vi
285: // * denotes an element of the vector defining H(i).
286: // *
287: // * =====================================================================
288: // *
289: // * .. Parameters ..
290: // * ..
291: // * .. Local Scalars ..
292: // * ..
293: // * .. External Subroutines ..
294: // * ..
295: // * .. External Functions ..
296: // * ..
297: // * .. Intrinsic Functions ..
298: // INTRINSIC MAX, MIN;
299: // * ..
300: // * .. Executable Statements ..
301: // *
302: // * Test the input parameters
303: // *
304:
305: #endregion
306:
307:
308: #region Body
309:
310: INFO = 0;
311: UPPER = this._lsame.Run(UPLO, "U");
312: if (!UPPER && !this._lsame.Run(UPLO, "L"))
313: {
314: INFO = - 1;
315: }
316: else
317: {
318: if (N < 0)
319: {
320: INFO = - 2;
321: }
322: else
323: {
324: if (LDA < Math.Max(1, N))
325: {
326: INFO = - 4;
327: }
328: }
329: }
330: if (INFO != 0)
331: {
332: this._xerbla.Run("DSYTD2", - INFO);
333: return;
334: }
335: // *
336: // * Quick return if possible
337: // *
338: if (N <= 0) return;
339: // *
340: if (UPPER)
341: {
342: // *
343: // * Reduce the upper triangle of A
344: // *
345: for (I = N - 1; I >= 1; I += - 1)
346: {
347: // *
348: // * Generate elementary reflector H(i) = I - tau * v * v'
349: // * to annihilate A(1:i-1,i+1)
350: // *
351: this._dlarfg.Run(I, ref A[I+(I + 1) * LDA + o_a], ref A, 1+(I + 1) * LDA + o_a, 1, ref TAUI);
352: E[I + o_e] = A[I+(I + 1) * LDA + o_a];
353: // *
354: if (TAUI != ZERO)
355: {
356: // *
357: // * Apply H(i) from both sides to A(1:i,1:i)
358: // *
359: A[I+(I + 1) * LDA + o_a] = ONE;
360: // *
361: // * Compute x := tau * A * v storing x in TAU(1:i)
362: // *
363: this._dsymv.Run(UPLO, I, TAUI, A, offset_a, LDA, A, 1+(I + 1) * LDA + o_a
364: , 1, ZERO, ref TAU, offset_tau, 1);
365: // *
366: // * Compute w := x - 1/2 * tau * (x'*v) * v
367: // *
368: ALPHA = - HALF * TAUI * this._ddot.Run(I, TAU, offset_tau, 1, A, 1+(I + 1) * LDA + o_a, 1);
369: this._daxpy.Run(I, ALPHA, A, 1+(I + 1) * LDA + o_a, 1, ref TAU, offset_tau, 1);
370: // *
371: // * Apply the transformation as a rank-2 update:
372: // * A := A - v * w' - w * v'
373: // *
374: this._dsyr2.Run(UPLO, I, - ONE, A, 1+(I + 1) * LDA + o_a, 1, TAU, offset_tau
375: , 1, ref A, offset_a, LDA);
376: // *
377: A[I+(I + 1) * LDA + o_a] = E[I + o_e];
378: }
379: D[I + 1 + o_d] = A[I + 1+(I + 1) * LDA + o_a];
380: TAU[I + o_tau] = TAUI;
381: }
382: D[1 + o_d] = A[1+1 * LDA + o_a];
383: }
384: else
385: {
386: // *
387: // * Reduce the lower triangle of A
388: // *
389: for (I = 1; I <= N - 1; I++)
390: {
391: // *
392: // * Generate elementary reflector H(i) = I - tau * v * v'
393: // * to annihilate A(i+2:n,i)
394: // *
395: this._dlarfg.Run(N - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, N)+I * LDA + o_a, 1, ref TAUI);
396: E[I + o_e] = A[I + 1+I * LDA + o_a];
397: // *
398: if (TAUI != ZERO)
399: {
400: // *
401: // * Apply H(i) from both sides to A(i+1:n,i+1:n)
402: // *
403: A[I + 1+I * LDA + o_a] = ONE;
404: // *
405: // * Compute x := tau * A * v storing y in TAU(i:n-1)
406: // *
407: this._dsymv.Run(UPLO, N - I, TAUI, A, I + 1+(I + 1) * LDA + o_a, LDA, A, I + 1+I * LDA + o_a
408: , 1, ZERO, ref TAU, I + o_tau, 1);
409: // *
410: // * Compute w := x - 1/2 * tau * (x'*v) * v
411: // *
412: ALPHA = - HALF * TAUI * this._ddot.Run(N - I, TAU, I + o_tau, 1, A, I + 1+I * LDA + o_a, 1);
413: this._daxpy.Run(N - I, ALPHA, A, I + 1+I * LDA + o_a, 1, ref TAU, I + o_tau, 1);
414: // *
415: // * Apply the transformation as a rank-2 update:
416: // * A := A - v * w' - w * v'
417: // *
418: this._dsyr2.Run(UPLO, N - I, - ONE, A, I + 1+I * LDA + o_a, 1, TAU, I + o_tau
419: , 1, ref A, I + 1+(I + 1) * LDA + o_a, LDA);
420: // *
421: A[I + 1+I * LDA + o_a] = E[I + o_e];
422: }
423: D[I + o_d] = A[I+I * LDA + o_a];
424: TAU[I + o_tau] = TAUI;
425: }
426: D[N + o_d] = A[N+N * LDA + o_a];
427: }
428: // *
429: return;
430: // *
431: // * End of DSYTD2
432: // *
433:
434: #endregion
435:
436: }
437: }
438: }