![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGETRI 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGETRI + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgetri.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgetri.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgetri.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGETRI( 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 * REAL A( LDA, * ), WORK( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SGETRI computes the inverse of a matrix using the LU factorization 00038 *> computed by SGETRF. 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 REAL array, dimension (LDA,N) 00056 *> On entry, the factors L and U from the factorization 00057 *> A = P*L*U as computed by SGETRF. 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 SGETRF; 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 REAL 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 realGEcomputational 00113 * 00114 * ===================================================================== 00115 SUBROUTINE SGETRI( 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 REAL A( LDA, * ), WORK( * ) 00128 * .. 00129 * 00130 * ===================================================================== 00131 * 00132 * .. Parameters .. 00133 REAL ZERO, ONE 00134 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00135 * .. 00136 * .. Local Scalars .. 00137 LOGICAL LQUERY 00138 INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB, 00139 $ NBMIN, NN 00140 * .. 00141 * .. External Functions .. 00142 INTEGER ILAENV 00143 EXTERNAL ILAENV 00144 * .. 00145 * .. External Subroutines .. 00146 EXTERNAL SGEMM, SGEMV, SSWAP, STRSM, STRTRI, XERBLA 00147 * .. 00148 * .. Intrinsic Functions .. 00149 INTRINSIC MAX, MIN 00150 * .. 00151 * .. Executable Statements .. 00152 * 00153 * Test the input parameters. 00154 * 00155 INFO = 0 00156 NB = ILAENV( 1, 'SGETRI', ' ', N, -1, -1, -1 ) 00157 LWKOPT = N*NB 00158 WORK( 1 ) = LWKOPT 00159 LQUERY = ( LWORK.EQ.-1 ) 00160 IF( N.LT.0 ) THEN 00161 INFO = -1 00162 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00163 INFO = -3 00164 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN 00165 INFO = -6 00166 END IF 00167 IF( INFO.NE.0 ) THEN 00168 CALL XERBLA( 'SGETRI', -INFO ) 00169 RETURN 00170 ELSE IF( LQUERY ) THEN 00171 RETURN 00172 END IF 00173 * 00174 * Quick return if possible 00175 * 00176 IF( N.EQ.0 ) 00177 $ RETURN 00178 * 00179 * Form inv(U). If INFO > 0 from STRTRI, then U is singular, 00180 * and the inverse is not computed. 00181 * 00182 CALL STRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO ) 00183 IF( INFO.GT.0 ) 00184 $ RETURN 00185 * 00186 NBMIN = 2 00187 LDWORK = N 00188 IF( NB.GT.1 .AND. NB.LT.N ) THEN 00189 IWS = MAX( LDWORK*NB, 1 ) 00190 IF( LWORK.LT.IWS ) THEN 00191 NB = LWORK / LDWORK 00192 NBMIN = MAX( 2, ILAENV( 2, 'SGETRI', ' ', N, -1, -1, -1 ) ) 00193 END IF 00194 ELSE 00195 IWS = N 00196 END IF 00197 * 00198 * Solve the equation inv(A)*L = inv(U) for inv(A). 00199 * 00200 IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN 00201 * 00202 * Use unblocked code. 00203 * 00204 DO 20 J = N, 1, -1 00205 * 00206 * Copy current column of L to WORK and replace with zeros. 00207 * 00208 DO 10 I = J + 1, N 00209 WORK( I ) = A( I, J ) 00210 A( I, J ) = ZERO 00211 10 CONTINUE 00212 * 00213 * Compute current column of inv(A). 00214 * 00215 IF( J.LT.N ) 00216 $ CALL SGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ), 00217 $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) 00218 20 CONTINUE 00219 ELSE 00220 * 00221 * Use blocked code. 00222 * 00223 NN = ( ( N-1 ) / NB )*NB + 1 00224 DO 50 J = NN, 1, -NB 00225 JB = MIN( NB, N-J+1 ) 00226 * 00227 * Copy current block column of L to WORK and replace with 00228 * zeros. 00229 * 00230 DO 40 JJ = J, J + JB - 1 00231 DO 30 I = JJ + 1, N 00232 WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) 00233 A( I, JJ ) = ZERO 00234 30 CONTINUE 00235 40 CONTINUE 00236 * 00237 * Compute current block column of inv(A). 00238 * 00239 IF( J+JB.LE.N ) 00240 $ CALL SGEMM( 'No transpose', 'No transpose', N, JB, 00241 $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA, 00242 $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA ) 00243 CALL STRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB, 00244 $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) 00245 50 CONTINUE 00246 END IF 00247 * 00248 * Apply column interchanges. 00249 * 00250 DO 60 J = N - 1, 1, -1 00251 JP = IPIV( J ) 00252 IF( JP.NE.J ) 00253 $ CALL SSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 ) 00254 60 CONTINUE 00255 * 00256 WORK( 1 ) = IWS 00257 RETURN 00258 * 00259 * End of SGETRI 00260 * 00261 END