LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ spftri()

subroutine spftri ( character  TRANSR,
character  UPLO,
integer  N,
real, dimension( 0: * )  A,
integer  INFO 
)

SPFTRI

Download SPFTRI + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SPFTRI computes the inverse of a real (symmetric) positive definite
 matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
 computed by SPFTRF.
Parameters
[in]TRANSR
          TRANSR is CHARACTER*1
          = 'N':  The Normal TRANSR of RFP A is stored;
          = 'T':  The Transpose TRANSR of RFP A is stored.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is REAL array, dimension ( N*(N+1)/2 )
          On entry, the symmetric matrix A in RFP format. RFP format is
          described by TRANSR, UPLO, and N as follows: If TRANSR = 'N'
          then RFP A is (0:N,0:k-1) when N is even; k=N/2. RFP A is
          (0:N-1,0:k) when N is odd; k=N/2. IF TRANSR = 'T' then RFP is
          the transpose of RFP A as defined when
          TRANSR = 'N'. The contents of RFP A are defined by UPLO as
          follows: If UPLO = 'U' the RFP A contains the nt elements of
          upper packed A. If UPLO = 'L' the RFP A contains the elements
          of lower packed A. The LDA of RFP A is (N+1)/2 when TRANSR =
          'T'. When TRANSR is 'N' the LDA is N+1 when N is even and N
          is odd. See the Note below for more details.

          On exit, the symmetric inverse of the original matrix, in the
          same storage format.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the (i,i) element of the factor U or L is
                zero, and the inverse could not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  We first consider Rectangular Full Packed (RFP) Format when N is
  even. We give an example where N = 6.

      AP is Upper             AP is Lower

   00 01 02 03 04 05       00
      11 12 13 14 15       10 11
         22 23 24 25       20 21 22
            33 34 35       30 31 32 33
               44 45       40 41 42 43 44
                  55       50 51 52 53 54 55


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
  the transpose of the first three columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
  the transpose of the last three columns of AP lower.
  This covers the case N even and TRANSR = 'N'.

         RFP A                   RFP A

        03 04 05                33 43 53
        13 14 15                00 44 54
        23 24 25                10 11 55
        33 34 35                20 21 22
        00 44 45                30 31 32
        01 11 55                40 41 42
        02 12 22                50 51 52

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:


           RFP A                   RFP A

     03 13 23 33 00 01 02    33 00 10 20 30 40 50
     04 14 24 34 44 11 12    43 44 11 21 31 41 51
     05 15 25 35 45 55 22    53 54 55 22 32 42 52


  We then consider Rectangular Full Packed (RFP) Format when N is
  odd. We give an example where N = 5.

     AP is Upper                 AP is Lower

   00 01 02 03 04              00
      11 12 13 14              10 11
         22 23 24              20 21 22
            33 34              30 31 32 33
               44              40 41 42 43 44


  Let TRANSR = 'N'. RFP holds AP as follows:
  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
  the transpose of the first two columns of AP upper.
  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
  the transpose of the last two columns of AP lower.
  This covers the case N odd and TRANSR = 'N'.

         RFP A                   RFP A

        02 03 04                00 33 43
        12 13 14                10 11 44
        22 23 24                20 21 22
        00 33 34                30 31 32
        01 11 44                40 41 42

  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
  transpose of RFP A above. One therefore gets:

           RFP A                   RFP A

     02 12 22 00 01             00 10 20 30 40 50
     03 13 23 33 11             33 11 21 31 41 51
     04 14 24 34 44             43 44 22 32 42 52

Definition at line 190 of file spftri.f.

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 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 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
Here is the call graph for this function:
Here is the caller graph for this function: