LAPACK 3.12.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 getri
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 REAL SROUNDUP_LWORK
141 EXTERNAL ilaenv, sroundup_lwork
142* ..
143* .. External Subroutines ..
144 EXTERNAL sgemm, sgemv, sswap, strsm, strtri, xerbla
145* ..
146* .. Intrinsic Functions ..
147 INTRINSIC max, min
148* ..
149* .. Executable Statements ..
150*
151* Test the input parameters.
152*
153 info = 0
154 nb = ilaenv( 1, 'SGETRI', ' ', n, -1, -1, -1 )
155 lwkopt = n*nb
156 work( 1 ) = sroundup_lwork(lwkopt)
157 lquery = ( lwork.EQ.-1 )
158 IF( n.LT.0 ) THEN
159 info = -1
160 ELSE IF( lda.LT.max( 1, n ) ) THEN
161 info = -3
162 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
163 info = -6
164 END IF
165 IF( info.NE.0 ) THEN
166 CALL xerbla( 'SGETRI', -info )
167 RETURN
168 ELSE IF( lquery ) THEN
169 RETURN
170 END IF
171*
172* Quick return if possible
173*
174 IF( n.EQ.0 )
175 $ RETURN
176*
177* Form inv(U). If INFO > 0 from STRTRI, then U is singular,
178* and the inverse is not computed.
179*
180 CALL strtri( 'Upper', 'Non-unit', n, a, lda, info )
181 IF( info.GT.0 )
182 $ RETURN
183*
184 nbmin = 2
185 ldwork = n
186 IF( nb.GT.1 .AND. nb.LT.n ) THEN
187 iws = max( ldwork*nb, 1 )
188 IF( lwork.LT.iws ) THEN
189 nb = lwork / ldwork
190 nbmin = max( 2, ilaenv( 2, 'SGETRI', ' ', n, -1, -1, -1 ) )
191 END IF
192 ELSE
193 iws = n
194 END IF
195*
196* Solve the equation inv(A)*L = inv(U) for inv(A).
197*
198 IF( nb.LT.nbmin .OR. nb.GE.n ) THEN
199*
200* Use unblocked code.
201*
202 DO 20 j = n, 1, -1
203*
204* Copy current column of L to WORK and replace with zeros.
205*
206 DO 10 i = j + 1, n
207 work( i ) = a( i, j )
208 a( i, j ) = zero
209 10 CONTINUE
210*
211* Compute current column of inv(A).
212*
213 IF( j.LT.n )
214 $ CALL sgemv( 'No transpose', n, n-j, -one, a( 1, j+1 ),
215 $ lda, work( j+1 ), 1, one, a( 1, j ), 1 )
216 20 CONTINUE
217 ELSE
218*
219* Use blocked code.
220*
221 nn = ( ( n-1 ) / nb )*nb + 1
222 DO 50 j = nn, 1, -nb
223 jb = min( nb, n-j+1 )
224*
225* Copy current block column of L to WORK and replace with
226* zeros.
227*
228 DO 40 jj = j, j + jb - 1
229 DO 30 i = jj + 1, n
230 work( i+( jj-j )*ldwork ) = a( i, jj )
231 a( i, jj ) = zero
232 30 CONTINUE
233 40 CONTINUE
234*
235* Compute current block column of inv(A).
236*
237 IF( j+jb.LE.n )
238 $ CALL sgemm( 'No transpose', 'No transpose', n, jb,
239 $ n-j-jb+1, -one, a( 1, j+jb ), lda,
240 $ work( j+jb ), ldwork, one, a( 1, j ), lda )
241 CALL strsm( 'Right', 'Lower', 'No transpose', 'Unit', n, jb,
242 $ one, work( j ), ldwork, a( 1, j ), lda )
243 50 CONTINUE
244 END IF
245*
246* Apply column interchanges.
247*
248 DO 60 j = n - 1, 1, -1
249 jp = ipiv( j )
250 IF( jp.NE.j )
251 $ CALL sswap( n, a( 1, j ), 1, a( 1, jp ), 1 )
252 60 CONTINUE
253*
254 work( 1 ) = sroundup_lwork(iws)
255 RETURN
256*
257* End of SGETRI
258*
259 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
subroutine sgetri(n, a, lda, ipiv, work, lwork, info)
SGETRI
Definition sgetri.f:114
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
subroutine strtri(uplo, diag, n, a, lda, info)
STRTRI
Definition strtri.f:109