1 SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, UPLO
10 INTEGER INFO, LDA, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DTRTRI computes the inverse of a real upper or lower triangular
20 * matrix A.
21 *
22 * This is the Level 3 BLAS version of the algorithm.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * = 'U': A is upper triangular;
29 * = 'L': A is lower triangular.
30 *
31 * DIAG (input) CHARACTER*1
32 * = 'N': A is non-unit triangular;
33 * = 'U': A is unit triangular.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
39 * On entry, the triangular matrix A. If UPLO = 'U', the
40 * leading N-by-N upper triangular part of the array A contains
41 * the upper triangular matrix, and the strictly lower
42 * triangular part of A is not referenced. If UPLO = 'L', the
43 * leading N-by-N lower triangular part of the array A contains
44 * the lower triangular matrix, and the strictly upper
45 * triangular part of A is not referenced. If DIAG = 'U', the
46 * diagonal elements of A are also not referenced and are
47 * assumed to be 1.
48 * On exit, the (triangular) inverse of the original matrix, in
49 * the same storage format.
50 *
51 * LDA (input) INTEGER
52 * The leading dimension of the array A. LDA >= max(1,N).
53 *
54 * INFO (output) INTEGER
55 * = 0: successful exit
56 * < 0: if INFO = -i, the i-th argument had an illegal value
57 * > 0: if INFO = i, A(i,i) is exactly zero. The triangular
58 * matrix is singular and its inverse can not be computed.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE, ZERO
64 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 LOGICAL NOUNIT, UPPER
68 INTEGER J, JB, NB, NN
69 * ..
70 * .. External Functions ..
71 LOGICAL LSAME
72 INTEGER ILAENV
73 EXTERNAL LSAME, ILAENV
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA
77 * ..
78 * .. Intrinsic Functions ..
79 INTRINSIC MAX, MIN
80 * ..
81 * .. Executable Statements ..
82 *
83 * Test the input parameters.
84 *
85 INFO = 0
86 UPPER = LSAME( UPLO, 'U' )
87 NOUNIT = LSAME( DIAG, 'N' )
88 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89 INFO = -1
90 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
91 INFO = -2
92 ELSE IF( N.LT.0 ) THEN
93 INFO = -3
94 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
95 INFO = -5
96 END IF
97 IF( INFO.NE.0 ) THEN
98 CALL XERBLA( 'DTRTRI', -INFO )
99 RETURN
100 END IF
101 *
102 * Quick return if possible
103 *
104 IF( N.EQ.0 )
105 $ RETURN
106 *
107 * Check for singularity if non-unit.
108 *
109 IF( NOUNIT ) THEN
110 DO 10 INFO = 1, N
111 IF( A( INFO, INFO ).EQ.ZERO )
112 $ RETURN
113 10 CONTINUE
114 INFO = 0
115 END IF
116 *
117 * Determine the block size for this environment.
118 *
119 NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 )
120 IF( NB.LE.1 .OR. NB.GE.N ) THEN
121 *
122 * Use unblocked code
123 *
124 CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
125 ELSE
126 *
127 * Use blocked code
128 *
129 IF( UPPER ) THEN
130 *
131 * Compute inverse of upper triangular matrix
132 *
133 DO 20 J = 1, N, NB
134 JB = MIN( NB, N-J+1 )
135 *
136 * Compute rows 1:j-1 of current block column
137 *
138 CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1,
139 $ JB, ONE, A, LDA, A( 1, J ), LDA )
140 CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1,
141 $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA )
142 *
143 * Compute inverse of current diagonal block
144 *
145 CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO )
146 20 CONTINUE
147 ELSE
148 *
149 * Compute inverse of lower triangular matrix
150 *
151 NN = ( ( N-1 ) / NB )*NB + 1
152 DO 30 J = NN, 1, -NB
153 JB = MIN( NB, N-J+1 )
154 IF( J+JB.LE.N ) THEN
155 *
156 * Compute rows j+jb:n of current block column
157 *
158 CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG,
159 $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA,
160 $ A( J+JB, J ), LDA )
161 CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG,
162 $ N-J-JB+1, JB, -ONE, A( J, J ), LDA,
163 $ A( J+JB, J ), LDA )
164 END IF
165 *
166 * Compute inverse of current diagonal block
167 *
168 CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO )
169 30 CONTINUE
170 END IF
171 END IF
172 *
173 RETURN
174 *
175 * End of DTRTRI
176 *
177 END
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER DIAG, UPLO
10 INTEGER INFO, LDA, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DTRTRI computes the inverse of a real upper or lower triangular
20 * matrix A.
21 *
22 * This is the Level 3 BLAS version of the algorithm.
23 *
24 * Arguments
25 * =========
26 *
27 * UPLO (input) CHARACTER*1
28 * = 'U': A is upper triangular;
29 * = 'L': A is lower triangular.
30 *
31 * DIAG (input) CHARACTER*1
32 * = 'N': A is non-unit triangular;
33 * = 'U': A is unit triangular.
34 *
35 * N (input) INTEGER
36 * The order of the matrix A. N >= 0.
37 *
38 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
39 * On entry, the triangular matrix A. If UPLO = 'U', the
40 * leading N-by-N upper triangular part of the array A contains
41 * the upper triangular matrix, and the strictly lower
42 * triangular part of A is not referenced. If UPLO = 'L', the
43 * leading N-by-N lower triangular part of the array A contains
44 * the lower triangular matrix, and the strictly upper
45 * triangular part of A is not referenced. If DIAG = 'U', the
46 * diagonal elements of A are also not referenced and are
47 * assumed to be 1.
48 * On exit, the (triangular) inverse of the original matrix, in
49 * the same storage format.
50 *
51 * LDA (input) INTEGER
52 * The leading dimension of the array A. LDA >= max(1,N).
53 *
54 * INFO (output) INTEGER
55 * = 0: successful exit
56 * < 0: if INFO = -i, the i-th argument had an illegal value
57 * > 0: if INFO = i, A(i,i) is exactly zero. The triangular
58 * matrix is singular and its inverse can not be computed.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ONE, ZERO
64 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 LOGICAL NOUNIT, UPPER
68 INTEGER J, JB, NB, NN
69 * ..
70 * .. External Functions ..
71 LOGICAL LSAME
72 INTEGER ILAENV
73 EXTERNAL LSAME, ILAENV
74 * ..
75 * .. External Subroutines ..
76 EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA
77 * ..
78 * .. Intrinsic Functions ..
79 INTRINSIC MAX, MIN
80 * ..
81 * .. Executable Statements ..
82 *
83 * Test the input parameters.
84 *
85 INFO = 0
86 UPPER = LSAME( UPLO, 'U' )
87 NOUNIT = LSAME( DIAG, 'N' )
88 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
89 INFO = -1
90 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
91 INFO = -2
92 ELSE IF( N.LT.0 ) THEN
93 INFO = -3
94 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
95 INFO = -5
96 END IF
97 IF( INFO.NE.0 ) THEN
98 CALL XERBLA( 'DTRTRI', -INFO )
99 RETURN
100 END IF
101 *
102 * Quick return if possible
103 *
104 IF( N.EQ.0 )
105 $ RETURN
106 *
107 * Check for singularity if non-unit.
108 *
109 IF( NOUNIT ) THEN
110 DO 10 INFO = 1, N
111 IF( A( INFO, INFO ).EQ.ZERO )
112 $ RETURN
113 10 CONTINUE
114 INFO = 0
115 END IF
116 *
117 * Determine the block size for this environment.
118 *
119 NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 )
120 IF( NB.LE.1 .OR. NB.GE.N ) THEN
121 *
122 * Use unblocked code
123 *
124 CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
125 ELSE
126 *
127 * Use blocked code
128 *
129 IF( UPPER ) THEN
130 *
131 * Compute inverse of upper triangular matrix
132 *
133 DO 20 J = 1, N, NB
134 JB = MIN( NB, N-J+1 )
135 *
136 * Compute rows 1:j-1 of current block column
137 *
138 CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1,
139 $ JB, ONE, A, LDA, A( 1, J ), LDA )
140 CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1,
141 $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA )
142 *
143 * Compute inverse of current diagonal block
144 *
145 CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO )
146 20 CONTINUE
147 ELSE
148 *
149 * Compute inverse of lower triangular matrix
150 *
151 NN = ( ( N-1 ) / NB )*NB + 1
152 DO 30 J = NN, 1, -NB
153 JB = MIN( NB, N-J+1 )
154 IF( J+JB.LE.N ) THEN
155 *
156 * Compute rows j+jb:n of current block column
157 *
158 CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG,
159 $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA,
160 $ A( J+JB, J ), LDA )
161 CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG,
162 $ N-J-JB+1, JB, -ONE, A( J, J ), LDA,
163 $ A( J+JB, J ), LDA )
164 END IF
165 *
166 * Compute inverse of current diagonal block
167 *
168 CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO )
169 30 CONTINUE
170 END IF
171 END IF
172 *
173 RETURN
174 *
175 * End of DTRTRI
176 *
177 END