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

DTFTRI

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

Purpose: 
 DTFTRI 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 DOUBLE PRECISION array, dimension (0:nt-1);
          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 dtftri.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  DOUBLE PRECISION a( 0: * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  DOUBLE PRECISION one
221  parameter ( one = 1.0d+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, dtrmm, dtrtri
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( 'DTFTRI', -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 dtrtri( 'L', diag, n1, a( 0 ), n, info )
302  IF( info.GT.0 )
303  $ RETURN
304  CALL dtrmm( 'R', 'L', 'N', diag, n2, n1, -one, a( 0 ),
305  $ n, a( n1 ), n )
306  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'L', diag, n1, a( n2 ), n, info )
321  IF( info.GT.0 )
322  $ RETURN
323  CALL dtrmm( 'L', 'L', 'T', diag, n1, n2, -one, a( n2 ),
324  $ n, a( 0 ), n )
325  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'U', diag, n1, a( 0 ), n1, info )
345  IF( info.GT.0 )
346  $ RETURN
347  CALL dtrmm( 'L', 'U', 'N', diag, n1, n2, -one, a( 0 ),
348  $ n1, a( n1*n1 ), n1 )
349  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'U', diag, n1, a( n2*n2 ), n2, info )
363  IF( info.GT.0 )
364  $ RETURN
365  CALL dtrmm( 'R', 'U', 'T', diag, n2, n1, -one,
366  $ a( n2*n2 ), n2, a( 0 ), n2 )
367  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'L', diag, k, a( 1 ), n+1, info )
393  IF( info.GT.0 )
394  $ RETURN
395  CALL dtrmm( 'R', 'L', 'N', diag, k, k, -one, a( 1 ),
396  $ n+1, a( k+1 ), n+1 )
397  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'L', diag, k, a( k+1 ), n+1, info )
412  IF( info.GT.0 )
413  $ RETURN
414  CALL dtrmm( 'L', 'L', 'T', diag, k, k, -one, a( k+1 ),
415  $ n+1, a( 0 ), n+1 )
416  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'U', diag, k, a( k ), k, info )
435  IF( info.GT.0 )
436  $ RETURN
437  CALL dtrmm( 'L', 'U', 'N', diag, k, k, -one, a( k ), k,
438  $ a( k*( k+1 ) ), k )
439  CALL dtrtri( '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 dtrmm( '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 dtrtri( 'U', diag, k, a( k*( k+1 ) ), k, info )
453  IF( info.GT.0 )
454  $ RETURN
455  CALL dtrmm( 'R', 'U', 'T', diag, k, k, -one,
456  $ a( k*( k+1 ) ), k, a( 0 ), k )
457  CALL dtrtri( '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 dtrmm( '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 DTFTRI
471 *
dtrmm
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dtrtri
subroutine dtrtri(UPLO, DIAG, N, A, LDA, INFO)
DTRTRI
Definition: dtrtri.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: