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: /// DGETRI computes the inverse of a matrix using the LU factorization
27: /// computed by DGETRF.
28: ///
29: /// This method inverts U and then computes inv(A) by solving the system
30: /// inv(A)*L = inv(U) for inv(A).
31: ///
32: ///</summary>
33: public class DGETRI
34: {
35:
36:
37: #region Dependencies
38:
39: ILAENV _ilaenv; DGEMM _dgemm; DGEMV _dgemv; DSWAP _dswap; DTRSM _dtrsm; DTRTRI _dtrtri; XERBLA _xerbla;
40:
41: #endregion
42:
43:
44: #region Fields
45:
46: const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; int I = 0; int IWS = 0; int J = 0; int JB = 0;
47: int JJ = 0;int JP = 0; int LDWORK = 0; int LWKOPT = 0; int NB = 0; int NBMIN = 0; int NN = 0;
48:
49: #endregion
50:
51: public DGETRI(ILAENV ilaenv, DGEMM dgemm, DGEMV dgemv, DSWAP dswap, DTRSM dtrsm, DTRTRI dtrtri, XERBLA xerbla)
52: {
53:
54:
55: #region Set Dependencies
56:
57: this._ilaenv = ilaenv; this._dgemm = dgemm; this._dgemv = dgemv; this._dswap = dswap; this._dtrsm = dtrsm;
58: this._dtrtri = dtrtri;this._xerbla = xerbla;
59:
60: #endregion
61:
62: }
63:
64: public DGETRI()
65: {
66:
67:
68: #region Dependencies (Initialization)
69:
70: IEEECK ieeeck = new IEEECK();
71: IPARMQ iparmq = new IPARMQ();
72: LSAME lsame = new LSAME();
73: XERBLA xerbla = new XERBLA();
74: DSWAP dswap = new DSWAP();
75: DSCAL dscal = new DSCAL();
76: ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
77: DGEMM dgemm = new DGEMM(lsame, xerbla);
78: DGEMV dgemv = new DGEMV(lsame, xerbla);
79: DTRSM dtrsm = new DTRSM(lsame, xerbla);
80: DTRMM dtrmm = new DTRMM(lsame, xerbla);
81: DTRMV dtrmv = new DTRMV(lsame, xerbla);
82: DTRTI2 dtrti2 = new DTRTI2(lsame, dscal, dtrmv, xerbla);
83: DTRTRI dtrtri = new DTRTRI(lsame, ilaenv, dtrmm, dtrsm, dtrti2, xerbla);
84:
85: #endregion
86:
87:
88: #region Set Dependencies
89:
90: this._ilaenv = ilaenv; this._dgemm = dgemm; this._dgemv = dgemv; this._dswap = dswap; this._dtrsm = dtrsm;
91: this._dtrtri = dtrtri;this._xerbla = xerbla;
92:
93: #endregion
94:
95: }
96: /// <summary>
97: /// Purpose
98: /// =======
99: ///
100: /// DGETRI computes the inverse of a matrix using the LU factorization
101: /// computed by DGETRF.
102: ///
103: /// This method inverts U and then computes inv(A) by solving the system
104: /// inv(A)*L = inv(U) for inv(A).
105: ///
106: ///</summary>
107: /// <param name="N">
108: /// (input) INTEGER
109: /// The order of the matrix A. N .GE. 0.
110: ///</param>
111: /// <param name="A">
112: /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
113: /// On entry, the factors L and U from the factorization
114: /// A = P*L*U as computed by DGETRF.
115: /// On exit, if INFO = 0, the inverse of the original matrix A.
116: ///</param>
117: /// <param name="LDA">
118: /// (input) INTEGER
119: /// The leading dimension of the array A. LDA .GE. max(1,N).
120: ///</param>
121: /// <param name="IPIV">
122: /// (input) INTEGER array, dimension (N)
123: /// The pivot indices from DGETRF; for 1.LE.i.LE.N, row i of the
124: /// matrix was interchanged with row IPIV(i).
125: ///</param>
126: /// <param name="WORK">
127: /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
128: /// On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
129: ///</param>
130: /// <param name="LWORK">
131: /// (input) INTEGER
132: /// The dimension of the array WORK. LWORK .GE. max(1,N).
133: /// For optimal performance LWORK .GE. N*NB, where NB is
134: /// the optimal blocksize returned by ILAENV.
135: ///
136: /// If LWORK = -1, then a workspace query is assumed; the routine
137: /// only calculates the optimal size of the WORK array, returns
138: /// this value as the first entry of the WORK array, and no error
139: /// message related to LWORK is issued by XERBLA.
140: ///</param>
141: /// <param name="INFO">
142: /// (output) INTEGER
143: /// = 0: successful exit
144: /// .LT. 0: if INFO = -i, the i-th argument had an illegal value
145: /// .GT. 0: if INFO = i, U(i,i) is exactly zero; the matrix is
146: /// singular and its inverse could not be computed.
147: ///</param>
148: public void Run(int N, ref double[] A, int offset_a, int LDA, int[] IPIV, int offset_ipiv, ref double[] WORK, int offset_work, int LWORK
149: , ref int INFO)
150: {
151:
152: #region Array Index Correction
153:
154: int o_a = -1 - LDA + offset_a; int o_ipiv = -1 + offset_ipiv; int o_work = -1 + offset_work;
155:
156: #endregion
157:
158:
159: #region Prolog
160:
161: // *
162: // * -- LAPACK routine (version 3.1) --
163: // * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
164: // * November 2006
165: // *
166: // * .. Scalar Arguments ..
167: // * ..
168: // * .. Array Arguments ..
169: // * ..
170: // *
171: // * Purpose
172: // * =======
173: // *
174: // * DGETRI computes the inverse of a matrix using the LU factorization
175: // * computed by DGETRF.
176: // *
177: // * This method inverts U and then computes inv(A) by solving the system
178: // * inv(A)*L = inv(U) for inv(A).
179: // *
180: // * Arguments
181: // * =========
182: // *
183: // * N (input) INTEGER
184: // * The order of the matrix A. N >= 0.
185: // *
186: // * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
187: // * On entry, the factors L and U from the factorization
188: // * A = P*L*U as computed by DGETRF.
189: // * On exit, if INFO = 0, the inverse of the original matrix A.
190: // *
191: // * LDA (input) INTEGER
192: // * The leading dimension of the array A. LDA >= max(1,N).
193: // *
194: // * IPIV (input) INTEGER array, dimension (N)
195: // * The pivot indices from DGETRF; for 1<=i<=N, row i of the
196: // * matrix was interchanged with row IPIV(i).
197: // *
198: // * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
199: // * On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
200: // *
201: // * LWORK (input) INTEGER
202: // * The dimension of the array WORK. LWORK >= max(1,N).
203: // * For optimal performance LWORK >= N*NB, where NB is
204: // * the optimal blocksize returned by ILAENV.
205: // *
206: // * If LWORK = -1, then a workspace query is assumed; the routine
207: // * only calculates the optimal size of the WORK array, returns
208: // * this value as the first entry of the WORK array, and no error
209: // * message related to LWORK is issued by XERBLA.
210: // *
211: // * INFO (output) INTEGER
212: // * = 0: successful exit
213: // * < 0: if INFO = -i, the i-th argument had an illegal value
214: // * > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
215: // * singular and its inverse could not be computed.
216: // *
217: // * =====================================================================
218: // *
219: // * .. Parameters ..
220: // * ..
221: // * .. Local Scalars ..
222: // * ..
223: // * .. External Functions ..
224: // * ..
225: // * .. External Subroutines ..
226: // * ..
227: // * .. Intrinsic Functions ..
228: // INTRINSIC MAX, MIN;
229: // * ..
230: // * .. Executable Statements ..
231: // *
232: // * Test the input parameters.
233: // *
234:
235: #endregion
236:
237:
238: #region Body
239:
240: INFO = 0;
241: NB = this._ilaenv.Run(1, "DGETRI", " ", N, - 1, - 1, - 1);
242: LWKOPT = N * NB;
243: WORK[1 + o_work] = LWKOPT;
244: LQUERY = (LWORK == - 1);
245: if (N < 0)
246: {
247: INFO = - 1;
248: }
249: else
250: {
251: if (LDA < Math.Max(1, N))
252: {
253: INFO = - 3;
254: }
255: else
256: {
257: if (LWORK < Math.Max(1, N) && !LQUERY)
258: {
259: INFO = - 6;
260: }
261: }
262: }
263: if (INFO != 0)
264: {
265: this._xerbla.Run("DGETRI", - INFO);
266: return;
267: }
268: else
269: {
270: if (LQUERY)
271: {
272: return;
273: }
274: }
275: // *
276: // * Quick return if possible
277: // *
278: if (N == 0) return;
279: // *
280: // * Form inv(U). If INFO > 0 from DTRTRI, then U is singular,
281: // * and the inverse is not computed.
282: // *
283: this._dtrtri.Run("Upper", "Non-unit", N, ref A, offset_a, LDA, ref INFO);
284: if (INFO > 0) return;
285: // *
286: NBMIN = 2;
287: LDWORK = N;
288: if (NB > 1 && NB < N)
289: {
290: IWS = Math.Max(LDWORK * NB, 1);
291: if (LWORK < IWS)
292: {
293: NB = LWORK / LDWORK;
294: NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGETRI", " ", N, - 1, - 1, - 1));
295: }
296: }
297: else
298: {
299: IWS = N;
300: }
301: // *
302: // * Solve the equation inv(A)*L = inv(U) for inv(A).
303: // *
304: if (NB < NBMIN || NB >= N)
305: {
306: // *
307: // * Use unblocked code.
308: // *
309: for (J = N; J >= 1; J += - 1)
310: {
311: // *
312: // * Copy current column of L to WORK and replace with zeros.
313: // *
314: for (I = J + 1; I <= N; I++)
315: {
316: WORK[I + o_work] = A[I+J * LDA + o_a];
317: A[I+J * LDA + o_a] = ZERO;
318: }
319: // *
320: // * Compute current column of inv(A).
321: // *
322: if (J < N)
323: {
324: this._dgemv.Run("No transpose", N, N - J, - ONE, A, 1+(J + 1) * LDA + o_a, LDA
325: , WORK, J + 1 + o_work, 1, ONE, ref A, 1+J * LDA + o_a, 1);
326: }
327: }
328: }
329: else
330: {
331: // *
332: // * Use blocked code.
333: // *
334: NN = ((N - 1) / NB) * NB + 1;
335: for (J = NN; ( - NB >= 0) ? (J <= 1) : (J >= 1); J += - NB)
336: {
337: JB = Math.Min(NB, N - J + 1);
338: // *
339: // * Copy current block column of L to WORK and replace with
340: // * zeros.
341: // *
342: for (JJ = J; JJ <= J + JB - 1; JJ++)
343: {
344: for (I = JJ + 1; I <= N; I++)
345: {
346: WORK[I + (JJ - J) * LDWORK + o_work] = A[I+JJ * LDA + o_a];
347: A[I+JJ * LDA + o_a] = ZERO;
348: }
349: }
350: // *
351: // * Compute current block column of inv(A).
352: // *
353: if (J + JB <= N)
354: {
355: this._dgemm.Run("No transpose", "No transpose", N, JB, N - J - JB + 1, - ONE
356: , A, 1+(J + JB) * LDA + o_a, LDA, WORK, J + JB + o_work, LDWORK, ONE, ref A, 1+J * LDA + o_a
357: , LDA);
358: }
359: this._dtrsm.Run("Right", "Lower", "No transpose", "Unit", N, JB
360: , ONE, WORK, J + o_work, LDWORK, ref A, 1+J * LDA + o_a, LDA);
361: }
362: }
363: // *
364: // * Apply column interchanges.
365: // *
366: for (J = N - 1; J >= 1; J += - 1)
367: {
368: JP = IPIV[J + o_ipiv];
369: if (JP != J) this._dswap.Run(N, ref A, 1+J * LDA + o_a, 1, ref A, 1+JP * LDA + o_a, 1);
370: }
371: // *
372: WORK[1 + o_work] = IWS;
373: return;
374: // *
375: // * End of DGETRI
376: // *
377:
378: #endregion
379:
380: }
381: }
382: }