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