LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 realOTHERcomputational
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)
XERBLA
Definition: xerbla.f:60
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
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:169