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: /// DSYTRD reduces a real symmetric matrix A to real symmetric
27: /// tridiagonal form T by an orthogonal similarity transformation:
28: /// Q**T * A * Q = T.
29: ///
30: ///</summary>
31: public class DSYTRD
32: {
33:
34:
35: #region Dependencies
36:
37: DLATRD _dlatrd; DSYR2K _dsyr2k; DSYTD2 _dsytd2; XERBLA _xerbla; LSAME _lsame; ILAENV _ilaenv;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double ONE = 1.0E+0; bool LQUERY = false; bool UPPER = false; int I = 0; int IINFO = 0; int IWS = 0; int J = 0;
45: int KK = 0;int LDWORK = 0; int LWKOPT = 0; int NB = 0; int NBMIN = 0; int NX = 0;
46:
47: #endregion
48:
49: public DSYTRD(DLATRD dlatrd, DSYR2K dsyr2k, DSYTD2 dsytd2, XERBLA xerbla, LSAME lsame, ILAENV ilaenv)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._dlatrd = dlatrd; this._dsyr2k = dsyr2k; this._dsytd2 = dsytd2; this._xerbla = xerbla; this._lsame = lsame;
56: this._ilaenv = ilaenv;
57:
58: #endregion
59:
60: }
61:
62: public DSYTRD()
63: {
64:
65:
66: #region Dependencies (Initialization)
67:
68: DAXPY daxpy = new DAXPY();
69: LSAME lsame = new LSAME();
70: XERBLA xerbla = new XERBLA();
71: DLAMC3 dlamc3 = new DLAMC3();
72: DLAPY2 dlapy2 = new DLAPY2();
73: DNRM2 dnrm2 = new DNRM2();
74: DSCAL dscal = new DSCAL();
75: DDOT ddot = new DDOT();
76: IEEECK ieeeck = new IEEECK();
77: IPARMQ iparmq = new IPARMQ();
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: DSYMV dsymv = new DSYMV(lsame, xerbla);
86: DLATRD dlatrd = new DLATRD(daxpy, dgemv, dlarfg, dscal, dsymv, lsame, ddot);
87: DSYR2K dsyr2k = new DSYR2K(lsame, xerbla);
88: DSYR2 dsyr2 = new DSYR2(lsame, xerbla);
89: DSYTD2 dsytd2 = new DSYTD2(daxpy, dlarfg, dsymv, dsyr2, xerbla, lsame, ddot);
90: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
91:
92: #endregion
93:
94:
95: #region Set Dependencies
96:
97: this._dlatrd = dlatrd; this._dsyr2k = dsyr2k; this._dsytd2 = dsytd2; this._xerbla = xerbla; this._lsame = lsame;
98: this._ilaenv = ilaenv;
99:
100: #endregion
101:
102: }
103: /// <summary>
104: /// Purpose
105: /// =======
106: ///
107: /// DSYTRD reduces a real symmetric matrix A to real symmetric
108: /// tridiagonal form T by an orthogonal similarity transformation:
109: /// Q**T * A * Q = T.
110: ///
111: ///</summary>
112: /// <param name="UPLO">
113: /// (input) CHARACTER*1
114: /// = 'U': Upper triangle of A is stored;
115: /// = 'L': Lower triangle of A is stored.
116: ///</param>
117: /// <param name="N">
118: /// (input) INTEGER
119: /// The order of the matrix A. N .GE. 0.
120: ///</param>
121: /// <param name="A">
122: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
123: /// On entry, the symmetric matrix A. If UPLO = 'U', the leading
124: /// N-by-N upper triangular part of A contains the upper
125: /// triangular part of the matrix A, and the strictly lower
126: /// triangular part of A is not referenced. If UPLO = 'L', the
127: /// leading N-by-N lower triangular part of A contains the lower
128: /// triangular part of the matrix A, and the strictly upper
129: /// triangular part of A is not referenced.
130: /// On exit, if UPLO = 'U', the diagonal and first superdiagonal
131: /// of A are overwritten by the corresponding elements of the
132: /// tridiagonal matrix T, and the elements above the first
133: /// superdiagonal, with the array TAU, represent the orthogonal
134: /// matrix Q as a product of elementary reflectors; if UPLO
135: /// = 'L', the diagonal and first subdiagonal of A are over-
136: /// written by the corresponding elements of the tridiagonal
137: /// matrix T, and the elements below the first subdiagonal, with
138: /// the array TAU, represent the orthogonal matrix Q as a product
139: /// of elementary reflectors. See Further Details.
140: ///</param>
141: /// <param name="LDA">
142: /// (input) INTEGER
143: /// The leading dimension of the array A. LDA .GE. max(1,N).
144: ///</param>
145: /// <param name="D">
146: /// (output) DOUBLE PRECISION array, dimension (N)
147: /// The diagonal elements of the tridiagonal matrix T:
148: /// D(i) = A(i,i).
149: ///</param>
150: /// <param name="E">
151: /// (output) DOUBLE PRECISION array, dimension (N-1)
152: /// The off-diagonal elements of the tridiagonal matrix T:
153: /// E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
154: ///</param>
155: /// <param name="TAU">
156: /// (output) DOUBLE PRECISION array, dimension (N-1)
157: /// The scalar factors of the elementary reflectors (see Further
158: /// Details).
159: ///</param>
160: /// <param name="WORK">
161: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
162: /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
163: ///</param>
164: /// <param name="LWORK">
165: /// (input) INTEGER
166: /// The dimension of the array WORK. LWORK .GE. 1.
167: /// For optimum performance LWORK .GE. N*NB, where NB is the
168: /// optimal blocksize.
169: ///
170: /// If LWORK = -1, then a workspace query is assumed; the routine
171: /// only calculates the optimal size of the WORK array, returns
172: /// this value as the first entry of the WORK array, and no error
173: /// message related to LWORK is issued by XERBLA.
174: ///</param>
175: /// <param name="INFO">
176: /// (output) INTEGER
177: /// = 0: successful exit
178: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
179: ///</param>
180: 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
181: , ref double[] TAU, int offset_tau, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
182: {
183:
184: #region Array Index Correction
185:
186: 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;
187: int o_work = -1 + offset_work;
188:
189: #endregion
190:
191:
192: #region Strings
193:
194: UPLO = UPLO.Substring(0, 1);
195:
196: #endregion
197:
198:
199: #region Prolog
200:
201: // *
202: // * -- LAPACK routine (version 3.1) --
203: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
204: // * November 2006
205: // *
206: // * .. Scalar Arguments ..
207: // * ..
208: // * .. Array Arguments ..
209: // * ..
210: // *
211: // * Purpose
212: // * =======
213: // *
214: // * DSYTRD reduces a real symmetric matrix A to real symmetric
215: // * tridiagonal form T by an orthogonal similarity transformation:
216: // * Q**T * A * Q = T.
217: // *
218: // * Arguments
219: // * =========
220: // *
221: // * UPLO (input) CHARACTER*1
222: // * = 'U': Upper triangle of A is stored;
223: // * = 'L': Lower triangle of A is stored.
224: // *
225: // * N (input) INTEGER
226: // * The order of the matrix A. N >= 0.
227: // *
228: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
229: // * On entry, the symmetric matrix A. If UPLO = 'U', the leading
230: // * N-by-N upper triangular part of A contains the upper
231: // * triangular part of the matrix A, and the strictly lower
232: // * triangular part of A is not referenced. If UPLO = 'L', the
233: // * leading N-by-N lower triangular part of A contains the lower
234: // * triangular part of the matrix A, and the strictly upper
235: // * triangular part of A is not referenced.
236: // * On exit, if UPLO = 'U', the diagonal and first superdiagonal
237: // * of A are overwritten by the corresponding elements of the
238: // * tridiagonal matrix T, and the elements above the first
239: // * superdiagonal, with the array TAU, represent the orthogonal
240: // * matrix Q as a product of elementary reflectors; if UPLO
241: // * = 'L', the diagonal and first subdiagonal of A are over-
242: // * written by the corresponding elements of the tridiagonal
243: // * matrix T, and the elements below the first subdiagonal, with
244: // * the array TAU, represent the orthogonal matrix Q as a product
245: // * of elementary reflectors. See Further Details.
246: // *
247: // * LDA (input) INTEGER
248: // * The leading dimension of the array A. LDA >= max(1,N).
249: // *
250: // * D (output) DOUBLE PRECISION array, dimension (N)
251: // * The diagonal elements of the tridiagonal matrix T:
252: // * D(i) = A(i,i).
253: // *
254: // * E (output) DOUBLE PRECISION array, dimension (N-1)
255: // * The off-diagonal elements of the tridiagonal matrix T:
256: // * E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
257: // *
258: // * TAU (output) DOUBLE PRECISION array, dimension (N-1)
259: // * The scalar factors of the elementary reflectors (see Further
260: // * Details).
261: // *
262: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
263: // * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
264: // *
265: // * LWORK (input) INTEGER
266: // * The dimension of the array WORK. LWORK >= 1.
267: // * For optimum performance LWORK >= N*NB, where NB is the
268: // * optimal blocksize.
269: // *
270: // * If LWORK = -1, then a workspace query is assumed; the routine
271: // * only calculates the optimal size of the WORK array, returns
272: // * this value as the first entry of the WORK array, and no error
273: // * message related to LWORK is issued by XERBLA.
274: // *
275: // * INFO (output) INTEGER
276: // * = 0: successful exit
277: // * < 0: if INFO = -i, the i-th argument had an illegal value
278: // *
279: // * Further Details
280: // * ===============
281: // *
282: // * If UPLO = 'U', the matrix Q is represented as a product of elementary
283: // * reflectors
284: // *
285: // * Q = H(n-1) . . . H(2) H(1).
286: // *
287: // * Each H(i) has the form
288: // *
289: // * H(i) = I - tau * v * v'
290: // *
291: // * where tau is a real scalar, and v is a real vector with
292: // * v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
293: // * A(1:i-1,i+1), and tau in TAU(i).
294: // *
295: // * If UPLO = 'L', the matrix Q is represented as a product of elementary
296: // * reflectors
297: // *
298: // * Q = H(1) H(2) . . . H(n-1).
299: // *
300: // * Each H(i) has the form
301: // *
302: // * H(i) = I - tau * v * v'
303: // *
304: // * where tau is a real scalar, and v is a real vector with
305: // * v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
306: // * and tau in TAU(i).
307: // *
308: // * The contents of A on exit are illustrated by the following examples
309: // * with n = 5:
310: // *
311: // * if UPLO = 'U': if UPLO = 'L':
312: // *
313: // * ( d e v2 v3 v4 ) ( d )
314: // * ( d e v3 v4 ) ( e d )
315: // * ( d e v4 ) ( v1 e d )
316: // * ( d e ) ( v1 v2 e d )
317: // * ( d ) ( v1 v2 v3 e d )
318: // *
319: // * where d and e denote diagonal and off-diagonal elements of T, and vi
320: // * denotes an element of the vector defining H(i).
321: // *
322: // * =====================================================================
323: // *
324: // * .. Parameters ..
325: // * ..
326: // * .. Local Scalars ..
327: // * ..
328: // * .. External Subroutines ..
329: // * ..
330: // * .. Intrinsic Functions ..
331: // INTRINSIC MAX;
332: // * ..
333: // * .. External Functions ..
334: // * ..
335: // * .. Executable Statements ..
336: // *
337: // * Test the input parameters
338: // *
339:
340: #endregion
341:
342:
343: #region Body
344:
345: INFO = 0;
346: UPPER = this._lsame.Run(UPLO, "U");
347: LQUERY = (LWORK == - 1);
348: if (!UPPER && !this._lsame.Run(UPLO, "L"))
349: {
350: INFO = - 1;
351: }
352: else
353: {
354: if (N < 0)
355: {
356: INFO = - 2;
357: }
358: else
359: {
360: if (LDA < Math.Max(1, N))
361: {
362: INFO = - 4;
363: }
364: else
365: {
366: if (LWORK < 1 && !LQUERY)
367: {
368: INFO = - 9;
369: }
370: }
371: }
372: }
373: // *
374: if (INFO == 0)
375: {
376: // *
377: // * Determine the block size.
378: // *
379: NB = this._ilaenv.Run(1, "DSYTRD", UPLO, N, - 1, - 1, - 1);
380: LWKOPT = N * NB;
381: WORK[1 + o_work] = LWKOPT;
382: }
383: // *
384: if (INFO != 0)
385: {
386: this._xerbla.Run("DSYTRD", - INFO);
387: return;
388: }
389: else
390: {
391: if (LQUERY)
392: {
393: return;
394: }
395: }
396: // *
397: // * Quick return if possible
398: // *
399: if (N == 0)
400: {
401: WORK[1 + o_work] = 1;
402: return;
403: }
404: // *
405: NX = N;
406: IWS = 1;
407: if (NB > 1 && NB < N)
408: {
409: // *
410: // * Determine when to cross over from blocked to unblocked code
411: // * (last block is always handled by unblocked code).
412: // *
413: NX = Math.Max(NB, this._ilaenv.Run(3, "DSYTRD", UPLO, N, - 1, - 1, - 1));
414: if (NX < N)
415: {
416: // *
417: // * Determine if workspace is large enough for blocked code.
418: // *
419: LDWORK = N;
420: IWS = LDWORK * NB;
421: if (LWORK < IWS)
422: {
423: // *
424: // * Not enough workspace to use optimal NB: determine the
425: // * minimum value of NB, and reduce NB or force use of
426: // * unblocked code by setting NX = N.
427: // *
428: NB = Math.Max(LWORK / LDWORK, 1);
429: NBMIN = this._ilaenv.Run(2, "DSYTRD", UPLO, N, - 1, - 1, - 1);
430: if (NB < NBMIN) NX = N;
431: }
432: }
433: else
434: {
435: NX = N;
436: }
437: }
438: else
439: {
440: NB = 1;
441: }
442: // *
443: if (UPPER)
444: {
445: // *
446: // * Reduce the upper triangle of A.
447: // * Columns 1:kk are handled by the unblocked method.
448: // *
449: KK = N - ((N - NX + NB - 1) / NB) * NB;
450: for (I = N - NB + 1; ( - NB >= 0) ? (I <= KK + 1) : (I >= KK + 1); I += - NB)
451: {
452: // *
453: // * Reduce columns i:i+nb-1 to tridiagonal form and form the
454: // * matrix W which is needed to update the unreduced part of
455: // * the matrix
456: // *
457: this._dlatrd.Run(UPLO, I + NB - 1, NB, ref A, offset_a, LDA, ref E, offset_e
458: , ref TAU, offset_tau, ref WORK, offset_work, LDWORK);
459: // *
460: // * Update the unreduced submatrix A(1:i-1,1:i-1), using an
461: // * update of the form: A := A - V*W' - W*V'
462: // *
463: this._dsyr2k.Run(UPLO, "No transpose", I - 1, NB, - ONE, A, 1+I * LDA + o_a
464: , LDA, WORK, offset_work, LDWORK, ONE, ref A, offset_a, LDA);
465: // *
466: // * Copy superdiagonal elements back into A, and diagonal
467: // * elements into D
468: // *
469: for (J = I; J <= I + NB - 1; J++)
470: {
471: A[J - 1+J * LDA + o_a] = E[J - 1 + o_e];
472: D[J + o_d] = A[J+J * LDA + o_a];
473: }
474: }
475: // *
476: // * Use unblocked code to reduce the last or only block
477: // *
478: this._dsytd2.Run(UPLO, KK, ref A, offset_a, LDA, ref D, offset_d, ref E, offset_e
479: , ref TAU, offset_tau, ref IINFO);
480: }
481: else
482: {
483: // *
484: // * Reduce the lower triangle of A
485: // *
486: for (I = 1; (NB >= 0) ? (I <= N - NX) : (I >= N - NX); I += NB)
487: {
488: // *
489: // * Reduce columns i:i+nb-1 to tridiagonal form and form the
490: // * matrix W which is needed to update the unreduced part of
491: // * the matrix
492: // *
493: this._dlatrd.Run(UPLO, N - I + 1, NB, ref A, I+I * LDA + o_a, LDA, ref E, I + o_e
494: , ref TAU, I + o_tau, ref WORK, offset_work, LDWORK);
495: // *
496: // * Update the unreduced submatrix A(i+ib:n,i+ib:n), using
497: // * an update of the form: A := A - V*W' - W*V'
498: // *
499: this._dsyr2k.Run(UPLO, "No transpose", N - I - NB + 1, NB, - ONE, A, I + NB+I * LDA + o_a
500: , LDA, WORK, NB + 1 + o_work, LDWORK, ONE, ref A, I + NB+(I + NB) * LDA + o_a, LDA);
501: // *
502: // * Copy subdiagonal elements back into A, and diagonal
503: // * elements into D
504: // *
505: for (J = I; J <= I + NB - 1; J++)
506: {
507: A[J + 1+J * LDA + o_a] = E[J + o_e];
508: D[J + o_d] = A[J+J * LDA + o_a];
509: }
510: }
511: // *
512: // * Use unblocked code to reduce the last or only block
513: // *
514: this._dsytd2.Run(UPLO, N - I + 1, ref A, I+I * LDA + o_a, LDA, ref D, I + o_d, ref E, I + o_e
515: , ref TAU, I + o_tau, ref IINFO);
516: }
517: // *
518: WORK[1 + o_work] = LWKOPT;
519: return;
520: // *
521: // * End of DSYTRD
522: // *
523:
524: #endregion
525:
526: }
527: }
528: }