1 SUBROUTINE DTPT01( UPLO, DIAG, N, AP, AINVP, RCOND, WORK, RESID )
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER DIAG, UPLO
9 INTEGER N
10 DOUBLE PRECISION RCOND, RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AINVP( * ), AP( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DTPT01 computes the residual for a triangular matrix A times its
20 * inverse when A is stored in packed format:
21 * RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
22 * where EPS is the machine epsilon.
23 *
24 * Arguments
25 * ==========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the matrix A is upper or lower triangular.
29 * = 'U': Upper triangular
30 * = 'L': Lower triangular
31 *
32 * DIAG (input) CHARACTER*1
33 * Specifies whether or not the matrix A is unit triangular.
34 * = 'N': Non-unit triangular
35 * = 'U': Unit triangular
36 *
37 * N (input) INTEGER
38 * The order of the matrix A. N >= 0.
39 *
40 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
41 * The original upper or lower triangular matrix A, packed
42 * columnwise in a linear array. The j-th column of A is stored
43 * in the array AP as follows:
44 * if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
45 * if UPLO = 'L',
46 * AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
47 *
48 * AINVP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
49 * On entry, the (triangular) inverse of the matrix A, packed
50 * columnwise in a linear array as in AP.
51 * On exit, the contents of AINVP are destroyed.
52 *
53 * RCOND (output) DOUBLE PRECISION
54 * The reciprocal condition number of A, computed as
55 * 1/(norm(A) * norm(AINV)).
56 *
57 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
58 *
59 * RESID (output) DOUBLE PRECISION
60 * norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
61 *
62 * =====================================================================
63 *
64 * .. Parameters ..
65 DOUBLE PRECISION ZERO, ONE
66 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
67 * ..
68 * .. Local Scalars ..
69 LOGICAL UNITD
70 INTEGER J, JC
71 DOUBLE PRECISION AINVNM, ANORM, EPS
72 * ..
73 * .. External Functions ..
74 LOGICAL LSAME
75 DOUBLE PRECISION DLAMCH, DLANTP
76 EXTERNAL LSAME, DLAMCH, DLANTP
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL DTPMV
80 * ..
81 * .. Intrinsic Functions ..
82 INTRINSIC DBLE
83 * ..
84 * .. Executable Statements ..
85 *
86 * Quick exit if N = 0.
87 *
88 IF( N.LE.0 ) THEN
89 RCOND = ONE
90 RESID = ZERO
91 RETURN
92 END IF
93 *
94 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
95 *
96 EPS = DLAMCH( 'Epsilon' )
97 ANORM = DLANTP( '1', UPLO, DIAG, N, AP, WORK )
98 AINVNM = DLANTP( '1', UPLO, DIAG, N, AINVP, WORK )
99 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
100 RCOND = ZERO
101 RESID = ONE / EPS
102 RETURN
103 END IF
104 RCOND = ( ONE / ANORM ) / AINVNM
105 *
106 * Compute A * AINV, overwriting AINV.
107 *
108 UNITD = LSAME( DIAG, 'U' )
109 IF( LSAME( UPLO, 'U' ) ) THEN
110 JC = 1
111 DO 10 J = 1, N
112 IF( UNITD )
113 $ AINVP( JC+J-1 ) = ONE
114 *
115 * Form the j-th column of A*AINV
116 *
117 CALL DTPMV( 'Upper', 'No transpose', DIAG, J, AP,
118 $ AINVP( JC ), 1 )
119 *
120 * Subtract 1 from the diagonal
121 *
122 AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE
123 JC = JC + J
124 10 CONTINUE
125 ELSE
126 JC = 1
127 DO 20 J = 1, N
128 IF( UNITD )
129 $ AINVP( JC ) = ONE
130 *
131 * Form the j-th column of A*AINV
132 *
133 CALL DTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ),
134 $ AINVP( JC ), 1 )
135 *
136 * Subtract 1 from the diagonal
137 *
138 AINVP( JC ) = AINVP( JC ) - ONE
139 JC = JC + N - J + 1
140 20 CONTINUE
141 END IF
142 *
143 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
144 *
145 RESID = DLANTP( '1', UPLO, 'Non-unit', N, AINVP, WORK )
146 *
147 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS
148 *
149 RETURN
150 *
151 * End of DTPT01
152 *
153 END
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER DIAG, UPLO
9 INTEGER N
10 DOUBLE PRECISION RCOND, RESID
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AINVP( * ), AP( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DTPT01 computes the residual for a triangular matrix A times its
20 * inverse when A is stored in packed format:
21 * RESID = norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS ),
22 * where EPS is the machine epsilon.
23 *
24 * Arguments
25 * ==========
26 *
27 * UPLO (input) CHARACTER*1
28 * Specifies whether the matrix A is upper or lower triangular.
29 * = 'U': Upper triangular
30 * = 'L': Lower triangular
31 *
32 * DIAG (input) CHARACTER*1
33 * Specifies whether or not the matrix A is unit triangular.
34 * = 'N': Non-unit triangular
35 * = 'U': Unit triangular
36 *
37 * N (input) INTEGER
38 * The order of the matrix A. N >= 0.
39 *
40 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
41 * The original upper or lower triangular matrix A, packed
42 * columnwise in a linear array. The j-th column of A is stored
43 * in the array AP as follows:
44 * if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
45 * if UPLO = 'L',
46 * AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
47 *
48 * AINVP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
49 * On entry, the (triangular) inverse of the matrix A, packed
50 * columnwise in a linear array as in AP.
51 * On exit, the contents of AINVP are destroyed.
52 *
53 * RCOND (output) DOUBLE PRECISION
54 * The reciprocal condition number of A, computed as
55 * 1/(norm(A) * norm(AINV)).
56 *
57 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
58 *
59 * RESID (output) DOUBLE PRECISION
60 * norm(A*AINV - I) / ( N * norm(A) * norm(AINV) * EPS )
61 *
62 * =====================================================================
63 *
64 * .. Parameters ..
65 DOUBLE PRECISION ZERO, ONE
66 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
67 * ..
68 * .. Local Scalars ..
69 LOGICAL UNITD
70 INTEGER J, JC
71 DOUBLE PRECISION AINVNM, ANORM, EPS
72 * ..
73 * .. External Functions ..
74 LOGICAL LSAME
75 DOUBLE PRECISION DLAMCH, DLANTP
76 EXTERNAL LSAME, DLAMCH, DLANTP
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL DTPMV
80 * ..
81 * .. Intrinsic Functions ..
82 INTRINSIC DBLE
83 * ..
84 * .. Executable Statements ..
85 *
86 * Quick exit if N = 0.
87 *
88 IF( N.LE.0 ) THEN
89 RCOND = ONE
90 RESID = ZERO
91 RETURN
92 END IF
93 *
94 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
95 *
96 EPS = DLAMCH( 'Epsilon' )
97 ANORM = DLANTP( '1', UPLO, DIAG, N, AP, WORK )
98 AINVNM = DLANTP( '1', UPLO, DIAG, N, AINVP, WORK )
99 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
100 RCOND = ZERO
101 RESID = ONE / EPS
102 RETURN
103 END IF
104 RCOND = ( ONE / ANORM ) / AINVNM
105 *
106 * Compute A * AINV, overwriting AINV.
107 *
108 UNITD = LSAME( DIAG, 'U' )
109 IF( LSAME( UPLO, 'U' ) ) THEN
110 JC = 1
111 DO 10 J = 1, N
112 IF( UNITD )
113 $ AINVP( JC+J-1 ) = ONE
114 *
115 * Form the j-th column of A*AINV
116 *
117 CALL DTPMV( 'Upper', 'No transpose', DIAG, J, AP,
118 $ AINVP( JC ), 1 )
119 *
120 * Subtract 1 from the diagonal
121 *
122 AINVP( JC+J-1 ) = AINVP( JC+J-1 ) - ONE
123 JC = JC + J
124 10 CONTINUE
125 ELSE
126 JC = 1
127 DO 20 J = 1, N
128 IF( UNITD )
129 $ AINVP( JC ) = ONE
130 *
131 * Form the j-th column of A*AINV
132 *
133 CALL DTPMV( 'Lower', 'No transpose', DIAG, N-J+1, AP( JC ),
134 $ AINVP( JC ), 1 )
135 *
136 * Subtract 1 from the diagonal
137 *
138 AINVP( JC ) = AINVP( JC ) - ONE
139 JC = JC + N - J + 1
140 20 CONTINUE
141 END IF
142 *
143 * Compute norm(A*AINV - I) / (N * norm(A) * norm(AINV) * EPS)
144 *
145 RESID = DLANTP( '1', UPLO, 'Non-unit', N, AINVP, WORK )
146 *
147 RESID = ( ( RESID*RCOND ) / DBLE( N ) ) / EPS
148 *
149 RETURN
150 *
151 * End of DTPT01
152 *
153 END