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: /// DSBTRD reduces a real symmetric band matrix A to symmetric
27: /// tridiagonal form T by an orthogonal similarity transformation:
28: /// Q**T * A * Q = T.
29: ///
30: ///</summary>
31: public class DSBTRD
32: {
33:
34:
35: #region Dependencies
36:
37: DLAR2V _dlar2v; DLARGV _dlargv; DLARTG _dlartg; DLARTV _dlartv; DLASET _dlaset; DROT _drot; XERBLA _xerbla; LSAME _lsame;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool INITQ = false; bool UPPER = false; bool WANTQ = false;
45: int I = 0;int I2 = 0; int IBL = 0; int INCA = 0; int INCX = 0; int IQAEND = 0; int IQB = 0; int IQEND = 0; int J = 0;
46: int J1 = 0;int J1END = 0; int J1INC = 0; int J2 = 0; int JEND = 0; int JIN = 0; int JINC = 0; int K = 0; int KD1 = 0;
47: int KDM1 = 0;int KDN = 0; int L = 0; int LAST = 0; int LEND = 0; int NQ = 0; int NR = 0; int NRT = 0; double TEMP = 0;
48:
49: #endregion
50:
51: public DSBTRD(DLAR2V dlar2v, DLARGV dlargv, DLARTG dlartg, DLARTV dlartv, DLASET dlaset, DROT drot, XERBLA xerbla, LSAME lsame)
52: {
53:
54:
55: #region Set Dependencies
56:
57: this._dlar2v = dlar2v; this._dlargv = dlargv; this._dlartg = dlartg; this._dlartv = dlartv; this._dlaset = dlaset;
58: this._drot = drot;this._xerbla = xerbla; this._lsame = lsame;
59:
60: #endregion
61:
62: }
63:
64: public DSBTRD()
65: {
66:
67:
68: #region Dependencies (Initialization)
69:
70: DLAR2V dlar2v = new DLAR2V();
71: DLARGV dlargv = new DLARGV();
72: LSAME lsame = new LSAME();
73: DLAMC3 dlamc3 = new DLAMC3();
74: DLARTV dlartv = new DLARTV();
75: DROT drot = new DROT();
76: XERBLA xerbla = new XERBLA();
77: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
78: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
79: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
80: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
81: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
82: DLARTG dlartg = new DLARTG(dlamch);
83: DLASET dlaset = new DLASET(lsame);
84:
85: #endregion
86:
87:
88: #region Set Dependencies
89:
90: this._dlar2v = dlar2v; this._dlargv = dlargv; this._dlartg = dlartg; this._dlartv = dlartv; this._dlaset = dlaset;
91: this._drot = drot;this._xerbla = xerbla; this._lsame = lsame;
92:
93: #endregion
94:
95: }
96: /// <summary>
97: /// Purpose
98: /// =======
99: ///
100: /// DSBTRD reduces a real symmetric band matrix A to symmetric
101: /// tridiagonal form T by an orthogonal similarity transformation:
102: /// Q**T * A * Q = T.
103: ///
104: ///</summary>
105: /// <param name="VECT">
106: /// (input) CHARACTER*1
107: /// = 'N': do not form Q;
108: /// = 'V': form Q;
109: /// = 'U': update a matrix X, by forming X*Q.
110: ///</param>
111: /// <param name="UPLO">
112: /// (input) CHARACTER*1
113: /// = 'U': Upper triangle of A is stored;
114: /// = 'L': Lower triangle of A is stored.
115: ///</param>
116: /// <param name="N">
117: /// (input) INTEGER
118: /// The order of the matrix A. N .GE. 0.
119: ///</param>
120: /// <param name="KD">
121: /// (input) INTEGER
122: /// The number of superdiagonals of the matrix A if UPLO = 'U',
123: /// or the number of subdiagonals if UPLO = 'L'. KD .GE. 0.
124: ///</param>
125: /// <param name="AB">
126: /// (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
127: /// On entry, the upper or lower triangle of the symmetric band
128: /// matrix A, stored in the first KD+1 rows of the array. The
129: /// j-th column of A is stored in the j-th column of the array AB
130: /// as follows:
131: /// if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd).LE.i.LE.j;
132: /// if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j.LE.i.LE.min(n,j+kd).
133: /// On exit, the diagonal elements of AB are overwritten by the
134: /// diagonal elements of the tridiagonal matrix T; if KD .GT. 0, the
135: /// elements on the first superdiagonal (if UPLO = 'U') or the
136: /// first subdiagonal (if UPLO = 'L') are overwritten by the
137: /// off-diagonal elements of T; the rest of AB is overwritten by
138: /// values generated during the reduction.
139: ///</param>
140: /// <param name="LDAB">
141: /// (input) INTEGER
142: /// The leading dimension of the array AB. LDAB .GE. KD+1.
143: ///</param>
144: /// <param name="D">
145: /// (output) DOUBLE PRECISION array, dimension (N)
146: /// The diagonal elements of the tridiagonal matrix T.
147: ///</param>
148: /// <param name="E">
149: /// (output) DOUBLE PRECISION array, dimension (N-1)
150: /// The off-diagonal elements of the tridiagonal matrix T:
151: /// E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
152: ///</param>
153: /// <param name="Q">
154: /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
155: /// On entry, if VECT = 'U', then Q must contain an N-by-N
156: /// matrix X; if VECT = 'N' or 'V', then Q need not be set.
157: ///
158: /// On exit:
159: /// if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
160: /// if VECT = 'U', Q contains the product X*Q;
161: /// if VECT = 'N', the array Q is not referenced.
162: ///</param>
163: /// <param name="LDQ">
164: /// (input) INTEGER
165: /// The leading dimension of the array Q.
166: /// LDQ .GE. 1, and LDQ .GE. N if VECT = 'V' or 'U'.
167: ///</param>
168: /// <param name="WORK">
169: /// (workspace) DOUBLE PRECISION array, dimension (N)
170: ///</param>
171: /// <param name="INFO">
172: /// (output) INTEGER
173: /// = 0: successful exit
174: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
175: ///</param>
176: public void Run(string VECT, string UPLO, int N, int KD, ref double[] AB, int offset_ab, int LDAB
177: , ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] Q, int offset_q, int LDQ, ref double[] WORK, int offset_work, ref int INFO)
178: {
179:
180: #region Array Index Correction
181:
182: int o_ab = -1 - LDAB + offset_ab; int o_d = -1 + offset_d; int o_e = -1 + offset_e; int o_q = -1 - LDQ + offset_q;
183: int o_work = -1 + offset_work;
184:
185: #endregion
186:
187:
188: #region Strings
189:
190: VECT = VECT.Substring(0, 1); UPLO = UPLO.Substring(0, 1);
191:
192: #endregion
193:
194:
195: #region Prolog
196:
197: // *
198: // * -- LAPACK routine (version 3.1) --
199: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
200: // * November 2006
201: // *
202: // * .. Scalar Arguments ..
203: // * ..
204: // * .. Array Arguments ..
205: // * ..
206: // *
207: // * Purpose
208: // * =======
209: // *
210: // * DSBTRD reduces a real symmetric band matrix A to symmetric
211: // * tridiagonal form T by an orthogonal similarity transformation:
212: // * Q**T * A * Q = T.
213: // *
214: // * Arguments
215: // * =========
216: // *
217: // * VECT (input) CHARACTER*1
218: // * = 'N': do not form Q;
219: // * = 'V': form Q;
220: // * = 'U': update a matrix X, by forming X*Q.
221: // *
222: // * UPLO (input) CHARACTER*1
223: // * = 'U': Upper triangle of A is stored;
224: // * = 'L': Lower triangle of A is stored.
225: // *
226: // * N (input) INTEGER
227: // * The order of the matrix A. N >= 0.
228: // *
229: // * KD (input) INTEGER
230: // * The number of superdiagonals of the matrix A if UPLO = 'U',
231: // * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
232: // *
233: // * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
234: // * On entry, the upper or lower triangle of the symmetric band
235: // * matrix A, stored in the first KD+1 rows of the array. The
236: // * j-th column of A is stored in the j-th column of the array AB
237: // * as follows:
238: // * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
239: // * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
240: // * On exit, the diagonal elements of AB are overwritten by the
241: // * diagonal elements of the tridiagonal matrix T; if KD > 0, the
242: // * elements on the first superdiagonal (if UPLO = 'U') or the
243: // * first subdiagonal (if UPLO = 'L') are overwritten by the
244: // * off-diagonal elements of T; the rest of AB is overwritten by
245: // * values generated during the reduction.
246: // *
247: // * LDAB (input) INTEGER
248: // * The leading dimension of the array AB. LDAB >= KD+1.
249: // *
250: // * D (output) DOUBLE PRECISION array, dimension (N)
251: // * The diagonal elements of the tridiagonal matrix T.
252: // *
253: // * E (output) DOUBLE PRECISION array, dimension (N-1)
254: // * The off-diagonal elements of the tridiagonal matrix T:
255: // * E(i) = T(i,i+1) if UPLO = 'U'; E(i) = T(i+1,i) if UPLO = 'L'.
256: // *
257: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
258: // * On entry, if VECT = 'U', then Q must contain an N-by-N
259: // * matrix X; if VECT = 'N' or 'V', then Q need not be set.
260: // *
261: // * On exit:
262: // * if VECT = 'V', Q contains the N-by-N orthogonal matrix Q;
263: // * if VECT = 'U', Q contains the product X*Q;
264: // * if VECT = 'N', the array Q is not referenced.
265: // *
266: // * LDQ (input) INTEGER
267: // * The leading dimension of the array Q.
268: // * LDQ >= 1, and LDQ >= N if VECT = 'V' or 'U'.
269: // *
270: // * WORK (workspace) DOUBLE PRECISION array, dimension (N)
271: // *
272: // * INFO (output) INTEGER
273: // * = 0: successful exit
274: // * < 0: if INFO = -i, the i-th argument had an illegal value
275: // *
276: // * Further Details
277: // * ===============
278: // *
279: // * Modified by Linda Kaufman, Bell Labs.
280: // *
281: // * =====================================================================
282: // *
283: // * .. Parameters ..
284: // * ..
285: // * .. Local Scalars ..
286: // * ..
287: // * .. External Subroutines ..
288: // * ..
289: // * .. Intrinsic Functions ..
290: // INTRINSIC MAX, MIN;
291: // * ..
292: // * .. External Functions ..
293: // * ..
294: // * .. Executable Statements ..
295: // *
296: // * Test the input parameters
297: // *
298:
299: #endregion
300:
301:
302: #region Body
303:
304: INITQ = this._lsame.Run(VECT, "V");
305: WANTQ = INITQ || this._lsame.Run(VECT, "U");
306: UPPER = this._lsame.Run(UPLO, "U");
307: KD1 = KD + 1;
308: KDM1 = KD - 1;
309: INCX = LDAB - 1;
310: IQEND = 1;
311: // *
312: INFO = 0;
313: if (!WANTQ && !this._lsame.Run(VECT, "N"))
314: {
315: INFO = - 1;
316: }
317: else
318: {
319: if (!UPPER && !this._lsame.Run(UPLO, "L"))
320: {
321: INFO = - 2;
322: }
323: else
324: {
325: if (N < 0)
326: {
327: INFO = - 3;
328: }
329: else
330: {
331: if (KD < 0)
332: {
333: INFO = - 4;
334: }
335: else
336: {
337: if (LDAB < KD1)
338: {
339: INFO = - 6;
340: }
341: else
342: {
343: if (LDQ < Math.Max(1, N) && WANTQ)
344: {
345: INFO = - 10;
346: }
347: }
348: }
349: }
350: }
351: }
352: if (INFO != 0)
353: {
354: this._xerbla.Run("DSBTRD", - INFO);
355: return;
356: }
357: // *
358: // * Quick return if possible
359: // *
360: if (N == 0) return;
361: // *
362: // * Initialize Q to the unit matrix, if needed
363: // *
364: if (INITQ)
365: {
366: this._dlaset.Run("Full", N, N, ZERO, ONE, ref Q, offset_q
367: , LDQ);
368: }
369: // *
370: // * Wherever possible, plane rotations are generated and applied in
371: // * vector operations of length NR over the index set J1:J2:KD1.
372: // *
373: // * The cosines and sines of the plane rotations are stored in the
374: // * arrays D and WORK.
375: // *
376: INCA = KD1 * LDAB;
377: KDN = Math.Min(N - 1, KD);
378: if (UPPER)
379: {
380: // *
381: if (KD > 1)
382: {
383: // *
384: // * Reduce to tridiagonal form, working with upper triangle
385: // *
386: NR = 0;
387: J1 = KDN + 2;
388: J2 = 1;
389: // *
390: for (I = 1; I <= N - 2; I++)
391: {
392: // *
393: // * Reduce i-th row of matrix to tridiagonal form
394: // *
395: for (K = KDN + 1; K >= 2; K += - 1)
396: {
397: J1 = J1 + KDN;
398: J2 = J2 + KDN;
399: // *
400: if (NR > 0)
401: {
402: // *
403: // * generate plane rotations to annihilate nonzero
404: // * elements which have been created outside the band
405: // *
406: this._dlargv.Run(NR, ref AB, 1+(J1 - 1) * LDAB + o_ab, INCA, ref WORK, J1 + o_work, KD1, ref D, J1 + o_d
407: , KD1);
408: // *
409: // * apply rotations from the right
410: // *
411: // *
412: // * Dependent on the the number of diagonals either
413: // * DLARTV or DROT is used
414: // *
415: if (NR >= 2 * KD - 1)
416: {
417: for (L = 1; L <= KD - 1; L++)
418: {
419: this._dlartv.Run(NR, ref AB, L + 1+(J1 - 1) * LDAB + o_ab, INCA, ref AB, L+J1 * LDAB + o_ab, INCA, D, J1 + o_d
420: , WORK, J1 + o_work, KD1);
421: }
422: // *
423: }
424: else
425: {
426: JEND = J1 + (NR - 1) * KD1;
427: for (JINC = J1; (KD1 >= 0) ? (JINC <= JEND) : (JINC >= JEND); JINC += KD1)
428: {
429: this._drot.Run(KDM1, ref AB, 2+(JINC - 1) * LDAB + o_ab, 1, ref AB, 1+JINC * LDAB + o_ab, 1, D[JINC + o_d]
430: , WORK[JINC + o_work]);
431: }
432: }
433: }
434: // *
435: // *
436: if (K > 2)
437: {
438: if (K <= N - I + 1)
439: {
440: // *
441: // * generate plane rotation to annihilate a(i,i+k-1)
442: // * within the band
443: // *
444: this._dlartg.Run(AB[KD - K + 3+(I + K - 2) * LDAB + o_ab], AB[KD - K + 2+(I + K - 1) * LDAB + o_ab], ref D[I + K - 1 + o_d], ref WORK[I + K - 1 + o_work], ref TEMP);
445: AB[KD - K + 3+(I + K - 2) * LDAB + o_ab] = TEMP;
446: // *
447: // * apply rotation from the right
448: // *
449: this._drot.Run(K - 3, ref AB, KD - K + 4+(I + K - 2) * LDAB + o_ab, 1, ref AB, KD - K + 3+(I + K - 1) * LDAB + o_ab, 1, D[I + K - 1 + o_d]
450: , WORK[I + K - 1 + o_work]);
451: }
452: NR = NR + 1;
453: J1 = J1 - KDN - 1;
454: }
455: // *
456: // * apply plane rotations from both sides to diagonal
457: // * blocks
458: // *
459: if (NR > 0)
460: {
461: this._dlar2v.Run(NR, ref AB, KD1+(J1 - 1) * LDAB + o_ab, ref AB, KD1+J1 * LDAB + o_ab, ref AB, KD+J1 * LDAB + o_ab, INCA, D, J1 + o_d
462: , WORK, J1 + o_work, KD1);
463: }
464: // *
465: // * apply plane rotations from the left
466: // *
467: if (NR > 0)
468: {
469: if (2 * KD - 1 < NR)
470: {
471: // *
472: // * Dependent on the the number of diagonals either
473: // * DLARTV or DROT is used
474: // *
475: for (L = 1; L <= KD - 1; L++)
476: {
477: if (J2 + L > N)
478: {
479: NRT = NR - 1;
480: }
481: else
482: {
483: NRT = NR;
484: }
485: if (NRT > 0)
486: {
487: this._dlartv.Run(NRT, ref AB, KD - L+(J1 + L) * LDAB + o_ab, INCA, ref AB, KD - L + 1+(J1 + L) * LDAB + o_ab, INCA, D, J1 + o_d
488: , WORK, J1 + o_work, KD1);
489: }
490: }
491: }
492: else
493: {
494: J1END = J1 + KD1 * (NR - 2);
495: if (J1END >= J1)
496: {
497: for (JIN = J1; (KD1 >= 0) ? (JIN <= J1END) : (JIN >= J1END); JIN += KD1)
498: {
499: this._drot.Run(KD - 1, ref AB, KD - 1+(JIN + 1) * LDAB + o_ab, INCX, ref AB, KD+(JIN + 1) * LDAB + o_ab, INCX, D[JIN + o_d]
500: , WORK[JIN + o_work]);
501: }
502: }
503: LEND = Math.Min(KDM1, N - J2);
504: LAST = J1END + KD1;
505: if (LEND > 0)
506: {
507: this._drot.Run(LEND, ref AB, KD - 1+(LAST + 1) * LDAB + o_ab, INCX, ref AB, KD+(LAST + 1) * LDAB + o_ab, INCX, D[LAST + o_d]
508: , WORK[LAST + o_work]);
509: }
510: }
511: }
512: // *
513: if (WANTQ)
514: {
515: // *
516: // * accumulate product of plane rotations in Q
517: // *
518: if (INITQ)
519: {
520: // *
521: // * take advantage of the fact that Q was
522: // * initially the Identity matrix
523: // *
524: IQEND = Math.Max(IQEND, J2);
525: I2 = Math.Max(0, K - 3);
526: IQAEND = 1 + I * KD;
527: if (K == 2) IQAEND = IQAEND + KD;
528: IQAEND = Math.Min(IQAEND, IQEND);
529: for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
530: {
531: IBL = I - I2 / KDM1;
532: I2 = I2 + 1;
533: IQB = Math.Max(1, J - IBL);
534: NQ = 1 + IQAEND - IQB;
535: IQAEND = Math.Min(IQAEND + KD, IQEND);
536: this._drot.Run(NQ, ref Q, IQB+(J - 1) * LDQ + o_q, 1, ref Q, IQB+J * LDQ + o_q, 1, D[J + o_d]
537: , WORK[J + o_work]);
538: }
539: }
540: else
541: {
542: // *
543: for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
544: {
545: this._drot.Run(N, ref Q, 1+(J - 1) * LDQ + o_q, 1, ref Q, 1+J * LDQ + o_q, 1, D[J + o_d]
546: , WORK[J + o_work]);
547: }
548: }
549: // *
550: }
551: // *
552: if (J2 + KDN > N)
553: {
554: // *
555: // * adjust J2 to keep within the bounds of the matrix
556: // *
557: NR = NR - 1;
558: J2 = J2 - KDN - 1;
559: }
560: // *
561: for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
562: {
563: // *
564: // * create nonzero element a(j-1,j+kd) outside the band
565: // * and store it in WORK
566: // *
567: WORK[J + KD + o_work] = WORK[J + o_work] * AB[1+(J + KD) * LDAB + o_ab];
568: AB[1+(J + KD) * LDAB + o_ab] = D[J + o_d] * AB[1+(J + KD) * LDAB + o_ab];
569: }
570: }
571: }
572: }
573: // *
574: if (KD > 0)
575: {
576: // *
577: // * copy off-diagonal elements to E
578: // *
579: for (I = 1; I <= N - 1; I++)
580: {
581: E[I + o_e] = AB[KD+(I + 1) * LDAB + o_ab];
582: }
583: }
584: else
585: {
586: // *
587: // * set E to zero if original matrix was diagonal
588: // *
589: for (I = 1; I <= N - 1; I++)
590: {
591: E[I + o_e] = ZERO;
592: }
593: }
594: // *
595: // * copy diagonal elements to D
596: // *
597: for (I = 1; I <= N; I++)
598: {
599: D[I + o_d] = AB[KD1+I * LDAB + o_ab];
600: }
601: // *
602: }
603: else
604: {
605: // *
606: if (KD > 1)
607: {
608: // *
609: // * Reduce to tridiagonal form, working with lower triangle
610: // *
611: NR = 0;
612: J1 = KDN + 2;
613: J2 = 1;
614: // *
615: for (I = 1; I <= N - 2; I++)
616: {
617: // *
618: // * Reduce i-th column of matrix to tridiagonal form
619: // *
620: for (K = KDN + 1; K >= 2; K += - 1)
621: {
622: J1 = J1 + KDN;
623: J2 = J2 + KDN;
624: // *
625: if (NR > 0)
626: {
627: // *
628: // * generate plane rotations to annihilate nonzero
629: // * elements which have been created outside the band
630: // *
631: this._dlargv.Run(NR, ref AB, KD1+(J1 - KD1) * LDAB + o_ab, INCA, ref WORK, J1 + o_work, KD1, ref D, J1 + o_d
632: , KD1);
633: // *
634: // * apply plane rotations from one side
635: // *
636: // *
637: // * Dependent on the the number of diagonals either
638: // * DLARTV or DROT is used
639: // *
640: if (NR > 2 * KD - 1)
641: {
642: for (L = 1; L <= KD - 1; L++)
643: {
644: this._dlartv.Run(NR, ref AB, KD1 - L+(J1 - KD1 + L) * LDAB + o_ab, INCA, ref AB, KD1 - L + 1+(J1 - KD1 + L) * LDAB + o_ab, INCA, D, J1 + o_d
645: , WORK, J1 + o_work, KD1);
646: }
647: }
648: else
649: {
650: JEND = J1 + KD1 * (NR - 1);
651: for (JINC = J1; (KD1 >= 0) ? (JINC <= JEND) : (JINC >= JEND); JINC += KD1)
652: {
653: this._drot.Run(KDM1, ref AB, KD+(JINC - KD) * LDAB + o_ab, INCX, ref AB, KD1+(JINC - KD) * LDAB + o_ab, INCX, D[JINC + o_d]
654: , WORK[JINC + o_work]);
655: }
656: }
657: // *
658: }
659: // *
660: if (K > 2)
661: {
662: if (K <= N - I + 1)
663: {
664: // *
665: // * generate plane rotation to annihilate a(i+k-1,i)
666: // * within the band
667: // *
668: this._dlartg.Run(AB[K - 1+I * LDAB + o_ab], AB[K+I * LDAB + o_ab], ref D[I + K - 1 + o_d], ref WORK[I + K - 1 + o_work], ref TEMP);
669: AB[K - 1+I * LDAB + o_ab] = TEMP;
670: // *
671: // * apply rotation from the left
672: // *
673: this._drot.Run(K - 3, ref AB, K - 2+(I + 1) * LDAB + o_ab, LDAB - 1, ref AB, K - 1+(I + 1) * LDAB + o_ab, LDAB - 1, D[I + K - 1 + o_d]
674: , WORK[I + K - 1 + o_work]);
675: }
676: NR = NR + 1;
677: J1 = J1 - KDN - 1;
678: }
679: // *
680: // * apply plane rotations from both sides to diagonal
681: // * blocks
682: // *
683: if (NR > 0)
684: {
685: this._dlar2v.Run(NR, ref AB, 1+(J1 - 1) * LDAB + o_ab, ref AB, 1+J1 * LDAB + o_ab, ref AB, 2+(J1 - 1) * LDAB + o_ab, INCA, D, J1 + o_d
686: , WORK, J1 + o_work, KD1);
687: }
688: // *
689: // * apply plane rotations from the right
690: // *
691: // *
692: // * Dependent on the the number of diagonals either
693: // * DLARTV or DROT is used
694: // *
695: if (NR > 0)
696: {
697: if (NR > 2 * KD - 1)
698: {
699: for (L = 1; L <= KD - 1; L++)
700: {
701: if (J2 + L > N)
702: {
703: NRT = NR - 1;
704: }
705: else
706: {
707: NRT = NR;
708: }
709: if (NRT > 0)
710: {
711: this._dlartv.Run(NRT, ref AB, L + 2+(J1 - 1) * LDAB + o_ab, INCA, ref AB, L + 1+J1 * LDAB + o_ab, INCA, D, J1 + o_d
712: , WORK, J1 + o_work, KD1);
713: }
714: }
715: }
716: else
717: {
718: J1END = J1 + KD1 * (NR - 2);
719: if (J1END >= J1)
720: {
721: for (J1INC = J1; (KD1 >= 0) ? (J1INC <= J1END) : (J1INC >= J1END); J1INC += KD1)
722: {
723: this._drot.Run(KDM1, ref AB, 3+(J1INC - 1) * LDAB + o_ab, 1, ref AB, 2+J1INC * LDAB + o_ab, 1, D[J1INC + o_d]
724: , WORK[J1INC + o_work]);
725: }
726: }
727: LEND = Math.Min(KDM1, N - J2);
728: LAST = J1END + KD1;
729: if (LEND > 0)
730: {
731: this._drot.Run(LEND, ref AB, 3+(LAST - 1) * LDAB + o_ab, 1, ref AB, 2+LAST * LDAB + o_ab, 1, D[LAST + o_d]
732: , WORK[LAST + o_work]);
733: }
734: }
735: }
736: // *
737: // *
738: // *
739: if (WANTQ)
740: {
741: // *
742: // * accumulate product of plane rotations in Q
743: // *
744: if (INITQ)
745: {
746: // *
747: // * take advantage of the fact that Q was
748: // * initially the Identity matrix
749: // *
750: IQEND = Math.Max(IQEND, J2);
751: I2 = Math.Max(0, K - 3);
752: IQAEND = 1 + I * KD;
753: if (K == 2) IQAEND = IQAEND + KD;
754: IQAEND = Math.Min(IQAEND, IQEND);
755: for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
756: {
757: IBL = I - I2 / KDM1;
758: I2 = I2 + 1;
759: IQB = Math.Max(1, J - IBL);
760: NQ = 1 + IQAEND - IQB;
761: IQAEND = Math.Min(IQAEND + KD, IQEND);
762: this._drot.Run(NQ, ref Q, IQB+(J - 1) * LDQ + o_q, 1, ref Q, IQB+J * LDQ + o_q, 1, D[J + o_d]
763: , WORK[J + o_work]);
764: }
765: }
766: else
767: {
768: // *
769: for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
770: {
771: this._drot.Run(N, ref Q, 1+(J - 1) * LDQ + o_q, 1, ref Q, 1+J * LDQ + o_q, 1, D[J + o_d]
772: , WORK[J + o_work]);
773: }
774: }
775: }
776: // *
777: if (J2 + KDN > N)
778: {
779: // *
780: // * adjust J2 to keep within the bounds of the matrix
781: // *
782: NR = NR - 1;
783: J2 = J2 - KDN - 1;
784: }
785: // *
786: for (J = J1; (KD1 >= 0) ? (J <= J2) : (J >= J2); J += KD1)
787: {
788: // *
789: // * create nonzero element a(j+kd,j-1) outside the
790: // * band and store it in WORK
791: // *
792: WORK[J + KD + o_work] = WORK[J + o_work] * AB[KD1+J * LDAB + o_ab];
793: AB[KD1+J * LDAB + o_ab] = D[J + o_d] * AB[KD1+J * LDAB + o_ab];
794: }
795: }
796: }
797: }
798: // *
799: if (KD > 0)
800: {
801: // *
802: // * copy off-diagonal elements to E
803: // *
804: for (I = 1; I <= N - 1; I++)
805: {
806: E[I + o_e] = AB[2+I * LDAB + o_ab];
807: }
808: }
809: else
810: {
811: // *
812: // * set E to zero if original matrix was diagonal
813: // *
814: for (I = 1; I <= N - 1; I++)
815: {
816: E[I + o_e] = ZERO;
817: }
818: }
819: // *
820: // * copy diagonal elements to D
821: // *
822: for (I = 1; I <= N; I++)
823: {
824: D[I + o_d] = AB[1+I * LDAB + o_ab];
825: }
826: }
827: // *
828: return;
829: // *
830: // * End of DSBTRD
831: // *
832:
833: #endregion
834:
835: }
836: }
837: }