LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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)
Definition cblat2.f:3285
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
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
Here is the call graph for this function:
Here is the caller graph for this function: