LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine stftri ( character  TRANSR,
character  UPLO,
character  DIAG,
integer  N,
real, dimension( 0: * )  A,
integer  INFO 
)

STFTRI

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

Purpose: 
 STFTRI computes the inverse of a triangular matrix A stored in RFP
 format.

 This is a Level 3 BLAS version of the algorithm.
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':  A is upper triangular;
          = 'L':  A is lower triangular.
[in]DIAG
          DIAG is CHARACTER*1
          = 'N':  A is non-unit triangular;
          = 'U':  A is unit triangular.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in,out]A
          A is REAL array, dimension (NT);
          NT=N*(N+1)/2. On entry, the triangular factor of a Hermitian
          Positive Definite 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 nt
          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 (triangular) 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, A(i,i) is exactly zero.  The triangular
               matrix is singular and its inverse can not be computed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011
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 203 of file stftri.f.

203 *
204 * -- LAPACK computational routine (version 3.4.0) --
205 * -- LAPACK is a software package provided by Univ. of Tennessee, --
206 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
207 * November 2011
208 *
209 * .. Scalar Arguments ..
210  CHARACTER transr, uplo, diag
211  INTEGER info, n
212 * ..
213 * .. Array Arguments ..
214  REAL a( 0: * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  REAL one
221  parameter ( one = 1.0e+0 )
222 * ..
223 * .. Local Scalars ..
224  LOGICAL lower, nisodd, normaltransr
225  INTEGER n1, n2, k
226 * ..
227 * .. External Functions ..
228  LOGICAL lsame
229  EXTERNAL lsame
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL xerbla, strmm, strtri
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC mod
236 * ..
237 * .. Executable Statements ..
238 *
239 * Test the input parameters.
240 *
241  info = 0
242  normaltransr = lsame( transr, 'N' )
243  lower = lsame( uplo, 'L' )
244  IF( .NOT.normaltransr .AND. .NOT.lsame( transr, 'T' ) ) THEN
245  info = -1
246  ELSE IF( .NOT.lower .AND. .NOT.lsame( uplo, 'U' ) ) THEN
247  info = -2
248  ELSE IF( .NOT.lsame( diag, 'N' ) .AND. .NOT.lsame( diag, 'U' ) )
249  $ THEN
250  info = -3
251  ELSE IF( n.LT.0 ) THEN
252  info = -4
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'STFTRI', -info )
256  RETURN
257  END IF
258 *
259 * Quick return if possible
260 *
261  IF( n.EQ.0 )
262  $ RETURN
263 *
264 * If N is odd, set NISODD = .TRUE.
265 * If N is even, set K = N/2 and NISODD = .FALSE.
266 *
267  IF( mod( n, 2 ).EQ.0 ) THEN
268  k = n / 2
269  nisodd = .false.
270  ELSE
271  nisodd = .true.
272  END IF
273 *
274 * Set N1 and N2 depending on LOWER
275 *
276  IF( lower ) THEN
277  n2 = n / 2
278  n1 = n - n2
279  ELSE
280  n1 = n / 2
281  n2 = n - n1
282  END IF
283 *
284 *
285 * start execution: there are eight cases
286 *
287  IF( nisodd ) THEN
288 *
289 * N is odd
290 *
291  IF( normaltransr ) THEN
292 *
293 * N is odd and TRANSR = 'N'
294 *
295  IF( lower ) THEN
296 *
297 * SRPA for LOWER, NORMAL and N is odd ( a(0:n-1,0:n1-1) )
298 * T1 -> a(0,0), T2 -> a(0,1), S -> a(n1,0)
299 * T1 -> a(0), T2 -> a(n), S -> a(n1)
300 *
301  CALL strtri( 'L', diag, n1, a( 0 ), n, info )
302  IF( info.GT.0 )
303  $ RETURN
304  CALL strmm( 'R', 'L', 'N', diag, n2, n1, -one, a( 0 ),
305  $ n, a( n1 ), n )
306  CALL strtri( 'U', diag, n2, a( n ), n, info )
307  IF( info.GT.0 )
308  $ info = info + n1
309  IF( info.GT.0 )
310  $ RETURN
311  CALL strmm( 'L', 'U', 'T', diag, n2, n1, one, a( n ), n,
312  $ a( n1 ), n )
313 *
314  ELSE
315 *
316 * SRPA for UPPER, NORMAL and N is odd ( a(0:n-1,0:n2-1)
317 * T1 -> a(n1+1,0), T2 -> a(n1,0), S -> a(0,0)
318 * T1 -> a(n2), T2 -> a(n1), S -> a(0)
319 *
320  CALL strtri( 'L', diag, n1, a( n2 ), n, info )
321  IF( info.GT.0 )
322  $ RETURN
323  CALL strmm( 'L', 'L', 'T', diag, n1, n2, -one, a( n2 ),
324  $ n, a( 0 ), n )
325  CALL strtri( 'U', diag, n2, a( n1 ), n, info )
326  IF( info.GT.0 )
327  $ info = info + n1
328  IF( info.GT.0 )
329  $ RETURN
330  CALL strmm( 'R', 'U', 'N', diag, n1, n2, one, a( n1 ),
331  $ n, a( 0 ), n )
332 *
333  END IF
334 *
335  ELSE
336 *
337 * N is odd and TRANSR = 'T'
338 *
339  IF( lower ) THEN
340 *
341 * SRPA for LOWER, TRANSPOSE and N is odd
342 * T1 -> a(0), T2 -> a(1), S -> a(0+n1*n1)
343 *
344  CALL strtri( 'U', diag, n1, a( 0 ), n1, info )
345  IF( info.GT.0 )
346  $ RETURN
347  CALL strmm( 'L', 'U', 'N', diag, n1, n2, -one, a( 0 ),
348  $ n1, a( n1*n1 ), n1 )
349  CALL strtri( 'L', diag, n2, a( 1 ), n1, info )
350  IF( info.GT.0 )
351  $ info = info + n1
352  IF( info.GT.0 )
353  $ RETURN
354  CALL strmm( 'R', 'L', 'T', diag, n1, n2, one, a( 1 ),
355  $ n1, a( n1*n1 ), n1 )
356 *
357  ELSE
358 *
359 * SRPA for UPPER, TRANSPOSE and N is odd
360 * T1 -> a(0+n2*n2), T2 -> a(0+n1*n2), S -> a(0)
361 *
362  CALL strtri( 'U', diag, n1, a( n2*n2 ), n2, info )
363  IF( info.GT.0 )
364  $ RETURN
365  CALL strmm( 'R', 'U', 'T', diag, n2, n1, -one,
366  $ a( n2*n2 ), n2, a( 0 ), n2 )
367  CALL strtri( 'L', diag, n2, a( n1*n2 ), n2, info )
368  IF( info.GT.0 )
369  $ info = info + n1
370  IF( info.GT.0 )
371  $ RETURN
372  CALL strmm( 'L', 'L', 'N', diag, n2, n1, one,
373  $ a( n1*n2 ), n2, a( 0 ), n2 )
374  END IF
375 *
376  END IF
377 *
378  ELSE
379 *
380 * N is even
381 *
382  IF( normaltransr ) THEN
383 *
384 * N is even and TRANSR = 'N'
385 *
386  IF( lower ) THEN
387 *
388 * SRPA for LOWER, NORMAL, and N is even ( a(0:n,0:k-1) )
389 * T1 -> a(1,0), T2 -> a(0,0), S -> a(k+1,0)
390 * T1 -> a(1), T2 -> a(0), S -> a(k+1)
391 *
392  CALL strtri( 'L', diag, k, a( 1 ), n+1, info )
393  IF( info.GT.0 )
394  $ RETURN
395  CALL strmm( 'R', 'L', 'N', diag, k, k, -one, a( 1 ),
396  $ n+1, a( k+1 ), n+1 )
397  CALL strtri( 'U', diag, k, a( 0 ), n+1, info )
398  IF( info.GT.0 )
399  $ info = info + k
400  IF( info.GT.0 )
401  $ RETURN
402  CALL strmm( 'L', 'U', 'T', diag, k, k, one, a( 0 ), n+1,
403  $ a( k+1 ), n+1 )
404 *
405  ELSE
406 *
407 * SRPA for UPPER, NORMAL, and N is even ( a(0:n,0:k-1) )
408 * T1 -> a(k+1,0) , T2 -> a(k,0), S -> a(0,0)
409 * T1 -> a(k+1), T2 -> a(k), S -> a(0)
410 *
411  CALL strtri( 'L', diag, k, a( k+1 ), n+1, info )
412  IF( info.GT.0 )
413  $ RETURN
414  CALL strmm( 'L', 'L', 'T', diag, k, k, -one, a( k+1 ),
415  $ n+1, a( 0 ), n+1 )
416  CALL strtri( 'U', diag, k, a( k ), n+1, info )
417  IF( info.GT.0 )
418  $ info = info + k
419  IF( info.GT.0 )
420  $ RETURN
421  CALL strmm( 'R', 'U', 'N', diag, k, k, one, a( k ), n+1,
422  $ a( 0 ), n+1 )
423  END IF
424  ELSE
425 *
426 * N is even and TRANSR = 'T'
427 *
428  IF( lower ) THEN
429 *
430 * SRPA for LOWER, TRANSPOSE and N is even (see paper)
431 * T1 -> B(0,1), T2 -> B(0,0), S -> B(0,k+1)
432 * T1 -> a(0+k), T2 -> a(0+0), S -> a(0+k*(k+1)); lda=k
433 *
434  CALL strtri( 'U', diag, k, a( k ), k, info )
435  IF( info.GT.0 )
436  $ RETURN
437  CALL strmm( 'L', 'U', 'N', diag, k, k, -one, a( k ), k,
438  $ a( k*( k+1 ) ), k )
439  CALL strtri( 'L', diag, k, a( 0 ), k, info )
440  IF( info.GT.0 )
441  $ info = info + k
442  IF( info.GT.0 )
443  $ RETURN
444  CALL strmm( 'R', 'L', 'T', diag, k, k, one, a( 0 ), k,
445  $ a( k*( k+1 ) ), k )
446  ELSE
447 *
448 * SRPA for UPPER, TRANSPOSE and N is even (see paper)
449 * T1 -> B(0,k+1), T2 -> B(0,k), S -> B(0,0)
450 * T1 -> a(0+k*(k+1)), T2 -> a(0+k*k), S -> a(0+0)); lda=k
451 *
452  CALL strtri( 'U', diag, k, a( k*( k+1 ) ), k, info )
453  IF( info.GT.0 )
454  $ RETURN
455  CALL strmm( 'R', 'U', 'T', diag, k, k, -one,
456  $ a( k*( k+1 ) ), k, a( 0 ), k )
457  CALL strtri( 'L', diag, k, a( k*k ), k, info )
458  IF( info.GT.0 )
459  $ info = info + k
460  IF( info.GT.0 )
461  $ RETURN
462  CALL strmm( 'L', 'L', 'N', diag, k, k, one, a( k*k ), k,
463  $ a( 0 ), k )
464  END IF
465  END IF
466  END IF
467 *
468  RETURN
469 *
470 * End of STFTRI
471 *
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
strmm
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:179
strtri
subroutine strtri(UPLO, DIAG, N, A, LDA, INFO)
STRTRI
Definition: strtri.f:111
lsame
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: