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();
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:  }