1 SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, LWORK, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 COMPLEX*16 A( LDA, * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZGETRI computes the inverse of a matrix using the LU factorization
20 * computed by ZGETRF.
21 *
22 * This method inverts U and then computes inv(A) by solving the system
23 * inv(A)*L = inv(U) for inv(A).
24 *
25 * Arguments
26 * =========
27 *
28 * N (input) INTEGER
29 * The order of the matrix A. N >= 0.
30 *
31 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
32 * On entry, the factors L and U from the factorization
33 * A = P*L*U as computed by ZGETRF.
34 * On exit, if INFO = 0, the inverse of the original matrix A.
35 *
36 * LDA (input) INTEGER
37 * The leading dimension of the array A. LDA >= max(1,N).
38 *
39 * IPIV (input) INTEGER array, dimension (N)
40 * The pivot indices from ZGETRF; for 1<=i<=N, row i of the
41 * matrix was interchanged with row IPIV(i).
42 *
43 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
44 * On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
45 *
46 * LWORK (input) INTEGER
47 * The dimension of the array WORK. LWORK >= max(1,N).
48 * For optimal performance LWORK >= N*NB, where NB is
49 * the optimal blocksize returned by ILAENV.
50 *
51 * If LWORK = -1, then a workspace query is assumed; the routine
52 * only calculates the optimal size of the WORK array, returns
53 * this value as the first entry of the WORK array, and no error
54 * message related to LWORK is issued by XERBLA.
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 * > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
60 * singular and its inverse could not be computed.
61 *
62 * =====================================================================
63 *
64 * .. Parameters ..
65 COMPLEX*16 ZERO, ONE
66 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
67 $ ONE = ( 1.0D+0, 0.0D+0 ) )
68 * ..
69 * .. Local Scalars ..
70 LOGICAL LQUERY
71 INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
72 $ NBMIN, NN
73 * ..
74 * .. External Functions ..
75 INTEGER ILAENV
76 EXTERNAL ILAENV
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL XERBLA, ZGEMM, ZGEMV, ZSWAP, ZTRSM, ZTRTRI
80 * ..
81 * .. Intrinsic Functions ..
82 INTRINSIC MAX, MIN
83 * ..
84 * .. Executable Statements ..
85 *
86 * Test the input parameters.
87 *
88 INFO = 0
89 NB = ILAENV( 1, 'ZGETRI', ' ', N, -1, -1, -1 )
90 LWKOPT = N*NB
91 WORK( 1 ) = LWKOPT
92 LQUERY = ( LWORK.EQ.-1 )
93 IF( N.LT.0 ) THEN
94 INFO = -1
95 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
96 INFO = -3
97 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
98 INFO = -6
99 END IF
100 IF( INFO.NE.0 ) THEN
101 CALL XERBLA( 'ZGETRI', -INFO )
102 RETURN
103 ELSE IF( LQUERY ) THEN
104 RETURN
105 END IF
106 *
107 * Quick return if possible
108 *
109 IF( N.EQ.0 )
110 $ RETURN
111 *
112 * Form inv(U). If INFO > 0 from ZTRTRI, then U is singular,
113 * and the inverse is not computed.
114 *
115 CALL ZTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
116 IF( INFO.GT.0 )
117 $ RETURN
118 *
119 NBMIN = 2
120 LDWORK = N
121 IF( NB.GT.1 .AND. NB.LT.N ) THEN
122 IWS = MAX( LDWORK*NB, 1 )
123 IF( LWORK.LT.IWS ) THEN
124 NB = LWORK / LDWORK
125 NBMIN = MAX( 2, ILAENV( 2, 'ZGETRI', ' ', N, -1, -1, -1 ) )
126 END IF
127 ELSE
128 IWS = N
129 END IF
130 *
131 * Solve the equation inv(A)*L = inv(U) for inv(A).
132 *
133 IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
134 *
135 * Use unblocked code.
136 *
137 DO 20 J = N, 1, -1
138 *
139 * Copy current column of L to WORK and replace with zeros.
140 *
141 DO 10 I = J + 1, N
142 WORK( I ) = A( I, J )
143 A( I, J ) = ZERO
144 10 CONTINUE
145 *
146 * Compute current column of inv(A).
147 *
148 IF( J.LT.N )
149 $ CALL ZGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
150 $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
151 20 CONTINUE
152 ELSE
153 *
154 * Use blocked code.
155 *
156 NN = ( ( N-1 ) / NB )*NB + 1
157 DO 50 J = NN, 1, -NB
158 JB = MIN( NB, N-J+1 )
159 *
160 * Copy current block column of L to WORK and replace with
161 * zeros.
162 *
163 DO 40 JJ = J, J + JB - 1
164 DO 30 I = JJ + 1, N
165 WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
166 A( I, JJ ) = ZERO
167 30 CONTINUE
168 40 CONTINUE
169 *
170 * Compute current block column of inv(A).
171 *
172 IF( J+JB.LE.N )
173 $ CALL ZGEMM( 'No transpose', 'No transpose', N, JB,
174 $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
175 $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
176 CALL ZTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
177 $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
178 50 CONTINUE
179 END IF
180 *
181 * Apply column interchanges.
182 *
183 DO 60 J = N - 1, 1, -1
184 JP = IPIV( J )
185 IF( JP.NE.J )
186 $ CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
187 60 CONTINUE
188 *
189 WORK( 1 ) = IWS
190 RETURN
191 *
192 * End of ZGETRI
193 *
194 END
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, LWORK, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 COMPLEX*16 A( LDA, * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZGETRI computes the inverse of a matrix using the LU factorization
20 * computed by ZGETRF.
21 *
22 * This method inverts U and then computes inv(A) by solving the system
23 * inv(A)*L = inv(U) for inv(A).
24 *
25 * Arguments
26 * =========
27 *
28 * N (input) INTEGER
29 * The order of the matrix A. N >= 0.
30 *
31 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
32 * On entry, the factors L and U from the factorization
33 * A = P*L*U as computed by ZGETRF.
34 * On exit, if INFO = 0, the inverse of the original matrix A.
35 *
36 * LDA (input) INTEGER
37 * The leading dimension of the array A. LDA >= max(1,N).
38 *
39 * IPIV (input) INTEGER array, dimension (N)
40 * The pivot indices from ZGETRF; for 1<=i<=N, row i of the
41 * matrix was interchanged with row IPIV(i).
42 *
43 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
44 * On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
45 *
46 * LWORK (input) INTEGER
47 * The dimension of the array WORK. LWORK >= max(1,N).
48 * For optimal performance LWORK >= N*NB, where NB is
49 * the optimal blocksize returned by ILAENV.
50 *
51 * If LWORK = -1, then a workspace query is assumed; the routine
52 * only calculates the optimal size of the WORK array, returns
53 * this value as the first entry of the WORK array, and no error
54 * message related to LWORK is issued by XERBLA.
55 *
56 * INFO (output) INTEGER
57 * = 0: successful exit
58 * < 0: if INFO = -i, the i-th argument had an illegal value
59 * > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
60 * singular and its inverse could not be computed.
61 *
62 * =====================================================================
63 *
64 * .. Parameters ..
65 COMPLEX*16 ZERO, ONE
66 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
67 $ ONE = ( 1.0D+0, 0.0D+0 ) )
68 * ..
69 * .. Local Scalars ..
70 LOGICAL LQUERY
71 INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
72 $ NBMIN, NN
73 * ..
74 * .. External Functions ..
75 INTEGER ILAENV
76 EXTERNAL ILAENV
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL XERBLA, ZGEMM, ZGEMV, ZSWAP, ZTRSM, ZTRTRI
80 * ..
81 * .. Intrinsic Functions ..
82 INTRINSIC MAX, MIN
83 * ..
84 * .. Executable Statements ..
85 *
86 * Test the input parameters.
87 *
88 INFO = 0
89 NB = ILAENV( 1, 'ZGETRI', ' ', N, -1, -1, -1 )
90 LWKOPT = N*NB
91 WORK( 1 ) = LWKOPT
92 LQUERY = ( LWORK.EQ.-1 )
93 IF( N.LT.0 ) THEN
94 INFO = -1
95 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
96 INFO = -3
97 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
98 INFO = -6
99 END IF
100 IF( INFO.NE.0 ) THEN
101 CALL XERBLA( 'ZGETRI', -INFO )
102 RETURN
103 ELSE IF( LQUERY ) THEN
104 RETURN
105 END IF
106 *
107 * Quick return if possible
108 *
109 IF( N.EQ.0 )
110 $ RETURN
111 *
112 * Form inv(U). If INFO > 0 from ZTRTRI, then U is singular,
113 * and the inverse is not computed.
114 *
115 CALL ZTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
116 IF( INFO.GT.0 )
117 $ RETURN
118 *
119 NBMIN = 2
120 LDWORK = N
121 IF( NB.GT.1 .AND. NB.LT.N ) THEN
122 IWS = MAX( LDWORK*NB, 1 )
123 IF( LWORK.LT.IWS ) THEN
124 NB = LWORK / LDWORK
125 NBMIN = MAX( 2, ILAENV( 2, 'ZGETRI', ' ', N, -1, -1, -1 ) )
126 END IF
127 ELSE
128 IWS = N
129 END IF
130 *
131 * Solve the equation inv(A)*L = inv(U) for inv(A).
132 *
133 IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
134 *
135 * Use unblocked code.
136 *
137 DO 20 J = N, 1, -1
138 *
139 * Copy current column of L to WORK and replace with zeros.
140 *
141 DO 10 I = J + 1, N
142 WORK( I ) = A( I, J )
143 A( I, J ) = ZERO
144 10 CONTINUE
145 *
146 * Compute current column of inv(A).
147 *
148 IF( J.LT.N )
149 $ CALL ZGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
150 $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
151 20 CONTINUE
152 ELSE
153 *
154 * Use blocked code.
155 *
156 NN = ( ( N-1 ) / NB )*NB + 1
157 DO 50 J = NN, 1, -NB
158 JB = MIN( NB, N-J+1 )
159 *
160 * Copy current block column of L to WORK and replace with
161 * zeros.
162 *
163 DO 40 JJ = J, J + JB - 1
164 DO 30 I = JJ + 1, N
165 WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
166 A( I, JJ ) = ZERO
167 30 CONTINUE
168 40 CONTINUE
169 *
170 * Compute current block column of inv(A).
171 *
172 IF( J+JB.LE.N )
173 $ CALL ZGEMM( 'No transpose', 'No transpose', N, JB,
174 $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
175 $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
176 CALL ZTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
177 $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
178 50 CONTINUE
179 END IF
180 *
181 * Apply column interchanges.
182 *
183 DO 60 J = N - 1, 1, -1
184 JP = IPIV( J )
185 IF( JP.NE.J )
186 $ CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
187 60 CONTINUE
188 *
189 WORK( 1 ) = IWS
190 RETURN
191 *
192 * End of ZGETRI
193 *
194 END