LAPACK  3.8.0
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.
Date
December 2016
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 193 of file spftri.f.

193 *
194 * -- LAPACK computational routine (version 3.7.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 * December 2016
198 *
199 * .. Scalar Arguments ..
200  CHARACTER transr, uplo
201  INTEGER info, n
202 * .. Array Arguments ..
203  REAL a( 0: * )
204 * ..
205 *
206 * =====================================================================
207 *
208 * .. Parameters ..
209  REAL one
210  parameter( one = 1.0e+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, stftri, slauum, strmm, ssyrk
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( 'SPFTRI', -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 stftri( 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 slauum( 'L', n1, a( 0 ), n, info )
294  CALL ssyrk( 'L', 'T', n1, n2, one, a( n1 ), n, one,
295  $ a( 0 ), n )
296  CALL strmm( 'L', 'U', 'N', 'N', n2, n1, one, a( n ), n,
297  $ a( n1 ), n )
298  CALL slauum( '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 slauum( 'L', n1, a( n2 ), n, info )
307  CALL ssyrk( 'L', 'N', n1, n2, one, a( 0 ), n, one,
308  $ a( n2 ), n )
309  CALL strmm( 'R', 'U', 'T', 'N', n1, n2, one, a( n1 ), n,
310  $ a( 0 ), n )
311  CALL slauum( '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 slauum( 'U', n1, a( 0 ), n1, info )
325  CALL ssyrk( 'U', 'N', n1, n2, one, a( n1*n1 ), n1, one,
326  $ a( 0 ), n1 )
327  CALL strmm( 'R', 'L', 'N', 'N', n1, n2, one, a( 1 ), n1,
328  $ a( n1*n1 ), n1 )
329  CALL slauum( '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 slauum( 'U', n1, a( n2*n2 ), n2, info )
337  CALL ssyrk( 'U', 'T', n1, n2, one, a( 0 ), n2, one,
338  $ a( n2*n2 ), n2 )
339  CALL strmm( 'L', 'L', 'T', 'N', n2, n1, one, a( n1*n2 ),
340  $ n2, a( 0 ), n2 )
341  CALL slauum( '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 slauum( 'L', k, a( 1 ), n+1, info )
362  CALL ssyrk( 'L', 'T', k, k, one, a( k+1 ), n+1, one,
363  $ a( 1 ), n+1 )
364  CALL strmm( 'L', 'U', 'N', 'N', k, k, one, a( 0 ), n+1,
365  $ a( k+1 ), n+1 )
366  CALL slauum( '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 slauum( 'L', k, a( k+1 ), n+1, info )
375  CALL ssyrk( 'L', 'N', k, k, one, a( 0 ), n+1, one,
376  $ a( k+1 ), n+1 )
377  CALL strmm( 'R', 'U', 'T', 'N', k, k, one, a( k ), n+1,
378  $ a( 0 ), n+1 )
379  CALL slauum( '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 slauum( 'U', k, a( k ), k, info )
394  CALL ssyrk( 'U', 'N', k, k, one, a( k*( k+1 ) ), k, one,
395  $ a( k ), k )
396  CALL strmm( 'R', 'L', 'N', 'N', k, k, one, a( 0 ), k,
397  $ a( k*( k+1 ) ), k )
398  CALL slauum( '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 slauum( 'U', k, a( k*( k+1 ) ), k, info )
407  CALL ssyrk( 'U', 'T', k, k, one, a( 0 ), k, one,
408  $ a( k*( k+1 ) ), k )
409  CALL strmm( 'L', 'L', 'T', 'N', k, k, one, a( k*k ), k,
410  $ a( 0 ), k )
411  CALL slauum( '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 SPFTRI
422 *
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
Definition: ssyrk.f:171
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:104
subroutine stftri(TRANSR, UPLO, DIAG, N, A, INFO)
STFTRI
Definition: stftri.f:203
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
Here is the call graph for this function:
Here is the caller graph for this function: