LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \htmlonly
9 *> Download DTFTRI + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtftri.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtftri.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtftri.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTFTRI( TRANSR, UPLO, DIAG, N, A, INFO )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER TRANSR, UPLO, DIAG
25 * INTEGER INFO, N
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION A( 0: * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> DTFTRI 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 DOUBLE PRECISION array, dimension (0:nt-1);
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 *> \date November 2011
111 *
112 *> \ingroup doubleOTHERcomputational
113 *
114 *> \par Further Details:
115 * =====================
116 *>
117 *> \verbatim
118 *>
119 *> We first consider Rectangular Full Packed (RFP) Format when N is
120 *> even. We give an example where N = 6.
121 *>
122 *> AP is Upper AP is Lower
123 *>
124 *> 00 01 02 03 04 05 00
125 *> 11 12 13 14 15 10 11
126 *> 22 23 24 25 20 21 22
127 *> 33 34 35 30 31 32 33
128 *> 44 45 40 41 42 43 44
129 *> 55 50 51 52 53 54 55
130 *>
131 *>
132 *> Let TRANSR = 'N'. RFP holds AP as follows:
133 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
134 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
135 *> the transpose of the first three columns of AP upper.
136 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
137 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
138 *> the transpose of the last three columns of AP lower.
139 *> This covers the case N even and TRANSR = 'N'.
140 *>
141 *> RFP A RFP A
142 *>
143 *> 03 04 05 33 43 53
144 *> 13 14 15 00 44 54
145 *> 23 24 25 10 11 55
146 *> 33 34 35 20 21 22
147 *> 00 44 45 30 31 32
148 *> 01 11 55 40 41 42
149 *> 02 12 22 50 51 52
150 *>
151 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
152 *> transpose of RFP A above. One therefore gets:
153 *>
154 *>
155 *> RFP A RFP A
156 *>
157 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
158 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
159 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
160 *>
161 *>
162 *> We then consider Rectangular Full Packed (RFP) Format when N is
163 *> odd. We give an example where N = 5.
164 *>
165 *> AP is Upper AP is Lower
166 *>
167 *> 00 01 02 03 04 00
168 *> 11 12 13 14 10 11
169 *> 22 23 24 20 21 22
170 *> 33 34 30 31 32 33
171 *> 44 40 41 42 43 44
172 *>
173 *>
174 *> Let TRANSR = 'N'. RFP holds AP as follows:
175 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
176 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
177 *> the transpose of the first two columns of AP upper.
178 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
179 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
180 *> the transpose of the last two columns of AP lower.
181 *> This covers the case N odd and TRANSR = 'N'.
182 *>
183 *> RFP A RFP A
184 *>
185 *> 02 03 04 00 33 43
186 *> 12 13 14 10 11 44
187 *> 22 23 24 20 21 22
188 *> 00 33 34 30 31 32
189 *> 01 11 44 40 41 42
190 *>
191 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
192 *> transpose of RFP A above. One therefore gets:
193 *>
194 *> RFP A RFP A
195 *>
196 *> 02 12 22 00 01 00 10 20 30 40 50
197 *> 03 13 23 33 11 33 11 21 31 41 51
198 *> 04 14 24 34 44 43 44 22 32 42 52
199 *> \endverbatim
200 *>
201 * =====================================================================
202  SUBROUTINE dtftri( TRANSR, UPLO, DIAG, N, A, INFO )
203 *
204 * -- LAPACK computational routine (version 3.4.0) --
205 * -- LAPACK is a software package provided by Univ. of Tennessee, --
206 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
207 * November 2011
208 *
209 * .. Scalar Arguments ..
210  CHARACTER transr, uplo, diag
211  INTEGER info, n
212 * ..
213 * .. Array Arguments ..
214  DOUBLE PRECISION a( 0: * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  DOUBLE PRECISION one
221  parameter( one = 1.0d+0 )
222 * ..
223 * .. Local Scalars ..
224  LOGICAL lower, nisodd, normaltransr
225  INTEGER n1, n2, k
226 * ..
227 * .. External Functions ..
228  LOGICAL lsame
229  EXTERNAL lsame
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL xerbla, dtrmm, dtrtri
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC mod
236 * ..
237 * .. Executable Statements ..
238 *
239 * Test the input parameters.
240 *
241  info = 0
242  normaltransr = lsame( transr, 'N' )
243  lower = lsame( uplo, 'L' )
244  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
245  info = -1
246  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
247  info = -2
248  ELSE IF( .NOT.lsame( diag, 'N' ) .AND. .NOT.lsame( diag, 'U' ) )
249  $ THEN
250  info = -3
251  ELSE IF( n.LT.0 ) THEN
252  info = -4
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'DTFTRI', -info )
256  return
257  END IF
258 *
259 * Quick return if possible
260 *
261  IF( n.EQ.0 )
262  $ return
263 *
264 * If N is odd, set NISODD = .TRUE.
265 * If N is even, set K = N/2 and NISODD = .FALSE.
266 *
267  IF( mod( n, 2 ).EQ.0 ) THEN
268  k = n / 2
269  nisodd = .false.
270  ELSE
271  nisodd = .true.
272  END IF
273 *
274 * Set N1 and N2 depending on LOWER
275 *
276  IF( lower ) THEN
277  n2 = n / 2
278  n1 = n - n2
279  ELSE
280  n1 = n / 2
281  n2 = n - n1
282  END IF
283 *
284 *
285 * start execution: there are eight cases
286 *
287  IF( nisodd ) THEN
288 *
289 * N is odd
290 *
291  IF( normaltransr ) THEN
292 *
293 * N is odd and TRANSR = 'N'
294 *
295  IF( lower ) THEN
296 *
297 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
298 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
299 * T1 -> a(0), T2 -> a(n), S -> a(n1)
300 *
301  CALL dtrtri( 'L', diag, n1, a( 0 ), n, info )
302  IF( info.GT.0 )
303  $ return
304  CALL dtrmm( 'R', 'L', 'N', diag, n2, n1, -one, a( 0 ),
305  $ n, a( n1 ), n )
306  CALL dtrtri( 'U', diag, n2, a( n ), n, info )
307  IF( info.GT.0 )
308  $ info = info + n1
309  IF( info.GT.0 )
310  $ return
311  CALL dtrmm( 'L', 'U', 'T', diag, n2, n1, one, a( n ), n,
312  $ a( n1 ), n )
313 *
314  ELSE
315 *
316 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
317 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
318 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
319 *
320  CALL dtrtri( 'L', diag, n1, a( n2 ), n, info )
321  IF( info.GT.0 )
322  $ return
323  CALL dtrmm( 'L', 'L', 'T', diag, n1, n2, -one, a( n2 ),
324  $ n, a( 0 ), n )
325  CALL dtrtri( 'U', diag, n2, a( n1 ), n, info )
326  IF( info.GT.0 )
327  $ info = info + n1
328  IF( info.GT.0 )
329  $ return
330  CALL dtrmm( 'R', 'U', 'N', diag, n1, n2, one, a( n1 ),
331  $ n, a( 0 ), n )
332 *
333  END IF
334 *
335  ELSE
336 *
337 * N is odd and TRANSR = 'T'
338 *
339  IF( lower ) THEN
340 *
341 * SRPA for LOWER, TRANSPOSE and N is odd
342 * T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
343 *
344  CALL dtrtri( 'U', diag, n1, a( 0 ), n1, info )
345  IF( info.GT.0 )
346  $ return
347  CALL dtrmm( 'L', 'U', 'N', diag, n1, n2, -one, a( 0 ),
348  $ n1, a( n1*n1 ), n1 )
349  CALL dtrtri( 'L', diag, n2, a( 1 ), n1, info )
350  IF( info.GT.0 )
351  $ info = info + n1
352  IF( info.GT.0 )
353  $ return
354  CALL dtrmm( 'R', 'L', 'T', diag, n1, n2, one, a( 1 ),
355  $ n1, a( n1*n1 ), n1 )
356 *
357  ELSE
358 *
359 * SRPA for UPPER, TRANSPOSE and N is odd
360 * T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
361 *
362  CALL dtrtri( 'U', diag, n1, a( n2*n2 ), n2, info )
363  IF( info.GT.0 )
364  $ return
365  CALL dtrmm( 'R', 'U', 'T', diag, n2, n1, -one,
366  $ a( n2*n2 ), n2, a( 0 ), n2 )
367  CALL dtrtri( 'L', diag, n2, a( n1*n2 ), n2, info )
368  IF( info.GT.0 )
369  $ info = info + n1
370  IF( info.GT.0 )
371  $ return
372  CALL dtrmm( 'L', 'L', 'N', diag, n2, n1, one,
373  $ a( n1*n2 ), n2, a( 0 ), n2 )
374  END IF
375 *
376  END IF
377 *
378  ELSE
379 *
380 * N is even
381 *
382  IF( normaltransr ) THEN
383 *
384 * N is even and TRANSR = 'N'
385 *
386  IF( lower ) THEN
387 *
388 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
389 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
390 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
391 *
392  CALL dtrtri( 'L', diag, k, a( 1 ), n+1, info )
393  IF( info.GT.0 )
394  $ return
395  CALL dtrmm( 'R', 'L', 'N', diag, k, k, -one, a( 1 ),
396  $ n+1, a( k+1 ), n+1 )
397  CALL dtrtri( 'U', diag, k, a( 0 ), n+1, info )
398  IF( info.GT.0 )
399  $ info = info + k
400  IF( info.GT.0 )
401  $ return
402  CALL dtrmm( 'L', 'U', 'T', diag, k, k, one, a( 0 ), n+1,
403  $ a( k+1 ), n+1 )
404 *
405  ELSE
406 *
407 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
408 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
409 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
410 *
411  CALL dtrtri( 'L', diag, k, a( k+1 ), n+1, info )
412  IF( info.GT.0 )
413  $ return
414  CALL dtrmm( 'L', 'L', 'T', diag, k, k, -one, a( k+1 ),
415  $ n+1, a( 0 ), n+1 )
416  CALL dtrtri( 'U', diag, k, a( k ), n+1, info )
417  IF( info.GT.0 )
418  $ info = info + k
419  IF( info.GT.0 )
420  $ return
421  CALL dtrmm( 'R', 'U', 'N', diag, k, k, one, a( k ), 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 ), k,
438  $ a( k*( k+1 ) ), k )
439  CALL dtrtri( 'L', diag, k, a( 0 ), k, info )
440  IF( info.GT.0 )
441  $ info = info + k
442  IF( info.GT.0 )
443  $ return
444  CALL dtrmm( 'R', 'L', 'T', diag, k, k, one, a( 0 ), k,
445  $ a( k*( k+1 ) ), k )
446  ELSE
447 *
448 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
449 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
450 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
451 *
452  CALL dtrtri( 'U', diag, k, a( k*( k+1 ) ), k, info )
453  IF( info.GT.0 )
454  $ return
455  CALL dtrmm( 'R', 'U', 'T', diag, k, k, -one,
456  $ a( k*( k+1 ) ), k, a( 0 ), k )
457  CALL dtrtri( 'L', diag, k, a( k*k ), k, info )
458  IF( info.GT.0 )
459  $ info = info + k
460  IF( info.GT.0 )
461  $ return
462  CALL dtrmm( 'L', 'L', 'N', diag, k, k, one, a( k*k ), k,
463  $ a( 0 ), k )
464  END IF
465  END IF
466  END IF
467 *
468  return
469 *
470 * End of DTFTRI
471 *
472  END