LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgetri.f
Go to the documentation of this file.
1*> \brief \b SGETRI
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGETRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgetri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgetri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgetri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, LDA, LWORK, N
25* ..
26* .. Array Arguments ..
27* INTEGER IPIV( * )
28* REAL A( LDA, * ), WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SGETRI computes the inverse of a matrix using the LU factorization
38*> computed by SGETRF.
39*>
40*> This method inverts U and then computes inv(A) by solving the system
41*> inv(A)*L = inv(U) for inv(A).
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] N
48*> \verbatim
49*> N is INTEGER
50*> The order of the matrix A. N >= 0.
51*> \endverbatim
52*>
53*> \param[in,out] A
54*> \verbatim
55*> A is REAL array, dimension (LDA,N)
56*> On entry, the factors L and U from the factorization
57*> A = P*L*U as computed by SGETRF.
58*> On exit, if INFO = 0, the inverse of the original matrix A.
59*> \endverbatim
60*>
61*> \param[in] LDA
62*> \verbatim
63*> LDA is INTEGER
64*> The leading dimension of the array A. LDA >= max(1,N).
65*> \endverbatim
66*>
67*> \param[in] IPIV
68*> \verbatim
69*> IPIV is INTEGER array, dimension (N)
70*> The pivot indices from SGETRF; for 1<=i<=N, row i of the
71*> matrix was interchanged with row IPIV(i).
72*> \endverbatim
73*>
74*> \param[out] WORK
75*> \verbatim
76*> WORK is REAL array, dimension (MAX(1,LWORK))
77*> On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
78*> \endverbatim
79*>
80*> \param[in] LWORK
81*> \verbatim
82*> LWORK is INTEGER
83*> The dimension of the array WORK. LWORK >= max(1,N).
84*> For optimal performance LWORK >= N*NB, where NB is
85*> the optimal blocksize returned by ILAENV.
86*>
87*> If LWORK = -1, then a workspace query is assumed; the routine
88*> only calculates the optimal size of the WORK array, returns
89*> this value as the first entry of the WORK array, and no error
90*> message related to LWORK is issued by XERBLA.
91*> \endverbatim
92*>
93*> \param[out] INFO
94*> \verbatim
95*> INFO is INTEGER
96*> = 0: successful exit
97*> < 0: if INFO = -i, the i-th argument had an illegal value
98*> > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
99*> singular and its inverse could not be computed.
100*> \endverbatim
101*
102* Authors:
103* ========
104*
105*> \author Univ. of Tennessee
106*> \author Univ. of California Berkeley
107*> \author Univ. of Colorado Denver
108*> \author NAG Ltd.
109*
110*> \ingroup realGEcomputational
111*
112* =====================================================================
113 SUBROUTINE sgetri( N, A, LDA, IPIV, WORK, LWORK, INFO )
114*
115* -- LAPACK computational routine --
116* -- LAPACK is a software package provided by Univ. of Tennessee, --
117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118*
119* .. Scalar Arguments ..
120 INTEGER INFO, LDA, LWORK, N
121* ..
122* .. Array Arguments ..
123 INTEGER IPIV( * )
124 REAL A( LDA, * ), WORK( * )
125* ..
126*
127* =====================================================================
128*
129* .. Parameters ..
130 REAL ZERO, ONE
131 parameter( zero = 0.0e+0, one = 1.0e+0 )
132* ..
133* .. Local Scalars ..
134 LOGICAL LQUERY
135 INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
136 $ NBMIN, NN
137* ..
138* .. External Functions ..
139 INTEGER ILAENV
140 EXTERNAL ilaenv
141* ..
142* .. External Subroutines ..
143 EXTERNAL sgemm, sgemv, sswap, strsm, strtri, xerbla
144* ..
145* .. Intrinsic Functions ..
146 INTRINSIC max, min
147* ..
148* .. Executable Statements ..
149*
150* Test the input parameters.
151*
152 info = 0
153 nb = ilaenv( 1, 'SGETRI', ' ', n, -1, -1, -1 )
154 lwkopt = n*nb
155 work( 1 ) = lwkopt
156 lquery = ( lwork.EQ.-1 )
157 IF( n.LT.0 ) THEN
158 info = -1
159 ELSE IF( lda.LT.max( 1, n ) ) THEN
160 info = -3
161 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
162 info = -6
163 END IF
164 IF( info.NE.0 ) THEN
165 CALL xerbla( 'SGETRI', -info )
166 RETURN
167 ELSE IF( lquery ) THEN
168 RETURN
169 END IF
170*
171* Quick return if possible
172*
173 IF( n.EQ.0 )
174 $ RETURN
175*
176* Form inv(U). If INFO > 0 from STRTRI, then U is singular,
177* and the inverse is not computed.
178*
179 CALL strtri( 'Upper', 'Non-unit', n, a, lda, info )
180 IF( info.GT.0 )
181 $ RETURN
182*
183 nbmin = 2
184 ldwork = n
185 IF( nb.GT.1 .AND. nb.LT.n ) THEN
186 iws = max( ldwork*nb, 1 )
187 IF( lwork.LT.iws ) THEN
188 nb = lwork / ldwork
189 nbmin = max( 2, ilaenv( 2, 'SGETRI', ' ', n, -1, -1, -1 ) )
190 END IF
191 ELSE
192 iws = n
193 END IF
194*
195* Solve the equation inv(A)*L = inv(U) for inv(A).
196*
197 IF( nb.LT.nbmin .OR. nb.GE.n ) THEN
198*
199* Use unblocked code.
200*
201 DO 20 j = n, 1, -1
202*
203* Copy current column of L to WORK and replace with zeros.
204*
205 DO 10 i = j + 1, n
206 work( i ) = a( i, j )
207 a( i, j ) = zero
208 10 CONTINUE
209*
210* Compute current column of inv(A).
211*
212 IF( j.LT.n )
213 $ CALL sgemv( 'No transpose', n, n-j, -one, a( 1, j+1 ),
214 $ lda, work( j+1 ), 1, one, a( 1, j ), 1 )
215 20 CONTINUE
216 ELSE
217*
218* Use blocked code.
219*
220 nn = ( ( n-1 ) / nb )*nb + 1
221 DO 50 j = nn, 1, -nb
222 jb = min( nb, n-j+1 )
223*
224* Copy current block column of L to WORK and replace with
225* zeros.
226*
227 DO 40 jj = j, j + jb - 1
228 DO 30 i = jj + 1, n
229 work( i+( jj-j )*ldwork ) = a( i, jj )
230 a( i, jj ) = zero
231 30 CONTINUE
232 40 CONTINUE
233*
234* Compute current block column of inv(A).
235*
236 IF( j+jb.LE.n )
237 $ CALL sgemm( 'No transpose', 'No transpose', n, jb,
238 $ n-j-jb+1, -one, a( 1, j+jb ), lda,
239 $ work( j+jb ), ldwork, one, a( 1, j ), lda )
240 CALL strsm( 'Right', 'Lower', 'No transpose', 'Unit', n, jb,
241 $ one, work( j ), ldwork, a( 1, j ), lda )
242 50 CONTINUE
243 END IF
244*
245* Apply column interchanges.
246*
247 DO 60 j = n - 1, 1, -1
248 jp = ipiv( j )
249 IF( jp.NE.j )
250 $ CALL sswap( n, a( 1, j ), 1, a( 1, jp ), 1 )
251 60 CONTINUE
252*
253 work( 1 ) = iws
254 RETURN
255*
256* End of SGETRI
257*
258 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sgetri(N, A, LDA, IPIV, WORK, LWORK, INFO)
SGETRI
Definition: sgetri.f:114
subroutine strtri(UPLO, DIAG, N, A, LDA, INFO)
STRTRI
Definition: strtri.f:109
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:181
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187