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.0) --
21: /// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
22: /// Courant Institute, Argonne National Lab, and Rice University
23: /// February 29, 1992
24: /// Purpose
25: /// =======
26: ///
27: /// DLACON estimates the 1-norm of a square, real matrix A.
28: /// Reverse communication is used for evaluating matrix-vector products.
29: ///
30: ///</summary>
31: public class DLACON
32: {
33:
34:
35: #region Dependencies
36:
37: IDAMAX _idamax; DASUM _dasum; DCOPY _dcopy;
38:
39: #endregion
40:
41:
42: #region Fields
43:
44: const int ITMAX = 5; const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double TWO = 2.0E+0; int I = 0;
45: int ITER = 0;int J = 0; int JLAST = 0; int JUMP = 0; double ALTSGN = 0; double ESTOLD = 0; double TEMP = 0;
46:
47: #endregion
48:
49: public DLACON(IDAMAX idamax, DASUM dasum, DCOPY dcopy)
50: {
51:
52:
53: #region Set Dependencies
54:
55: this._idamax = idamax; this._dasum = dasum; this._dcopy = dcopy;
56:
57: #endregion
58:
59: }
60:
61: public DLACON()
62: {
63:
64:
65: #region Dependencies (Initialization)
66:
67: IDAMAX idamax = new IDAMAX();
68: DASUM dasum = new DASUM();
69: DCOPY dcopy = new DCOPY();
70:
71: #endregion
72:
73:
74: #region Set Dependencies
75:
76: this._idamax = idamax; this._dasum = dasum; this._dcopy = dcopy;
77:
78: #endregion
79:
80: }
81: /// <summary>
82: /// Purpose
83: /// =======
84: ///
85: /// DLACON estimates the 1-norm of a square, real matrix A.
86: /// Reverse communication is used for evaluating matrix-vector products.
87: ///
88: ///</summary>
89: /// <param name="N">
90: /// (input) INTEGER
91: /// The order of the matrix. N .GE. 1.
92: ///</param>
93: /// <param name="V">
94: /// (workspace) DOUBLE PRECISION array, dimension (N)
95: /// On the final return, V = A*W, where EST = norm(V)/norm(W)
96: /// (W is not returned).
97: ///</param>
98: /// <param name="X">
99: /// (input/output) DOUBLE PRECISION array, dimension (N)
100: /// On an intermediate return, X should be overwritten by
101: /// A * X, if KASE=1,
102: /// A' * X, if KASE=2,
103: /// and DLACON must be re-called with all the other parameters
104: /// unchanged.
105: ///</param>
106: /// <param name="ISGN">
107: /// (workspace) INTEGER array, dimension (N)
108: ///</param>
109: /// <param name="EST">
110: /// (output) DOUBLE PRECISION
111: /// An estimate (a lower bound) for norm(A).
112: ///</param>
113: /// <param name="KASE">
114: /// (input/output) INTEGER
115: /// On the initial call to DLACON, KASE should be 0.
116: /// On an intermediate return, KASE will be 1 or 2, indicating
117: /// whether X should be overwritten by A * X or A' * X.
118: /// On the final return from DLACON, KASE will again be 0.
119: ///</param>
120: public void Run(int N, ref double[] V, int offset_v, ref double[] X, int offset_x, ref int[] ISGN, int offset_isgn, ref double EST, ref int KASE)
121: {
122:
123: #region Array Index Correction
124:
125: int o_v = -1 + offset_v; int o_x = -1 + offset_x; int o_isgn = -1 + offset_isgn;
126:
127: #endregion
128:
129:
130: #region Prolog
131:
132: // *
133: // * -- LAPACK auxiliary routine (version 3.0) --
134: // * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
135: // * Courant Institute, Argonne National Lab, and Rice University
136: // * February 29, 1992
137: // *
138: // * .. Scalar Arguments ..
139: // * ..
140: // * .. Array Arguments ..
141: // * ..
142: // *
143: // * Purpose
144: // * =======
145: // *
146: // * DLACON estimates the 1-norm of a square, real matrix A.
147: // * Reverse communication is used for evaluating matrix-vector products.
148: // *
149: // * Arguments
150: // * =========
151: // *
152: // * N (input) INTEGER
153: // * The order of the matrix. N >= 1.
154: // *
155: // * V (workspace) DOUBLE PRECISION array, dimension (N)
156: // * On the final return, V = A*W, where EST = norm(V)/norm(W)
157: // * (W is not returned).
158: // *
159: // * X (input/output) DOUBLE PRECISION array, dimension (N)
160: // * On an intermediate return, X should be overwritten by
161: // * A * X, if KASE=1,
162: // * A' * X, if KASE=2,
163: // * and DLACON must be re-called with all the other parameters
164: // * unchanged.
165: // *
166: // * ISGN (workspace) INTEGER array, dimension (N)
167: // *
168: // * EST (output) DOUBLE PRECISION
169: // * An estimate (a lower bound) for norm(A).
170: // *
171: // * KASE (input/output) INTEGER
172: // * On the initial call to DLACON, KASE should be 0.
173: // * On an intermediate return, KASE will be 1 or 2, indicating
174: // * whether X should be overwritten by A * X or A' * X.
175: // * On the final return from DLACON, KASE will again be 0.
176: // *
177: // * Further Details
178: // * ======= =======
179: // *
180: // * Contributed by Nick Higham, University of Manchester.
181: // * Originally named SONEST, dated March 16, 1988.
182: // *
183: // * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
184: // * a real or complex matrix, with applications to condition estimation",
185: // * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
186: // *
187: // * =====================================================================
188: // *
189: // * .. Parameters ..
190: // * ..
191: // * .. Local Scalars ..
192: // * ..
193: // * .. External Functions ..
194: // * ..
195: // * .. External Subroutines ..
196: // * ..
197: // * .. Intrinsic Functions ..
198: // INTRINSIC ABS, DBLE, NINT, SIGN;
199: // * ..
200: // * .. Save statement ..
201: // * ..
202: // * .. Executable Statements ..
203: // *
204:
205: #endregion
206:
207:
208: #region Body
209:
210: if (KASE == 0)
211: {
212: for (I = 1; I <= N; I++)
213: {
214: X[I + o_x] = ONE / Convert.ToDouble(N);
215: }
216: KASE = 1;
217: JUMP = 1;
218: return;
219: }
220: // *
221: switch (JUMP)
222: {
223: case 1: goto LABEL20;
224: case 2: goto LABEL40;
225: case 3: goto LABEL70;
226: case 4: goto LABEL110;
227: case 5: goto LABEL140;
228: }
229: // *
230: // * ................ ENTRY (JUMP = 1)
231: // * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
232: // *
233: LABEL20:;
234: if (N == 1)
235: {
236: V[1 + o_v] = X[1 + o_x];
237: EST = Math.Abs(V[1 + o_v]);
238: // * ... QUIT
239: goto LABEL150;
240: }
241: EST = this._dasum.Run(N, X, offset_x, 1);
242: // *
243: for (I = 1; I <= N; I++)
244: {
245: X[I + o_x] = FortranLib.Sign(ONE,X[I + o_x]);
246: ISGN[I + o_isgn] = (int)Math.Round(X[I + o_x]);
247: }
248: KASE = 2;
249: JUMP = 2;
250: return;
251: // *
252: // * ................ ENTRY (JUMP = 2)
253: // * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
254: // *
255: LABEL40:;
256: J = this._idamax.Run(N, X, offset_x, 1);
257: ITER = 2;
258: // *
259: // * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
260: // *
261: LABEL50:;
262: for (I = 1; I <= N; I++)
263: {
264: X[I + o_x] = ZERO;
265: }
266: X[J + o_x] = ONE;
267: KASE = 1;
268: JUMP = 3;
269: return;
270: // *
271: // * ................ ENTRY (JUMP = 3)
272: // * X HAS BEEN OVERWRITTEN BY A*X.
273: // *
274: LABEL70:;
275: this._dcopy.Run(N, X, offset_x, 1, ref V, offset_v, 1);
276: ESTOLD = EST;
277: EST = this._dasum.Run(N, V, offset_v, 1);
278: for (I = 1; I <= N; I++)
279: {
280: if (Math.Round(FortranLib.Sign(ONE,X[I + o_x])) != ISGN[I + o_isgn]) goto LABEL90;
281: }
282: // * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
283: goto LABEL120;
284: // *
285: LABEL90:;
286: // * TEST FOR CYCLING.
287: if (EST <= ESTOLD) goto LABEL120;
288: // *
289: for (I = 1; I <= N; I++)
290: {
291: X[I + o_x] = FortranLib.Sign(ONE,X[I + o_x]);
292: ISGN[I + o_isgn] = (int)Math.Round(X[I + o_x]);
293: }
294: KASE = 2;
295: JUMP = 4;
296: return;
297: // *
298: // * ................ ENTRY (JUMP = 4)
299: // * X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
300: // *
301: LABEL110:;
302: JLAST = J;
303: J = this._idamax.Run(N, X, offset_x, 1);
304: if ((X[JLAST + o_x] != Math.Abs(X[J + o_x])) && (ITER < ITMAX))
305: {
306: ITER = ITER + 1;
307: goto LABEL50;
308: }
309: // *
310: // * ITERATION COMPLETE. FINAL STAGE.
311: // *
312: LABEL120:;
313: ALTSGN = ONE;
314: for (I = 1; I <= N; I++)
315: {
316: X[I + o_x] = ALTSGN * (ONE + Convert.ToDouble(I - 1) / Convert.ToDouble(N - 1));
317: ALTSGN = - ALTSGN;
318: }
319: KASE = 1;
320: JUMP = 5;
321: return;
322: // *
323: // * ................ ENTRY (JUMP = 5)
324: // * X HAS BEEN OVERWRITTEN BY A*X.
325: // *
326: LABEL140:;
327: TEMP = TWO * (this._dasum.Run(N, X, offset_x, 1) / Convert.ToDouble(3 * N));
328: if (TEMP > EST)
329: {
330: this._dcopy.Run(N, X, offset_x, 1, ref V, offset_v, 1);
331: EST = TEMP;
332: }
333: // *
334: LABEL150:;
335: KASE = 0;
336: return;
337: // *
338: // * End of DLACON
339: // *
340:
341: #endregion
342:
343: }
344: }
345: }