![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGETRI 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGETRI + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgetri.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgetri.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgetri.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, LWORK, N 00025 * .. 00026 * .. Array Arguments .. 00027 * INTEGER IPIV( * ) 00028 * COMPLEX A( LDA, * ), WORK( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> CGETRI computes the inverse of a matrix using the LU factorization 00038 *> computed by CGETRF. 00039 *> 00040 *> This method inverts U and then computes inv(A) by solving the system 00041 *> inv(A)*L = inv(U) for inv(A). 00042 *> \endverbatim 00043 * 00044 * Arguments: 00045 * ========== 00046 * 00047 *> \param[in] N 00048 *> \verbatim 00049 *> N is INTEGER 00050 *> The order of the matrix A. N >= 0. 00051 *> \endverbatim 00052 *> 00053 *> \param[in,out] A 00054 *> \verbatim 00055 *> A is COMPLEX array, dimension (LDA,N) 00056 *> On entry, the factors L and U from the factorization 00057 *> A = P*L*U as computed by CGETRF. 00058 *> On exit, if INFO = 0, the inverse of the original matrix A. 00059 *> \endverbatim 00060 *> 00061 *> \param[in] LDA 00062 *> \verbatim 00063 *> LDA is INTEGER 00064 *> The leading dimension of the array A. LDA >= max(1,N). 00065 *> \endverbatim 00066 *> 00067 *> \param[in] IPIV 00068 *> \verbatim 00069 *> IPIV is INTEGER array, dimension (N) 00070 *> The pivot indices from CGETRF; for 1<=i<=N, row i of the 00071 *> matrix was interchanged with row IPIV(i). 00072 *> \endverbatim 00073 *> 00074 *> \param[out] WORK 00075 *> \verbatim 00076 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00077 *> On exit, if INFO=0, then WORK(1) returns the optimal LWORK. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] LWORK 00081 *> \verbatim 00082 *> LWORK is INTEGER 00083 *> The dimension of the array WORK. LWORK >= max(1,N). 00084 *> For optimal performance LWORK >= N*NB, where NB is 00085 *> the optimal blocksize returned by ILAENV. 00086 *> 00087 *> If LWORK = -1, then a workspace query is assumed; the routine 00088 *> only calculates the optimal size of the WORK array, returns 00089 *> this value as the first entry of the WORK array, and no error 00090 *> message related to LWORK is issued by XERBLA. 00091 *> \endverbatim 00092 *> 00093 *> \param[out] INFO 00094 *> \verbatim 00095 *> INFO is INTEGER 00096 *> = 0: successful exit 00097 *> < 0: if INFO = -i, the i-th argument had an illegal value 00098 *> > 0: if INFO = i, U(i,i) is exactly zero; the matrix is 00099 *> singular and its inverse could not be computed. 00100 *> \endverbatim 00101 * 00102 * Authors: 00103 * ======== 00104 * 00105 *> \author Univ. of Tennessee 00106 *> \author Univ. of California Berkeley 00107 *> \author Univ. of Colorado Denver 00108 *> \author NAG Ltd. 00109 * 00110 *> \date November 2011 00111 * 00112 *> \ingroup complexGEcomputational 00113 * 00114 * ===================================================================== 00115 SUBROUTINE CGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) 00116 * 00117 * -- LAPACK computational routine (version 3.4.0) -- 00118 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00120 * November 2011 00121 * 00122 * .. Scalar Arguments .. 00123 INTEGER INFO, LDA, LWORK, N 00124 * .. 00125 * .. Array Arguments .. 00126 INTEGER IPIV( * ) 00127 COMPLEX A( LDA, * ), WORK( * ) 00128 * .. 00129 * 00130 * ===================================================================== 00131 * 00132 * .. Parameters .. 00133 COMPLEX ZERO, ONE 00134 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), 00135 $ ONE = ( 1.0E+0, 0.0E+0 ) ) 00136 * .. 00137 * .. Local Scalars .. 00138 LOGICAL LQUERY 00139 INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB, 00140 $ NBMIN, NN 00141 * .. 00142 * .. External Functions .. 00143 INTEGER ILAENV 00144 EXTERNAL ILAENV 00145 * .. 00146 * .. External Subroutines .. 00147 EXTERNAL CGEMM, CGEMV, CSWAP, CTRSM, CTRTRI, XERBLA 00148 * .. 00149 * .. Intrinsic Functions .. 00150 INTRINSIC MAX, MIN 00151 * .. 00152 * .. Executable Statements .. 00153 * 00154 * Test the input parameters. 00155 * 00156 INFO = 0 00157 NB = ILAENV( 1, 'CGETRI', ' ', N, -1, -1, -1 ) 00158 LWKOPT = N*NB 00159 WORK( 1 ) = LWKOPT 00160 LQUERY = ( LWORK.EQ.-1 ) 00161 IF( N.LT.0 ) THEN 00162 INFO = -1 00163 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00164 INFO = -3 00165 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00166 INFO = -6 00167 END IF 00168 IF( INFO.NE.0 ) THEN 00169 CALL XERBLA( 'CGETRI', -INFO ) 00170 RETURN 00171 ELSE IF( LQUERY ) THEN 00172 RETURN 00173 END IF 00174 * 00175 * Quick return if possible 00176 * 00177 IF( N.EQ.0 ) 00178 $ RETURN 00179 * 00180 * Form inv(U). If INFO > 0 from CTRTRI, then U is singular, 00181 * and the inverse is not computed. 00182 * 00183 CALL CTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO ) 00184 IF( INFO.GT.0 ) 00185 $ RETURN 00186 * 00187 NBMIN = 2 00188 LDWORK = N 00189 IF( NB.GT.1 .AND. NB.LT.N ) THEN 00190 IWS = MAX( LDWORK*NB, 1 ) 00191 IF( LWORK.LT.IWS ) THEN 00192 NB = LWORK / LDWORK 00193 NBMIN = MAX( 2, ILAENV( 2, 'CGETRI', ' ', N, -1, -1, -1 ) ) 00194 END IF 00195 ELSE 00196 IWS = N 00197 END IF 00198 * 00199 * Solve the equation inv(A)*L = inv(U) for inv(A). 00200 * 00201 IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN 00202 * 00203 * Use unblocked code. 00204 * 00205 DO 20 J = N, 1, -1 00206 * 00207 * Copy current column of L to WORK and replace with zeros. 00208 * 00209 DO 10 I = J + 1, N 00210 WORK( I ) = A( I, J ) 00211 A( I, J ) = ZERO 00212 10 CONTINUE 00213 * 00214 * Compute current column of inv(A). 00215 * 00216 IF( J.LT.N ) 00217 $ CALL CGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ), 00218 $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) 00219 20 CONTINUE 00220 ELSE 00221 * 00222 * Use blocked code. 00223 * 00224 NN = ( ( N-1 ) / NB )*NB + 1 00225 DO 50 J = NN, 1, -NB 00226 JB = MIN( NB, N-J+1 ) 00227 * 00228 * Copy current block column of L to WORK and replace with 00229 * zeros. 00230 * 00231 DO 40 JJ = J, J + JB - 1 00232 DO 30 I = JJ + 1, N 00233 WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) 00234 A( I, JJ ) = ZERO 00235 30 CONTINUE 00236 40 CONTINUE 00237 * 00238 * Compute current block column of inv(A). 00239 * 00240 IF( J+JB.LE.N ) 00241 $ CALL CGEMM( 'No transpose', 'No transpose', N, JB, 00242 $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA, 00243 $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA ) 00244 CALL CTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB, 00245 $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) 00246 50 CONTINUE 00247 END IF 00248 * 00249 * Apply column interchanges. 00250 * 00251 DO 60 J = N - 1, 1, -1 00252 JP = IPIV( J ) 00253 IF( JP.NE.J ) 00254 $ CALL CSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 ) 00255 60 CONTINUE 00256 * 00257 WORK( 1 ) = IWS 00258 RETURN 00259 * 00260 * End of CGETRI 00261 * 00262 END