LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dtftri()

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.
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 200 of file dtftri.f.

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