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: /// DTRTRI computes the inverse of a real upper or lower triangular
27: /// matrix A.
28: ///
29: /// This is the Level 3 BLAS version of the algorithm.
30: ///
31: ///</summary>
32: public class DTRTRI
33: {
34:
35:
36: #region Dependencies
37:
38: LSAME _lsame; ILAENV _ilaenv; DTRMM _dtrmm; DTRSM _dtrsm; DTRTI2 _dtrti2; XERBLA _xerbla;
39:
40: #endregion
41:
42:
43: #region Fields
44:
45: const double ONE = 1.0E+0; const double ZERO = 0.0E+0; bool NOUNIT = false; bool UPPER = false; int J = 0; int JB = 0;
46: int NB = 0;int NN = 0;
47:
48: #endregion
49:
50: public DTRTRI(LSAME lsame, ILAENV ilaenv, DTRMM dtrmm, DTRSM dtrsm, DTRTI2 dtrti2, XERBLA xerbla)
51: {
52:
53:
54: #region Set Dependencies
55:
56: this._lsame = lsame; this._ilaenv = ilaenv; this._dtrmm = dtrmm; this._dtrsm = dtrsm; this._dtrti2 = dtrti2;
57: this._xerbla = xerbla;
58:
59: #endregion
60:
61: }
62:
63: public DTRTRI()
64: {
65:
66:
67: #region Dependencies (Initialization)
68:
69: LSAME lsame = new LSAME();
70: IEEECK ieeeck = new IEEECK();
71: IPARMQ iparmq = new IPARMQ();
72: XERBLA xerbla = new XERBLA();
73: DSCAL dscal = new DSCAL();
74: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
75: DTRMM dtrmm = new DTRMM(lsame, xerbla);
76: DTRSM dtrsm = new DTRSM(lsame, xerbla);
77: DTRMV dtrmv = new DTRMV(lsame, xerbla);
78: DTRTI2 dtrti2 = new DTRTI2(lsame, dscal, dtrmv, xerbla);
79:
80: #endregion
81:
82:
83: #region Set Dependencies
84:
85: this._lsame = lsame; this._ilaenv = ilaenv; this._dtrmm = dtrmm; this._dtrsm = dtrsm; this._dtrti2 = dtrti2;
86: this._xerbla = xerbla;
87:
88: #endregion
89:
90: }
91: /// <summary>
92: /// Purpose
93: /// =======
94: ///
95: /// DTRTRI computes the inverse of a real upper or lower triangular
96: /// matrix A.
97: ///
98: /// This is the Level 3 BLAS version of the algorithm.
99: ///
100: ///</summary>
101: /// <param name="UPLO">
102: /// (input) CHARACTER*1
103: /// = 'U': A is upper triangular;
104: /// = 'L': A is lower triangular.
105: ///</param>
106: /// <param name="DIAG">
107: /// (input) CHARACTER*1
108: /// = 'N': A is non-unit triangular;
109: /// = 'U': A is unit triangular.
110: ///</param>
111: /// <param name="N">
112: /// (input) INTEGER
113: /// The order of the matrix A. N .GE. 0.
114: ///</param>
115: /// <param name="A">
116: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
117: /// On entry, the triangular matrix A. If UPLO = 'U', the
118: /// leading N-by-N upper triangular part of the array A contains
119: /// the upper triangular matrix, and the strictly lower
120: /// triangular part of A is not referenced. If UPLO = 'L', the
121: /// leading N-by-N lower triangular part of the array A contains
122: /// the lower triangular matrix, and the strictly upper
123: /// triangular part of A is not referenced. If DIAG = 'U', the
124: /// diagonal elements of A are also not referenced and are
125: /// assumed to be 1.
126: /// On exit, the (triangular) inverse of the original matrix, in
127: /// the same storage format.
128: ///</param>
129: /// <param name="LDA">
130: /// (input) INTEGER
131: /// The leading dimension of the array A. LDA .GE. max(1,N).
132: ///</param>
133: /// <param name="INFO">
134: /// (output) INTEGER
135: /// = 0: successful exit
136: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
137: /// .GT. 0: if INFO = i, A(i,i) is exactly zero. The triangular
138: /// matrix is singular and its inverse can not be computed.
139: ///</param>
140: public void Run(string UPLO, string DIAG, int N, ref double[] A, int offset_a, int LDA, ref int INFO)
141: {
142:
143: #region Array Index Correction
144:
145: int o_a = -1 - LDA + offset_a;
146:
147: #endregion
148:
149:
150: #region Strings
151:
152: UPLO = UPLO.Substring(0, 1); DIAG = DIAG.Substring(0, 1);
153:
154: #endregion
155:
156:
157: #region Prolog
158:
159: // *
160: // * -- LAPACK routine (version 3.1) --
161: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
162: // * November 2006
163: // *
164: // * .. Scalar Arguments ..
165: // * ..
166: // * .. Array Arguments ..
167: // * ..
168: // *
169: // * Purpose
170: // * =======
171: // *
172: // * DTRTRI computes the inverse of a real upper or lower triangular
173: // * matrix A.
174: // *
175: // * This is the Level 3 BLAS version of the algorithm.
176: // *
177: // * Arguments
178: // * =========
179: // *
180: // * UPLO (input) CHARACTER*1
181: // * = 'U': A is upper triangular;
182: // * = 'L': A is lower triangular.
183: // *
184: // * DIAG (input) CHARACTER*1
185: // * = 'N': A is non-unit triangular;
186: // * = 'U': A is unit triangular.
187: // *
188: // * N (input) INTEGER
189: // * The order of the matrix A. N >= 0.
190: // *
191: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
192: // * On entry, the triangular matrix A. If UPLO = 'U', the
193: // * leading N-by-N upper triangular part of the array A contains
194: // * the upper triangular matrix, and the strictly lower
195: // * triangular part of A is not referenced. If UPLO = 'L', the
196: // * leading N-by-N lower triangular part of the array A contains
197: // * the lower triangular matrix, and the strictly upper
198: // * triangular part of A is not referenced. If DIAG = 'U', the
199: // * diagonal elements of A are also not referenced and are
200: // * assumed to be 1.
201: // * On exit, the (triangular) inverse of the original matrix, in
202: // * the same storage format.
203: // *
204: // * LDA (input) INTEGER
205: // * The leading dimension of the array A. LDA >= max(1,N).
206: // *
207: // * INFO (output) INTEGER
208: // * = 0: successful exit
209: // * < 0: if INFO = -i, the i-th argument had an illegal value
210: // * > 0: if INFO = i, A(i,i) is exactly zero. The triangular
211: // * matrix is singular and its inverse can not be computed.
212: // *
213: // * =====================================================================
214: // *
215: // * .. Parameters ..
216: // * ..
217: // * .. Local Scalars ..
218: // * ..
219: // * .. External Functions ..
220: // * ..
221: // * .. External Subroutines ..
222: // * ..
223: // * .. Intrinsic Functions ..
224: // INTRINSIC MAX, MIN;
225: // * ..
226: // * .. Executable Statements ..
227: // *
228: // * Test the input parameters.
229: // *
230:
231: #endregion
232:
233:
234: #region Body
235:
236: INFO = 0;
237: UPPER = this._lsame.Run(UPLO, "U");
238: NOUNIT = this._lsame.Run(DIAG, "N");
239: if (!UPPER && !this._lsame.Run(UPLO, "L"))
240: {
241: INFO = - 1;
242: }
243: else
244: {
245: if (!NOUNIT && !this._lsame.Run(DIAG, "U"))
246: {
247: INFO = - 2;
248: }
249: else
250: {
251: if (N < 0)
252: {
253: INFO = - 3;
254: }
255: else
256: {
257: if (LDA < Math.Max(1, N))
258: {
259: INFO = - 5;
260: }
261: }
262: }
263: }
264: if (INFO != 0)
265: {
266: this._xerbla.Run("DTRTRI", - INFO);
267: return;
268: }
269: // *
270: // * Quick return if possible
271: // *
272: if (N == 0) return;
273: // *
274: // * Check for singularity if non-unit.
275: // *
276: if (NOUNIT)
277: {
278: for (INFO = 1; INFO <= N; INFO++)
279: {
280: if (A[INFO+INFO * LDA + o_a] == ZERO) return;
281: }
282: INFO = 0;
283: }
284: // *
285: // * Determine the block size for this environment.
286: // *
287: NB = this._ilaenv.Run(1, "DTRTRI", UPLO + DIAG, N, - 1, - 1, - 1);
288: if (NB <= 1 || NB >= N)
289: {
290: // *
291: // * Use unblocked code
292: // *
293: this._dtrti2.Run(UPLO, DIAG, N, ref A, offset_a, LDA, ref INFO);
294: }
295: else
296: {
297: // *
298: // * Use blocked code
299: // *
300: if (UPPER)
301: {
302: // *
303: // * Compute inverse of upper triangular matrix
304: // *
305: for (J = 1; (NB >= 0) ? (J <= N) : (J >= N); J += NB)
306: {
307: JB = Math.Min(NB, N - J + 1);
308: // *
309: // * Compute rows 1:j-1 of current block column
310: // *
311: this._dtrmm.Run("Left", "Upper", "No transpose", DIAG, J - 1, JB
312: , ONE, A, offset_a, LDA, ref A, 1+J * LDA + o_a, LDA);
313: this._dtrsm.Run("Right", "Upper", "No transpose", DIAG, J - 1, JB
314: , - ONE, A, J+J * LDA + o_a, LDA, ref A, 1+J * LDA + o_a, LDA);
315: // *
316: // * Compute inverse of current diagonal block
317: // *
318: this._dtrti2.Run("Upper", DIAG, JB, ref A, J+J * LDA + o_a, LDA, ref INFO);
319: }
320: }
321: else
322: {
323: // *
324: // * Compute inverse of lower triangular matrix
325: // *
326: NN = ((N - 1) / NB) * NB + 1;
327: for (J = NN; ( - NB >= 0) ? (J <= 1) : (J >= 1); J += - NB)
328: {
329: JB = Math.Min(NB, N - J + 1);
330: if (J + JB <= N)
331: {
332: // *
333: // * Compute rows j+jb:n of current block column
334: // *
335: this._dtrmm.Run("Left", "Lower", "No transpose", DIAG, N - J - JB + 1, JB
336: , ONE, A, J + JB+(J + JB) * LDA + o_a, LDA, ref A, J + JB+J * LDA + o_a, LDA);
337: this._dtrsm.Run("Right", "Lower", "No transpose", DIAG, N - J - JB + 1, JB
338: , - ONE, A, J+J * LDA + o_a, LDA, ref A, J + JB+J * LDA + o_a, LDA);
339: }
340: // *
341: // * Compute inverse of current diagonal block
342: // *
343: this._dtrti2.Run("Lower", DIAG, JB, ref A, J+J * LDA + o_a, LDA, ref INFO);
344: }
345: }
346: }
347: // *
348: return;
349: // *
350: // * End of DTRTRI
351: // *
352:
353: #endregion
354:
355: }
356: }
357: }