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.0) --
21: /// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
22: /// Courant Institute, Argonne National Lab, and Rice University
23: /// March 31, 1993
24: /// Purpose
25: /// =======
26: ///
27: /// DTRCON estimates the reciprocal of the condition number of a
28: /// triangular matrix A, in either the 1-norm or the infinity-norm.
29: ///
30: /// The norm of A is computed and an estimate is obtained for
31: /// norm(inv(A)), then the reciprocal of the condition number is
32: /// computed as
33: /// RCOND = 1 / ( norm(A) * norm(inv(A)) ).
34: ///
35: ///</summary>
36: public class DTRCON
37: {
38:
39:
40: #region Dependencies
41:
42: LSAME _lsame; IDAMAX _idamax; DLAMCH _dlamch; DLANTR _dlantr; DLACON _dlacon; DLATRS _dlatrs; DRSCL _drscl;
43: XERBLA _xerbla;
44:
45: #endregion
46:
47:
48: #region Fields
49:
50: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; bool NOUNIT = false; bool ONENRM = false; bool UPPER = false;
51: string NORMIN = new string(' ', 1);int IX = 0; int KASE = 0; int KASE1 = 0; double AINVNM = 0; double ANORM = 0;
52: double SCALE = 0;double SMLNUM = 0; double XNORM = 0;
53:
54: #endregion
55:
56: public DTRCON(LSAME lsame, IDAMAX idamax, DLAMCH dlamch, DLANTR dlantr, DLACON dlacon, DLATRS dlatrs, DRSCL drscl, XERBLA xerbla)
57: {
58:
59:
60: #region Set Dependencies
61:
62: this._lsame = lsame; this._idamax = idamax; this._dlamch = dlamch; this._dlantr = dlantr; this._dlacon = dlacon;
63: this._dlatrs = dlatrs;this._drscl = drscl; this._xerbla = xerbla;
64:
65: #endregion
66:
67: }
68:
69: public DTRCON()
70: {
71:
72:
73: #region Dependencies (Initialization)
74:
75: LSAME lsame = new LSAME();
76: IDAMAX idamax = new IDAMAX();
77: DLAMC3 dlamc3 = new DLAMC3();
78: DLASSQ dlassq = new DLASSQ();
79: DASUM dasum = new DASUM();
80: DCOPY dcopy = new DCOPY();
81: DDOT ddot = new DDOT();
82: DAXPY daxpy = new DAXPY();
83: DSCAL dscal = new DSCAL();
84: XERBLA xerbla = new XERBLA();
85: DLABAD dlabad = new DLABAD();
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: DLANTR dlantr = new DLANTR(dlassq, lsame);
92: DLACON dlacon = new DLACON(idamax, dasum, dcopy);
93: DTRSV dtrsv = new DTRSV(lsame, xerbla);
94: DLATRS dlatrs = new DLATRS(lsame, idamax, dasum, ddot, dlamch, daxpy, dscal, dtrsv, xerbla);
95: DRSCL drscl = new DRSCL(dlamch, dscal, dlabad);
96:
97: #endregion
98:
99:
100: #region Set Dependencies
101:
102: this._lsame = lsame; this._idamax = idamax; this._dlamch = dlamch; this._dlantr = dlantr; this._dlacon = dlacon;
103: this._dlatrs = dlatrs;this._drscl = drscl; this._xerbla = xerbla;
104:
105: #endregion
106:
107: }
108: /// <summary>
109: /// Purpose
110: /// =======
111: ///
112: /// DTRCON estimates the reciprocal of the condition number of a
113: /// triangular matrix A, in either the 1-norm or the infinity-norm.
114: ///
115: /// The norm of A is computed and an estimate is obtained for
116: /// norm(inv(A)), then the reciprocal of the condition number is
117: /// computed as
118: /// RCOND = 1 / ( norm(A) * norm(inv(A)) ).
119: ///
120: ///</summary>
121: /// <param name="NORM">
122: /// (input) CHARACTER*1
123: /// Specifies whether the 1-norm condition number or the
124: /// infinity-norm condition number is required:
125: /// = '1' or 'O': 1-norm;
126: /// = 'I': Infinity-norm.
127: ///</param>
128: /// <param name="UPLO">
129: /// (input) CHARACTER*1
130: /// = 'U': A is upper triangular;
131: /// = 'L': A is lower triangular.
132: ///</param>
133: /// <param name="DIAG">
134: /// (input) CHARACTER*1
135: /// = 'N': A is non-unit triangular;
136: /// = 'U': A is unit triangular.
137: ///</param>
138: /// <param name="N">
139: /// (input) INTEGER
140: /// The order of the matrix A. N .GE. 0.
141: ///</param>
142: /// <param name="A">
143: /// (input) DOUBLE PRECISION array, dimension (LDA,N)
144: /// The triangular matrix A. If UPLO = 'U', the leading N-by-N
145: /// upper triangular part of the array A contains the upper
146: /// triangular matrix, and the strictly lower triangular part of
147: /// A is not referenced. If UPLO = 'L', the leading N-by-N lower
148: /// triangular part of the array A contains the lower triangular
149: /// matrix, and the strictly upper triangular part of A is not
150: /// referenced. If DIAG = 'U', the diagonal elements of A are
151: /// also not referenced and are assumed to be 1.
152: ///</param>
153: /// <param name="LDA">
154: /// (input) INTEGER
155: /// The leading dimension of the array A. LDA .GE. max(1,N).
156: ///</param>
157: /// <param name="RCOND">
158: /// (output) DOUBLE PRECISION
159: /// The reciprocal of the condition number of the matrix A,
160: /// computed as RCOND = 1/(norm(A) * norm(inv(A))).
161: ///</param>
162: /// <param name="WORK">
163: /// (workspace) DOUBLE PRECISION array, dimension (3*N)
164: ///</param>
165: /// <param name="IWORK">
166: /// (workspace) INTEGER array, dimension (N)
167: ///</param>
168: /// <param name="INFO">
169: /// (output) INTEGER
170: /// = 0: successful exit
171: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
172: ///</param>
173: public void Run(string NORM, string UPLO, string DIAG, int N, double[] A, int offset_a, int LDA
174: , ref double RCOND, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)
175: {
176:
177: #region Array Index Correction
178:
179: int o_a = -1 - LDA + offset_a; int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork;
180:
181: #endregion
182:
183:
184: #region Strings
185:
186: NORM = NORM.Substring(0, 1); UPLO = UPLO.Substring(0, 1); DIAG = DIAG.Substring(0, 1);
187:
188: #endregion
189:
190:
191: #region Prolog
192:
193: // *
194: // * -- LAPACK routine (version 3.0) --
195: // * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
196: // * Courant Institute, Argonne National Lab, and Rice University
197: // * March 31, 1993
198: // *
199: // * .. Scalar Arguments ..
200: // * ..
201: // * .. Array Arguments ..
202: // * ..
203: // *
204: // * Purpose
205: // * =======
206: // *
207: // * DTRCON estimates the reciprocal of the condition number of a
208: // * triangular matrix A, in either the 1-norm or the infinity-norm.
209: // *
210: // * The norm of A is computed and an estimate is obtained for
211: // * norm(inv(A)), then the reciprocal of the condition number is
212: // * computed as
213: // * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
214: // *
215: // * Arguments
216: // * =========
217: // *
218: // * NORM (input) CHARACTER*1
219: // * Specifies whether the 1-norm condition number or the
220: // * infinity-norm condition number is required:
221: // * = '1' or 'O': 1-norm;
222: // * = 'I': Infinity-norm.
223: // *
224: // * UPLO (input) CHARACTER*1
225: // * = 'U': A is upper triangular;
226: // * = 'L': A is lower triangular.
227: // *
228: // * DIAG (input) CHARACTER*1
229: // * = 'N': A is non-unit triangular;
230: // * = 'U': A is unit triangular.
231: // *
232: // * N (input) INTEGER
233: // * The order of the matrix A. N >= 0.
234: // *
235: // * A (input) DOUBLE PRECISION array, dimension (LDA,N)
236: // * The triangular matrix A. If UPLO = 'U', the leading N-by-N
237: // * upper triangular part of the array A contains the upper
238: // * triangular matrix, and the strictly lower triangular part of
239: // * A is not referenced. If UPLO = 'L', the leading N-by-N lower
240: // * triangular part of the array A contains the lower triangular
241: // * matrix, and the strictly upper triangular part of A is not
242: // * referenced. If DIAG = 'U', the diagonal elements of A are
243: // * also not referenced and are assumed to be 1.
244: // *
245: // * LDA (input) INTEGER
246: // * The leading dimension of the array A. LDA >= max(1,N).
247: // *
248: // * RCOND (output) DOUBLE PRECISION
249: // * The reciprocal of the condition number of the matrix A,
250: // * computed as RCOND = 1/(norm(A) * norm(inv(A))).
251: // *
252: // * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
253: // *
254: // * IWORK (workspace) INTEGER array, dimension (N)
255: // *
256: // * INFO (output) INTEGER
257: // * = 0: successful exit
258: // * < 0: if INFO = -i, the i-th argument had an illegal value
259: // *
260: // * =====================================================================
261: // *
262: // * .. Parameters ..
263: // * ..
264: // * .. Local Scalars ..
265: // * ..
266: // * .. External Functions ..
267: // * ..
268: // * .. External Subroutines ..
269: // * ..
270: // * .. Intrinsic Functions ..
271: // INTRINSIC ABS, DBLE, MAX;
272: // * ..
273: // * .. Executable Statements ..
274: // *
275: // * Test the input parameters.
276: // *
277:
278: #endregion
279:
280:
281: #region Body
282:
283: INFO = 0;
284: UPPER = this._lsame.Run(UPLO, "U");
285: ONENRM = NORM == "1" || this._lsame.Run(NORM, "O");
286: NOUNIT = this._lsame.Run(DIAG, "N");
287: // *
288: if (!ONENRM && !this._lsame.Run(NORM, "I"))
289: {
290: INFO = - 1;
291: }
292: else
293: {
294: if (!UPPER && !this._lsame.Run(UPLO, "L"))
295: {
296: INFO = - 2;
297: }
298: else
299: {
300: if (!NOUNIT && !this._lsame.Run(DIAG, "U"))
301: {
302: INFO = - 3;
303: }
304: else
305: {
306: if (N < 0)
307: {
308: INFO = - 4;
309: }
310: else
311: {
312: if (LDA < Math.Max(1, N))
313: {
314: INFO = - 6;
315: }
316: }
317: }
318: }
319: }
320: if (INFO != 0)
321: {
322: this._xerbla.Run("DTRCON", - INFO);
323: return;
324: }
325: // *
326: // * Quick return if possible
327: // *
328: if (N == 0)
329: {
330: RCOND = ONE;
331: return;
332: }
333: // *
334: RCOND = ZERO;
335: SMLNUM = this._dlamch.Run("Safe minimum") * Convert.ToDouble(Math.Max(1, N));
336: // *
337: // * Compute the norm of the triangular matrix A.
338: // *
339: ANORM = this._dlantr.Run(NORM, UPLO, DIAG, N, N, A, offset_a, LDA, ref WORK, offset_work);
340: // *
341: // * Continue only if ANORM > 0.
342: // *
343: if (ANORM > ZERO)
344: {
345: // *
346: // * Estimate the norm of the inverse of A.
347: // *
348: AINVNM = ZERO;
349: FortranLib.Copy(ref NORMIN , "N");
350: if (ONENRM)
351: {
352: KASE1 = 1;
353: }
354: else
355: {
356: KASE1 = 2;
357: }
358: KASE = 0;
359: LABEL10:;
360: this._dlacon.Run(N, ref WORK, N + 1 + o_work, ref WORK, offset_work, ref IWORK, offset_iwork, ref AINVNM, ref KASE);
361: if (KASE != 0)
362: {
363: if (KASE == KASE1)
364: {
365: // *
366: // * Multiply by inv(A).
367: // *
368: this._dlatrs.Run(UPLO, "No transpose", DIAG, NORMIN, N, A, offset_a
369: , LDA, ref WORK, offset_work, ref SCALE, ref WORK, 2 * N + 1 + o_work, ref INFO);
370: }
371: else
372: {
373: // *
374: // * Multiply by inv(A').
375: // *
376: this._dlatrs.Run(UPLO, "Transpose", DIAG, NORMIN, N, A, offset_a
377: , LDA, ref WORK, offset_work, ref SCALE, ref WORK, 2 * N + 1 + o_work, ref INFO);
378: }
379: FortranLib.Copy(ref NORMIN , "Y");
380: // *
381: // * Multiply by 1/SCALE if doing so will not cause overflow.
382: // *
383: if (SCALE != ONE)
384: {
385: IX = this._idamax.Run(N, WORK, offset_work, 1);
386: XNORM = Math.Abs(WORK[IX + o_work]);
387: if (SCALE < XNORM * SMLNUM || SCALE == ZERO) goto LABEL20;
388: this._drscl.Run(N, SCALE, ref WORK, offset_work, 1);
389: }
390: goto LABEL10;
391: }
392: // *
393: // * Compute the estimate of the reciprocal condition number.
394: // *
395: if (AINVNM != ZERO) RCOND = (ONE / ANORM) / AINVNM;
396: }
397: // *
398: LABEL20:;
399: return;
400: // *
401: // * End of DTRCON
402: // *
403:
404: #endregion
405:
406: }
407: }
408: }