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