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: /// DTREXC reorders the real Schur factorization of a real matrix
27: /// A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
28: /// moved to row ILST.
29: ///
30: /// The real Schur form T is reordered by an orthogonal similarity
31: /// transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
32: /// is updated by postmultiplying it with Z.
33: ///
34: /// T must be in Schur canonical form (as returned by DHSEQR), that is,
35: /// block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
36: /// 2-by-2 diagonal block has its diagonal elements equal and its
37: /// off-diagonal elements of opposite sign.
38: ///
39: ///</summary>
40: public class DTREXC
41: {
42:
43:
44: #region Dependencies
45:
46: LSAME _lsame; DLAEXC _dlaexc; XERBLA _xerbla;
47:
48: #endregion
49:
50:
51: #region Fields
52:
53: const double ZERO = 0.0E+0; bool WANTQ = false; int HERE = 0; int NBF = 0; int NBL = 0; int NBNEXT = 0;
54:
55: #endregion
56:
57: public DTREXC(LSAME lsame, DLAEXC dlaexc, XERBLA xerbla)
58: {
59:
60:
61: #region Set Dependencies
62:
63: this._lsame = lsame; this._dlaexc = dlaexc; this._xerbla = xerbla;
64:
65: #endregion
66:
67: }
68:
69: public DTREXC()
70: {
71:
72:
73: #region Dependencies (Initialization)
74:
75: LSAME lsame = new LSAME();
76: DLAMC3 dlamc3 = new DLAMC3();
77: DLASSQ dlassq = new DLASSQ();
78: DLAPY2 dlapy2 = new DLAPY2();
79: DNRM2 dnrm2 = new DNRM2();
80: DSCAL dscal = new DSCAL();
81: XERBLA xerbla = new XERBLA();
82: IDAMAX idamax = new IDAMAX();
83: DCOPY dcopy = new DCOPY();
84: DSWAP dswap = new DSWAP();
85: DROT drot = new DROT();
86: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
87: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
88: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
89: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
90: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
91: DLANGE dlange = new DLANGE(dlassq, lsame);
92: DLACPY dlacpy = new DLACPY(lsame);
93: DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
94: DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
95: DGEMV dgemv = new DGEMV(lsame, xerbla);
96: DGER dger = new DGER(xerbla);
97: DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
98: DLARTG dlartg = new DLARTG(dlamch);
99: DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
100: DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
101:
102: #endregion
103:
104:
105: #region Set Dependencies
106:
107: this._lsame = lsame; this._dlaexc = dlaexc; this._xerbla = xerbla;
108:
109: #endregion
110:
111: }
112: /// <summary>
113: /// Purpose
114: /// =======
115: ///
116: /// DTREXC reorders the real Schur factorization of a real matrix
117: /// A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
118: /// moved to row ILST.
119: ///
120: /// The real Schur form T is reordered by an orthogonal similarity
121: /// transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
122: /// is updated by postmultiplying it with Z.
123: ///
124: /// T must be in Schur canonical form (as returned by DHSEQR), that is,
125: /// block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
126: /// 2-by-2 diagonal block has its diagonal elements equal and its
127: /// off-diagonal elements of opposite sign.
128: ///
129: ///</summary>
130: /// <param name="COMPQ">
131: /// (input) CHARACTER*1
132: /// = 'V': update the matrix Q of Schur vectors;
133: /// = 'N': do not update Q.
134: ///</param>
135: /// <param name="N">
136: /// (input) INTEGER
137: /// The order of the matrix T. N .GE. 0.
138: ///</param>
139: /// <param name="T">
140: /// must be in Schur canonical form (as returned by DHSEQR), that is,
141: ///</param>
142: /// <param name="LDT">
143: /// (input) INTEGER
144: /// The leading dimension of the array T. LDT .GE. max(1,N).
145: ///</param>
146: /// <param name="Q">
147: /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
148: /// On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
149: /// On exit, if COMPQ = 'V', Q has been postmultiplied by the
150: /// orthogonal transformation matrix Z which reorders T.
151: /// If COMPQ = 'N', Q is not referenced.
152: ///</param>
153: /// <param name="LDQ">
154: /// (input) INTEGER
155: /// The leading dimension of the array Q. LDQ .GE. max(1,N).
156: ///</param>
157: /// <param name="IFST">
158: /// (input/output) INTEGER
159: ///</param>
160: /// <param name="ILST">
161: /// (input/output) INTEGER
162: /// Specify the reordering of the diagonal blocks of T.
163: /// The block with row index IFST is moved to row ILST, by a
164: /// sequence of transpositions between adjacent blocks.
165: /// On exit, if IFST pointed on entry to the second row of a
166: /// 2-by-2 block, it is changed to point to the first row; ILST
167: /// always points to the first row of the block in its final
168: /// position (which may differ from its input value by +1 or -1).
169: /// 1 .LE. IFST .LE. N; 1 .LE. ILST .LE. N.
170: ///</param>
171: /// <param name="WORK">
172: /// (workspace) DOUBLE PRECISION array, dimension (N)
173: ///</param>
174: /// <param name="INFO">
175: /// (output) INTEGER
176: /// = 0: successful exit
177: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
178: /// = 1: two adjacent blocks were too close to swap (the problem
179: /// is very ill-conditioned); T may have been partially
180: /// reordered, and ILST points to the first row of the
181: /// current position of the block being moved.
182: ///</param>
183: public void Run(string COMPQ, int N, ref double[] T, int offset_t, int LDT, ref double[] Q, int offset_q, int LDQ
184: , ref int IFST, ref int ILST, ref double[] WORK, int offset_work, ref int INFO)
185: {
186:
187: #region Array Index Correction
188:
189: int o_t = -1 - LDT + offset_t; int o_q = -1 - LDQ + offset_q; int o_work = -1 + offset_work;
190:
191: #endregion
192:
193:
194: #region Strings
195:
196: COMPQ = COMPQ.Substring(0, 1);
197:
198: #endregion
199:
200:
201: #region Prolog
202:
203: // *
204: // * -- LAPACK routine (version 3.1) --
205: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
206: // * November 2006
207: // *
208: // * .. Scalar Arguments ..
209: // * ..
210: // * .. Array Arguments ..
211: // * ..
212: // *
213: // * Purpose
214: // * =======
215: // *
216: // * DTREXC reorders the real Schur factorization of a real matrix
217: // * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
218: // * moved to row ILST.
219: // *
220: // * The real Schur form T is reordered by an orthogonal similarity
221: // * transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
222: // * is updated by postmultiplying it with Z.
223: // *
224: // * T must be in Schur canonical form (as returned by DHSEQR), that is,
225: // * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
226: // * 2-by-2 diagonal block has its diagonal elements equal and its
227: // * off-diagonal elements of opposite sign.
228: // *
229: // * Arguments
230: // * =========
231: // *
232: // * COMPQ (input) CHARACTER*1
233: // * = 'V': update the matrix Q of Schur vectors;
234: // * = 'N': do not update Q.
235: // *
236: // * N (input) INTEGER
237: // * The order of the matrix T. N >= 0.
238: // *
239: // * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
240: // * On entry, the upper quasi-triangular matrix T, in Schur
241: // * Schur canonical form.
242: // * On exit, the reordered upper quasi-triangular matrix, again
243: // * in Schur canonical form.
244: // *
245: // * LDT (input) INTEGER
246: // * The leading dimension of the array T. LDT >= max(1,N).
247: // *
248: // * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
249: // * On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
250: // * On exit, if COMPQ = 'V', Q has been postmultiplied by the
251: // * orthogonal transformation matrix Z which reorders T.
252: // * If COMPQ = 'N', Q is not referenced.
253: // *
254: // * LDQ (input) INTEGER
255: // * The leading dimension of the array Q. LDQ >= max(1,N).
256: // *
257: // * IFST (input/output) INTEGER
258: // * ILST (input/output) INTEGER
259: // * Specify the reordering of the diagonal blocks of T.
260: // * The block with row index IFST is moved to row ILST, by a
261: // * sequence of transpositions between adjacent blocks.
262: // * On exit, if IFST pointed on entry to the second row of a
263: // * 2-by-2 block, it is changed to point to the first row; ILST
264: // * always points to the first row of the block in its final
265: // * position (which may differ from its input value by +1 or -1).
266: // * 1 <= IFST <= N; 1 <= ILST <= N.
267: // *
268: // * WORK (workspace) DOUBLE PRECISION array, dimension (N)
269: // *
270: // * INFO (output) INTEGER
271: // * = 0: successful exit
272: // * < 0: if INFO = -i, the i-th argument had an illegal value
273: // * = 1: two adjacent blocks were too close to swap (the problem
274: // * is very ill-conditioned); T may have been partially
275: // * reordered, and ILST points to the first row of the
276: // * current position of the block being moved.
277: // *
278: // * =====================================================================
279: // *
280: // * .. Parameters ..
281: // * ..
282: // * .. Local Scalars ..
283: // * ..
284: // * .. External Functions ..
285: // * ..
286: // * .. External Subroutines ..
287: // * ..
288: // * .. Intrinsic Functions ..
289: // INTRINSIC MAX;
290: // * ..
291: // * .. Executable Statements ..
292: // *
293: // * Decode and test the input arguments.
294: // *
295:
296: #endregion
297:
298:
299: #region Body
300:
301: INFO = 0;
302: WANTQ = this._lsame.Run(COMPQ, "V");
303: if (!WANTQ && !this._lsame.Run(COMPQ, "N"))
304: {
305: INFO = - 1;
306: }
307: else
308: {
309: if (N < 0)
310: {
311: INFO = - 2;
312: }
313: else
314: {
315: if (LDT < Math.Max(1, N))
316: {
317: INFO = - 4;
318: }
319: else
320: {
321: if (LDQ < 1 || (WANTQ && LDQ < Math.Max(1, N)))
322: {
323: INFO = - 6;
324: }
325: else
326: {
327: if (IFST < 1 || IFST > N)
328: {
329: INFO = - 7;
330: }
331: else
332: {
333: if (ILST < 1 || ILST > N)
334: {
335: INFO = - 8;
336: }
337: }
338: }
339: }
340: }
341: }
342: if (INFO != 0)
343: {
344: this._xerbla.Run("DTREXC", - INFO);
345: return;
346: }
347: // *
348: // * Quick return if possible
349: // *
350: if (N <= 1) return;
351: // *
352: // * Determine the first row of specified block
353: // * and find out it is 1 by 1 or 2 by 2.
354: // *
355: if (IFST > 1)
356: {
357: if (T[IFST+(IFST - 1) * LDT + o_t] != ZERO) IFST = IFST - 1;
358: }
359: NBF = 1;
360: if (IFST < N)
361: {
362: if (T[IFST + 1+IFST * LDT + o_t] != ZERO) NBF = 2;
363: }
364: // *
365: // * Determine the first row of the final block
366: // * and find out it is 1 by 1 or 2 by 2.
367: // *
368: if (ILST > 1)
369: {
370: if (T[ILST+(ILST - 1) * LDT + o_t] != ZERO) ILST = ILST - 1;
371: }
372: NBL = 1;
373: if (ILST < N)
374: {
375: if (T[ILST + 1+ILST * LDT + o_t] != ZERO) NBL = 2;
376: }
377: // *
378: if (IFST == ILST) return;
379: // *
380: if (IFST < ILST)
381: {
382: // *
383: // * Update ILST
384: // *
385: if (NBF == 2 && NBL == 1) ILST = ILST - 1;
386: if (NBF == 1 && NBL == 2) ILST = ILST + 1;
387: // *
388: HERE = IFST;
389: // *
390: LABEL10:;
391: // *
392: // * Swap block with next one below
393: // *
394: if (NBF == 1 || NBF == 2)
395: {
396: // *
397: // * Current block either 1 by 1 or 2 by 2
398: // *
399: NBNEXT = 1;
400: if (HERE + NBF + 1 <= N)
401: {
402: if (T[HERE + NBF + 1+(HERE + NBF) * LDT + o_t] != ZERO) NBNEXT = 2;
403: }
404: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
405: , HERE, NBF, NBNEXT, ref WORK, offset_work, ref INFO);
406: if (INFO != 0)
407: {
408: ILST = HERE;
409: return;
410: }
411: HERE = HERE + NBNEXT;
412: // *
413: // * Test if 2 by 2 block breaks into two 1 by 1 blocks
414: // *
415: if (NBF == 2)
416: {
417: if (T[HERE + 1+HERE * LDT + o_t] == ZERO) NBF = 3;
418: }
419: // *
420: }
421: else
422: {
423: // *
424: // * Current block consists of two 1 by 1 blocks each of which
425: // * must be swapped individually
426: // *
427: NBNEXT = 1;
428: if (HERE + 3 <= N)
429: {
430: if (T[HERE + 3+(HERE + 2) * LDT + o_t] != ZERO) NBNEXT = 2;
431: }
432: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
433: , HERE + 1, 1, NBNEXT, ref WORK, offset_work, ref INFO);
434: if (INFO != 0)
435: {
436: ILST = HERE;
437: return;
438: }
439: if (NBNEXT == 1)
440: {
441: // *
442: // * Swap two 1 by 1 blocks, no problems possible
443: // *
444: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
445: , HERE, 1, NBNEXT, ref WORK, offset_work, ref INFO);
446: HERE = HERE + 1;
447: }
448: else
449: {
450: // *
451: // * Recompute NBNEXT in case 2 by 2 split
452: // *
453: if (T[HERE + 2+(HERE + 1) * LDT + o_t] == ZERO) NBNEXT = 1;
454: if (NBNEXT == 2)
455: {
456: // *
457: // * 2 by 2 Block did not split
458: // *
459: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
460: , HERE, 1, NBNEXT, ref WORK, offset_work, ref INFO);
461: if (INFO != 0)
462: {
463: ILST = HERE;
464: return;
465: }
466: HERE = HERE + 2;
467: }
468: else
469: {
470: // *
471: // * 2 by 2 Block did split
472: // *
473: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
474: , HERE, 1, 1, ref WORK, offset_work, ref INFO);
475: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
476: , HERE + 1, 1, 1, ref WORK, offset_work, ref INFO);
477: HERE = HERE + 2;
478: }
479: }
480: }
481: if (HERE < ILST) goto LABEL10;
482: // *
483: }
484: else
485: {
486: // *
487: HERE = IFST;
488: LABEL20:;
489: // *
490: // * Swap block with next one above
491: // *
492: if (NBF == 1 || NBF == 2)
493: {
494: // *
495: // * Current block either 1 by 1 or 2 by 2
496: // *
497: NBNEXT = 1;
498: if (HERE >= 3)
499: {
500: if (T[HERE - 1+(HERE - 2) * LDT + o_t] != ZERO) NBNEXT = 2;
501: }
502: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
503: , HERE - NBNEXT, NBNEXT, NBF, ref WORK, offset_work, ref INFO);
504: if (INFO != 0)
505: {
506: ILST = HERE;
507: return;
508: }
509: HERE = HERE - NBNEXT;
510: // *
511: // * Test if 2 by 2 block breaks into two 1 by 1 blocks
512: // *
513: if (NBF == 2)
514: {
515: if (T[HERE + 1+HERE * LDT + o_t] == ZERO) NBF = 3;
516: }
517: // *
518: }
519: else
520: {
521: // *
522: // * Current block consists of two 1 by 1 blocks each of which
523: // * must be swapped individually
524: // *
525: NBNEXT = 1;
526: if (HERE >= 3)
527: {
528: if (T[HERE - 1+(HERE - 2) * LDT + o_t] != ZERO) NBNEXT = 2;
529: }
530: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
531: , HERE - NBNEXT, NBNEXT, 1, ref WORK, offset_work, ref INFO);
532: if (INFO != 0)
533: {
534: ILST = HERE;
535: return;
536: }
537: if (NBNEXT == 1)
538: {
539: // *
540: // * Swap two 1 by 1 blocks, no problems possible
541: // *
542: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
543: , HERE, NBNEXT, 1, ref WORK, offset_work, ref INFO);
544: HERE = HERE - 1;
545: }
546: else
547: {
548: // *
549: // * Recompute NBNEXT in case 2 by 2 split
550: // *
551: if (T[HERE+(HERE - 1) * LDT + o_t] == ZERO) NBNEXT = 1;
552: if (NBNEXT == 2)
553: {
554: // *
555: // * 2 by 2 Block did not split
556: // *
557: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
558: , HERE - 1, 2, 1, ref WORK, offset_work, ref INFO);
559: if (INFO != 0)
560: {
561: ILST = HERE;
562: return;
563: }
564: HERE = HERE - 2;
565: }
566: else
567: {
568: // *
569: // * 2 by 2 Block did split
570: // *
571: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
572: , HERE, 1, 1, ref WORK, offset_work, ref INFO);
573: this._dlaexc.Run(WANTQ, N, ref T, offset_t, LDT, ref Q, offset_q, LDQ
574: , HERE - 1, 1, 1, ref WORK, offset_work, ref INFO);
575: HERE = HERE - 2;
576: }
577: }
578: }
579: if (HERE > ILST) goto LABEL20;
580: }
581: ILST = HERE;
582: // *
583: return;
584: // *
585: // * End of DTREXC
586: // *
587:
588: #endregion
589:
590: }
591: }
592: }