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