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: /// DGEBAL balances a general real matrix A. This involves, first,
27: /// permuting A by a similarity transformation to isolate eigenvalues
28: /// in the first 1 to ILO-1 and last IHI+1 to N elements on the
29: /// diagonal; and second, applying a diagonal similarity transformation
30: /// to rows and columns ILO to IHI to make the rows and columns as
31: /// close in norm as possible. Both steps are optional.
32: ///
33: /// Balancing may reduce the 1-norm of the matrix, and improve the
34: /// accuracy of the computed eigenvalues and/or eigenvectors.
35: ///
36: ///</summary>
37: public class DGEBAL
38: {
39:
40:
41: #region Dependencies
42:
43: LSAME _lsame; IDAMAX _idamax; DLAMCH _dlamch; DSCAL _dscal; DSWAP _dswap; XERBLA _xerbla;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double SCLFAC = 2.0E+0; const double FACTOR = 0.95E+0;
51: bool NOCONV = false;int I = 0; int ICA = 0; int IEXC = 0; int IRA = 0; int J = 0; int K = 0; int L = 0; int M = 0;
52: double C = 0;double CA = 0; double F = 0; double G = 0; double R = 0; double RA = 0; double S = 0; double SFMAX1 = 0;
53: double SFMAX2 = 0;double SFMIN1 = 0; double SFMIN2 = 0;
54:
55: #endregion
56:
57: public DGEBAL(LSAME lsame, IDAMAX idamax, DLAMCH dlamch, DSCAL dscal, DSWAP dswap, XERBLA xerbla)
58: {
59:
60:
61: #region Set Dependencies
62:
63: this._lsame = lsame; this._idamax = idamax; this._dlamch = dlamch; this._dscal = dscal; this._dswap = dswap;
64: this._xerbla = xerbla;
65:
66: #endregion
67:
68: }
69:
70: public DGEBAL()
71: {
72:
73:
74: #region Dependencies (Initialization)
75:
76: LSAME lsame = new LSAME();
77: IDAMAX idamax = new IDAMAX();
78: DLAMC3 dlamc3 = new DLAMC3();
79: DSCAL dscal = new DSCAL();
80: DSWAP dswap = new DSWAP();
81: XERBLA xerbla = new XERBLA();
82: DLAMC1 dlamc1 = new DLAMC1(dlamc3);
83: DLAMC4 dlamc4 = new DLAMC4(dlamc3);
84: DLAMC5 dlamc5 = new DLAMC5(dlamc3);
85: DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
86: DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
87:
88: #endregion
89:
90:
91: #region Set Dependencies
92:
93: this._lsame = lsame; this._idamax = idamax; this._dlamch = dlamch; this._dscal = dscal; this._dswap = dswap;
94: this._xerbla = xerbla;
95:
96: #endregion
97:
98: }
99: /// <summary>
100: /// Purpose
101: /// =======
102: ///
103: /// DGEBAL balances a general real matrix A. This involves, first,
104: /// permuting A by a similarity transformation to isolate eigenvalues
105: /// in the first 1 to ILO-1 and last IHI+1 to N elements on the
106: /// diagonal; and second, applying a diagonal similarity transformation
107: /// to rows and columns ILO to IHI to make the rows and columns as
108: /// close in norm as possible. Both steps are optional.
109: ///
110: /// Balancing may reduce the 1-norm of the matrix, and improve the
111: /// accuracy of the computed eigenvalues and/or eigenvectors.
112: ///
113: ///</summary>
114: /// <param name="JOB">
115: /// (input) CHARACTER*1
116: /// Specifies the operations to be performed on A:
117: /// = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
118: /// for i = 1,...,N;
119: /// = 'P': permute only;
120: /// = 'S': scale only;
121: /// = 'B': both permute and scale.
122: ///</param>
123: /// <param name="N">
124: /// (input) INTEGER
125: /// The order of the matrix A. N .GE. 0.
126: ///</param>
127: /// <param name="A">
128: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
129: /// On entry, the input matrix A.
130: /// On exit, A is overwritten by the balanced matrix.
131: /// If JOB = 'N', A is not referenced.
132: /// 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="ILO">
139: /// (output) INTEGER
140: ///</param>
141: /// <param name="IHI">
142: /// (output) INTEGER
143: /// ILO and IHI are set to integers such that on exit
144: /// A(i,j) = 0 if i .GT. j and j = 1,...,ILO-1 or I = IHI+1,...,N.
145: /// If JOB = 'N' or 'S', ILO = 1 and IHI = N.
146: ///</param>
147: /// <param name="SCALE">
148: /// (output) DOUBLE PRECISION array, dimension (N)
149: /// Details of the permutations and scaling factors applied to
150: /// A. If P(j) is the index of the row and column interchanged
151: /// with row and column j and D(j) is the scaling factor
152: /// applied to row and column j, then
153: /// SCALE(j) = P(j) for j = 1,...,ILO-1
154: /// = D(j) for j = ILO,...,IHI
155: /// = P(j) for j = IHI+1,...,N.
156: /// The order in which the interchanges are made is N to IHI+1,
157: /// then 1 to ILO-1.
158: ///</param>
159: /// <param name="INFO">
160: /// (output) INTEGER
161: /// = 0: successful exit.
162: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value.
163: ///</param>
164: public void Run(string JOB, int N, ref double[] A, int offset_a, int LDA, ref int ILO, ref int IHI
165: , ref double[] SCALE, int offset_scale, ref int INFO)
166: {
167:
168: #region Array Index Correction
169:
170: int o_a = -1 - LDA + offset_a; int o_scale = -1 + offset_scale;
171:
172: #endregion
173:
174:
175: #region Strings
176:
177: JOB = JOB.Substring(0, 1);
178:
179: #endregion
180:
181:
182: #region Prolog
183:
184: // *
185: // * -- LAPACK routine (version 3.1) --
186: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
187: // * November 2006
188: // *
189: // * .. Scalar Arguments ..
190: // * ..
191: // * .. Array Arguments ..
192: // * ..
193: // *
194: // * Purpose
195: // * =======
196: // *
197: // * DGEBAL balances a general real matrix A. This involves, first,
198: // * permuting A by a similarity transformation to isolate eigenvalues
199: // * in the first 1 to ILO-1 and last IHI+1 to N elements on the
200: // * diagonal; and second, applying a diagonal similarity transformation
201: // * to rows and columns ILO to IHI to make the rows and columns as
202: // * close in norm as possible. Both steps are optional.
203: // *
204: // * Balancing may reduce the 1-norm of the matrix, and improve the
205: // * accuracy of the computed eigenvalues and/or eigenvectors.
206: // *
207: // * Arguments
208: // * =========
209: // *
210: // * JOB (input) CHARACTER*1
211: // * Specifies the operations to be performed on A:
212: // * = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
213: // * for i = 1,...,N;
214: // * = 'P': permute only;
215: // * = 'S': scale only;
216: // * = 'B': both permute and scale.
217: // *
218: // * N (input) INTEGER
219: // * The order of the matrix A. N >= 0.
220: // *
221: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
222: // * On entry, the input matrix A.
223: // * On exit, A is overwritten by the balanced matrix.
224: // * If JOB = 'N', A is not referenced.
225: // * See Further Details.
226: // *
227: // * LDA (input) INTEGER
228: // * The leading dimension of the array A. LDA >= max(1,N).
229: // *
230: // * ILO (output) INTEGER
231: // * IHI (output) INTEGER
232: // * ILO and IHI are set to integers such that on exit
233: // * A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
234: // * If JOB = 'N' or 'S', ILO = 1 and IHI = N.
235: // *
236: // * SCALE (output) DOUBLE PRECISION array, dimension (N)
237: // * Details of the permutations and scaling factors applied to
238: // * A. If P(j) is the index of the row and column interchanged
239: // * with row and column j and D(j) is the scaling factor
240: // * applied to row and column j, then
241: // * SCALE(j) = P(j) for j = 1,...,ILO-1
242: // * = D(j) for j = ILO,...,IHI
243: // * = P(j) for j = IHI+1,...,N.
244: // * The order in which the interchanges are made is N to IHI+1,
245: // * then 1 to ILO-1.
246: // *
247: // * INFO (output) INTEGER
248: // * = 0: successful exit.
249: // * < 0: if INFO = -i, the i-th argument had an illegal value.
250: // *
251: // * Further Details
252: // * ===============
253: // *
254: // * The permutations consist of row and column interchanges which put
255: // * the matrix in the form
256: // *
257: // * ( T1 X Y )
258: // * P A P = ( 0 B Z )
259: // * ( 0 0 T2 )
260: // *
261: // * where T1 and T2 are upper triangular matrices whose eigenvalues lie
262: // * along the diagonal. The column indices ILO and IHI mark the starting
263: // * and ending columns of the submatrix B. Balancing consists of applying
264: // * a diagonal similarity transformation inv(D) * B * D to make the
265: // * 1-norms of each row of B and its corresponding column nearly equal.
266: // * The output matrix is
267: // *
268: // * ( T1 X*D Y )
269: // * ( 0 inv(D)*B*D inv(D)*Z ).
270: // * ( 0 0 T2 )
271: // *
272: // * Information about the permutations P and the diagonal matrix D is
273: // * returned in the vector SCALE.
274: // *
275: // * This subroutine is based on the EISPACK routine BALANC.
276: // *
277: // * Modified by Tzu-Yi Chen, Computer Science Division, University of
278: // * California at Berkeley, USA
279: // *
280: // * =====================================================================
281: // *
282: // * .. Parameters ..
283: // * ..
284: // * .. Local Scalars ..
285: // * ..
286: // * .. External Functions ..
287: // * ..
288: // * .. External Subroutines ..
289: // * ..
290: // * .. Intrinsic Functions ..
291: // INTRINSIC ABS, MAX, MIN;
292: // * ..
293: // * .. Executable Statements ..
294: // *
295: // * Test the input parameters
296: // *
297:
298: #endregion
299:
300:
301: #region Body
302:
303: INFO = 0;
304: if (!this._lsame.Run(JOB, "N") && !this._lsame.Run(JOB, "P") && !this._lsame.Run(JOB, "S") && !this._lsame.Run(JOB, "B"))
305: {
306: INFO = - 1;
307: }
308: else
309: {
310: if (N < 0)
311: {
312: INFO = - 2;
313: }
314: else
315: {
316: if (LDA < Math.Max(1, N))
317: {
318: INFO = - 4;
319: }
320: }
321: }
322: if (INFO != 0)
323: {
324: this._xerbla.Run("DGEBAL", - INFO);
325: return;
326: }
327: // *
328: K = 1;
329: L = N;
330: // *
331: if (N == 0) goto LABEL210;
332: // *
333: if (this._lsame.Run(JOB, "N"))
334: {
335: for (I = 1; I <= N; I++)
336: {
337: SCALE[I + o_scale] = ONE;
338: }
339: goto LABEL210;
340: }
341: // *
342: if (this._lsame.Run(JOB, "S")) goto LABEL120;
343: // *
344: // * Permutation to isolate eigenvalues if possible
345: // *
346: goto LABEL50;
347: // *
348: // * Row and column exchange.
349: // *
350: LABEL20:;
351: SCALE[M + o_scale] = J;
352: if (J == M) goto LABEL30;
353: // *
354: this._dswap.Run(L, ref A, 1+J * LDA + o_a, 1, ref A, 1+M * LDA + o_a, 1);
355: this._dswap.Run(N - K + 1, ref A, J+K * LDA + o_a, LDA, ref A, M+K * LDA + o_a, LDA);
356: // *
357: LABEL30:;
358: switch (IEXC)
359: {
360: case 1: goto LABEL40;
361: case 2: goto LABEL80;
362: }
363: // *
364: // * Search for rows isolating an eigenvalue and push them down.
365: // *
366: LABEL40:;
367: if (L == 1) goto LABEL210;
368: L = L - 1;
369: // *
370: LABEL50:;
371: for (J = L; J >= 1; J += - 1)
372: {
373: // *
374: for (I = 1; I <= L; I++)
375: {
376: if (I == J) goto LABEL60;
377: if (A[J+I * LDA + o_a] != ZERO) goto LABEL70;
378: LABEL60:;
379: }
380: // *
381: M = L;
382: IEXC = 1;
383: goto LABEL20;
384: LABEL70:;
385: }
386: // *
387: goto LABEL90;
388: // *
389: // * Search for columns isolating an eigenvalue and push them left.
390: // *
391: LABEL80:;
392: K = K + 1;
393: // *
394: LABEL90:;
395: for (J = K; J <= L; J++)
396: {
397: // *
398: for (I = K; I <= L; I++)
399: {
400: if (I == J) goto LABEL100;
401: if (A[I+J * LDA + o_a] != ZERO) goto LABEL110;
402: LABEL100:;
403: }
404: // *
405: M = K;
406: IEXC = 2;
407: goto LABEL20;
408: LABEL110:;
409: }
410: // *
411: LABEL120:;
412: for (I = K; I <= L; I++)
413: {
414: SCALE[I + o_scale] = ONE;
415: }
416: // *
417: if (this._lsame.Run(JOB, "P")) goto LABEL210;
418: // *
419: // * Balance the submatrix in rows K to L.
420: // *
421: // * Iterative loop for norm reduction
422: // *
423: SFMIN1 = this._dlamch.Run("S") / this._dlamch.Run("P");
424: SFMAX1 = ONE / SFMIN1;
425: SFMIN2 = SFMIN1 * SCLFAC;
426: SFMAX2 = ONE / SFMIN2;
427: LABEL140:;
428: NOCONV = false;
429: // *
430: for (I = K; I <= L; I++)
431: {
432: C = ZERO;
433: R = ZERO;
434: // *
435: for (J = K; J <= L; J++)
436: {
437: if (J == I) goto LABEL150;
438: C = C + Math.Abs(A[J+I * LDA + o_a]);
439: R = R + Math.Abs(A[I+J * LDA + o_a]);
440: LABEL150:;
441: }
442: ICA = this._idamax.Run(L, A, 1+I * LDA + o_a, 1);
443: CA = Math.Abs(A[ICA+I * LDA + o_a]);
444: IRA = this._idamax.Run(N - K + 1, A, I+K * LDA + o_a, LDA);
445: RA = Math.Abs(A[I+(IRA + K - 1) * LDA + o_a]);
446: // *
447: // * Guard against zero C or R due to underflow.
448: // *
449: if (C == ZERO || R == ZERO) goto LABEL200;
450: G = R / SCLFAC;
451: F = ONE;
452: S = C + R;
453: LABEL160:;
454: if (C >= G || Math.Max(F, Math.Max(C, CA)) >= SFMAX2 || Math.Min(R, Math.Min(G, RA)) <= SFMIN2) goto LABEL170;
455: F = F * SCLFAC;
456: C = C * SCLFAC;
457: CA = CA * SCLFAC;
458: R = R / SCLFAC;
459: G = G / SCLFAC;
460: RA = RA / SCLFAC;
461: goto LABEL160;
462: // *
463: LABEL170:;
464: G = C / SCLFAC;
465: LABEL180:;
466: if (G < R || Math.Max(R, RA) >= SFMAX2 || Math.Min(F, Math.Min(C, Math.Min(G, CA))) <= SFMIN2) goto LABEL190;
467: F = F / SCLFAC;
468: C = C / SCLFAC;
469: G = G / SCLFAC;
470: CA = CA / SCLFAC;
471: R = R * SCLFAC;
472: RA = RA * SCLFAC;
473: goto LABEL180;
474: // *
475: // * Now balance.
476: // *
477: LABEL190:;
478: if ((C + R) >= FACTOR * S) goto LABEL200;
479: if (F < ONE && SCALE[I + o_scale] < ONE)
480: {
481: if (F * SCALE[I + o_scale] <= SFMIN1) goto LABEL200;
482: }
483: if (F > ONE && SCALE[I + o_scale] > ONE)
484: {
485: if (SCALE[I + o_scale] >= SFMAX1 / F) goto LABEL200;
486: }
487: G = ONE / F;
488: SCALE[I + o_scale] = SCALE[I + o_scale] * F;
489: NOCONV = true;
490: // *
491: this._dscal.Run(N - K + 1, G, ref A, I+K * LDA + o_a, LDA);
492: this._dscal.Run(L, F, ref A, 1+I * LDA + o_a, 1);
493: // *
494: LABEL200:;
495: }
496: // *
497: if (NOCONV) goto LABEL140;
498: // *
499: LABEL210:;
500: ILO = K;
501: IHI = L;
502: // *
503: return;
504: // *
505: // * End of DGEBAL
506: // *
507:
508: #endregion
509:
510: }
511: }
512: }