LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
spftri.f
Go to the documentation of this file.
1*> \brief \b SPFTRI
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SPFTRI + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spftri.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spftri.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spftri.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SPFTRI( TRANSR, UPLO, N, A, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER TRANSR, UPLO
25* INTEGER INFO, N
26* .. Array Arguments ..
27* REAL A( 0: * )
28* ..
29*
30*
31*> \par Purpose:
32* =============
33*>
34*> \verbatim
35*>
36*> SPFTRI computes the inverse of a real (symmetric) positive definite
37*> matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
38*> computed by SPFTRF.
39*> \endverbatim
40*
41* Arguments:
42* ==========
43*
44*> \param[in] TRANSR
45*> \verbatim
46*> TRANSR is CHARACTER*1
47*> = 'N': The Normal TRANSR of RFP A is stored;
48*> = 'T': The Transpose TRANSR of RFP A is stored.
49*> \endverbatim
50*>
51*> \param[in] UPLO
52*> \verbatim
53*> UPLO is CHARACTER*1
54*> = 'U': Upper triangle of A is stored;
55*> = 'L': Lower triangle of A is stored.
56*> \endverbatim
57*>
58*> \param[in] N
59*> \verbatim
60*> N is INTEGER
61*> The order of the matrix A. N >= 0.
62*> \endverbatim
63*>
64*> \param[in,out] A
65*> \verbatim
66*> A is REAL array, dimension ( N*(N+1)/2 )
67*> On entry, the symmetric matrix A in RFP format. RFP format is
68*> described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
69*> then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
70*> (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
71*> the transpose of RFP A as defined when
72*> TRANSR = 'N'. The contents of RFP A are defined by UPLO as
73*> follows: If UPLO = 'U' the RFP A contains the nt elements of
74*> upper packed A. If UPLO = 'L' the RFP A contains the elements
75*> of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
76*> 'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
77*> is odd. See the Note below for more details.
78*>
79*> On exit, the symmetric inverse of the original matrix, in the
80*> same storage format.
81*> \endverbatim
82*>
83*> \param[out] INFO
84*> \verbatim
85*> INFO is INTEGER
86*> = 0: successful exit
87*> < 0: if INFO = -i, the i-th argument had an illegal value
88*> > 0: if INFO = i, the (i,i) element of the factor U or L is
89*> zero, and the inverse could not be computed.
90*> \endverbatim
91*
92* Authors:
93* ========
94*
95*> \author Univ. of Tennessee
96*> \author Univ. of California Berkeley
97*> \author Univ. of Colorado Denver
98*> \author NAG Ltd.
99*
100*> \ingroup pftri
101*
102*> \par Further Details:
103* =====================
104*>
105*> \verbatim
106*>
107*> We first consider Rectangular Full Packed (RFP) Format when N is
108*> even. We give an example where N = 6.
109*>
110*> AP is Upper AP is Lower
111*>
112*> 00 01 02 03 04 05 00
113*> 11 12 13 14 15 10 11
114*> 22 23 24 25 20 21 22
115*> 33 34 35 30 31 32 33
116*> 44 45 40 41 42 43 44
117*> 55 50 51 52 53 54 55
118*>
119*>
120*> Let TRANSR = 'N'. RFP holds AP as follows:
121*> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
122*> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
123*> the transpose of the first three columns of AP upper.
124*> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
125*> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
126*> the transpose of the last three columns of AP lower.
127*> This covers the case N even and TRANSR = 'N'.
128*>
129*> RFP A RFP A
130*>
131*> 03 04 05 33 43 53
132*> 13 14 15 00 44 54
133*> 23 24 25 10 11 55
134*> 33 34 35 20 21 22
135*> 00 44 45 30 31 32
136*> 01 11 55 40 41 42
137*> 02 12 22 50 51 52
138*>
139*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
140*> transpose of RFP A above. One therefore gets:
141*>
142*>
143*> RFP A RFP A
144*>
145*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
146*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
147*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
148*>
149*>
150*> We then consider Rectangular Full Packed (RFP) Format when N is
151*> odd. We give an example where N = 5.
152*>
153*> AP is Upper AP is Lower
154*>
155*> 00 01 02 03 04 00
156*> 11 12 13 14 10 11
157*> 22 23 24 20 21 22
158*> 33 34 30 31 32 33
159*> 44 40 41 42 43 44
160*>
161*>
162*> Let TRANSR = 'N'. RFP holds AP as follows:
163*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
164*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
165*> the transpose of the first two columns of AP upper.
166*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
167*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
168*> the transpose of the last two columns of AP lower.
169*> This covers the case N odd and TRANSR = 'N'.
170*>
171*> RFP A RFP A
172*>
173*> 02 03 04 00 33 43
174*> 12 13 14 10 11 44
175*> 22 23 24 20 21 22
176*> 00 33 34 30 31 32
177*> 01 11 44 40 41 42
178*>
179*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
180*> transpose of RFP A above. One therefore gets:
181*>
182*> RFP A RFP A
183*>
184*> 02 12 22 00 01 00 10 20 30 40 50
185*> 03 13 23 33 11 33 11 21 31 41 51
186*> 04 14 24 34 44 43 44 22 32 42 52
187*> \endverbatim
188*>
189* =====================================================================
190 SUBROUTINE spftri( TRANSR, UPLO, N, A, INFO )
191*
192* -- LAPACK computational routine --
193* -- LAPACK is a software package provided by Univ. of Tennessee, --
194* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195*
196* .. Scalar Arguments ..
197 CHARACTER TRANSR, UPLO
198 INTEGER INFO, N
199* .. Array Arguments ..
200 REAL A( 0: * )
201* ..
202*
203* =====================================================================
204*
205* .. Parameters ..
206 REAL ONE
207 parameter( one = 1.0e+0 )
208* ..
209* .. Local Scalars ..
210 LOGICAL LOWER, NISODD, NORMALTRANSR
211 INTEGER N1, N2, K
212* ..
213* .. External Functions ..
214 LOGICAL LSAME
215 EXTERNAL lsame
216* ..
217* .. External Subroutines ..
218 EXTERNAL xerbla, stftri, slauum, strmm, ssyrk
219* ..
220* .. Intrinsic Functions ..
221 INTRINSIC mod
222* ..
223* .. Executable Statements ..
224*
225* Test the input parameters.
226*
227 info = 0
228 normaltransr = lsame( transr, 'N' )
229 lower = lsame( uplo, 'L' )
230 IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
231 info = -1
232 ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
233 info = -2
234 ELSE IF( n.LT.0 ) THEN
235 info = -3
236 END IF
237 IF( info.NE.0 ) THEN
238 CALL xerbla( 'SPFTRI', -info )
239 RETURN
240 END IF
241*
242* Quick return if possible
243*
244 IF( n.EQ.0 )
245 $ RETURN
246*
247* Invert the triangular Cholesky factor U or L.
248*
249 CALL stftri( transr, uplo, 'N', n, a, info )
250 IF( info.GT.0 )
251 $ RETURN
252*
253* If N is odd, set NISODD = .TRUE.
254* If N is even, set K = N/2 and NISODD = .FALSE.
255*
256 IF( mod( n, 2 ).EQ.0 ) THEN
257 k = n / 2
258 nisodd = .false.
259 ELSE
260 nisodd = .true.
261 END IF
262*
263* Set N1 and N2 depending on LOWER
264*
265 IF( lower ) THEN
266 n2 = n / 2
267 n1 = n - n2
268 ELSE
269 n1 = n / 2
270 n2 = n - n1
271 END IF
272*
273* Start execution of triangular matrix multiply: inv(U)*inv(U)^C or
274* inv(L)^C*inv(L). There are eight cases.
275*
276 IF( nisodd ) THEN
277*
278* N is odd
279*
280 IF( normaltransr ) THEN
281*
282* N is odd and TRANSR = 'N'
283*
284 IF( lower ) THEN
285*
286* SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:N1-1) )
287* T1 -> a(0,0), T2 -> a(0,1), S -> a(N1,0)
288* T1 -> a(0), T2 -> a(n), S -> a(N1)
289*
290 CALL slauum( 'L', n1, a( 0 ), n, info )
291 CALL ssyrk( 'L', 'T', n1, n2, one, a( n1 ), n, one,
292 $ a( 0 ), n )
293 CALL strmm( 'L', 'U', 'N', 'N', n2, n1, one, a( n ), n,
294 $ a( n1 ), n )
295 CALL slauum( 'U', n2, a( n ), n, info )
296*
297 ELSE
298*
299* SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:N2-1)
300* T1 -> a(N1+1,0), T2 -> a(N1,0), S -> a(0,0)
301* T1 -> a(N2), T2 -> a(N1), S -> a(0)
302*
303 CALL slauum( 'L', n1, a( n2 ), n, info )
304 CALL ssyrk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
305 $ a( n2 ), n )
306 CALL strmm( 'R', 'U', 'T', 'N', n1, n2, one, a( n1 ), n,
307 $ a( 0 ), n )
308 CALL slauum( 'U', n2, a( n1 ), n, info )
309*
310 END IF
311*
312 ELSE
313*
314* N is odd and TRANSR = 'T'
315*
316 IF( lower ) THEN
317*
318* SRPA for LOWER, TRANSPOSE, and N is odd
319* T1 -> a(0), T2 -> a(1), S -> a(0+N1*N1)
320*
321 CALL slauum( 'U', n1, a( 0 ), n1, info )
322 CALL ssyrk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
323 $ a( 0 ), n1 )
324 CALL strmm( 'R', 'L', 'N', 'N', n1, n2, one, a( 1 ), n1,
325 $ a( n1*n1 ), n1 )
326 CALL slauum( 'L', n2, a( 1 ), n1, info )
327*
328 ELSE
329*
330* SRPA for UPPER, TRANSPOSE, and N is odd
331* T1 -> a(0+N2*N2), T2 -> a(0+N1*N2), S -> a(0)
332*
333 CALL slauum( 'U', n1, a( n2*n2 ), n2, info )
334 CALL ssyrk( 'U', 'T', n1, n2, one, a( 0 ), n2, one,
335 $ a( n2*n2 ), n2 )
336 CALL strmm( 'L', 'L', 'T', 'N', n2, n1, one, a( n1*n2 ),
337 $ n2, a( 0 ), n2 )
338 CALL slauum( 'L', n2, a( n1*n2 ), n2, info )
339*
340 END IF
341*
342 END IF
343*
344 ELSE
345*
346* N is even
347*
348 IF( normaltransr ) THEN
349*
350* N is even and TRANSR = 'N'
351*
352 IF( lower ) THEN
353*
354* SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
355* T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
356* T1 -> a(1), T2 -> a(0), S -> a(k+1)
357*
358 CALL slauum( 'L', k, a( 1 ), n+1, info )
359 CALL ssyrk( 'L', 'T', k, k, one, a( k+1 ), n+1, one,
360 $ a( 1 ), n+1 )
361 CALL strmm( 'L', 'U', 'N', 'N', k, k, one, a( 0 ), n+1,
362 $ a( k+1 ), n+1 )
363 CALL slauum( 'U', k, a( 0 ), n+1, info )
364*
365 ELSE
366*
367* SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
368* T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
369* T1 -> a(k+1), T2 -> a(k), S -> a(0)
370*
371 CALL slauum( 'L', k, a( k+1 ), n+1, info )
372 CALL ssyrk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
373 $ a( k+1 ), n+1 )
374 CALL strmm( 'R', 'U', 'T', 'N', k, k, one, a( k ), n+1,
375 $ a( 0 ), n+1 )
376 CALL slauum( 'U', k, a( k ), n+1, info )
377*
378 END IF
379*
380 ELSE
381*
382* N is even and TRANSR = 'T'
383*
384 IF( lower ) THEN
385*
386* SRPA for LOWER, TRANSPOSE, and N is even (see paper)
387* T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1),
388* T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
389*
390 CALL slauum( 'U', k, a( k ), k, info )
391 CALL ssyrk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
392 $ a( k ), k )
393 CALL strmm( 'R', 'L', 'N', 'N', k, k, one, a( 0 ), k,
394 $ a( k*( k+1 ) ), k )
395 CALL slauum( 'L', k, a( 0 ), k, info )
396*
397 ELSE
398*
399* SRPA for UPPER, TRANSPOSE, and N is even (see paper)
400* T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0),
401* T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
402*
403 CALL slauum( 'U', k, a( k*( k+1 ) ), k, info )
404 CALL ssyrk( 'U', 'T', k, k, one, a( 0 ), k, one,
405 $ a( k*( k+1 ) ), k )
406 CALL strmm( 'L', 'L', 'T', 'N', k, k, one, a( k*k ), k,
407 $ a( 0 ), k )
408 CALL slauum( 'L', k, a( k*k ), k, info )
409*
410 END IF
411*
412 END IF
413*
414 END IF
415*
416 RETURN
417*
418* End of SPFTRI
419*
420 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
subroutine slauum(uplo, n, a, lda, info)
SLAUUM computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked...
Definition slauum.f:102
subroutine spftri(transr, uplo, n, a, info)
SPFTRI
Definition spftri.f:191
subroutine stftri(transr, uplo, diag, n, a, info)
STFTRI
Definition stftri.f:201
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177