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