1 SUBROUTINE ZPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
2 $ RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDWORK, N
11 DOUBLE PRECISION RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( * ), AINV( * ), WORK( LDWORK, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZPPT03 computes the residual for a Hermitian packed matrix times its
22 * inverse:
23 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the upper or lower triangular part of the
31 * Hermitian matrix A is stored:
32 * = 'U': Upper triangular
33 * = 'L': Lower triangular
34 *
35 * N (input) INTEGER
36 * The number of rows and columns of the matrix A. N >= 0.
37 *
38 * A (input) COMPLEX*16 array, dimension (N*(N+1)/2)
39 * The original Hermitian matrix A, stored as a packed
40 * triangular matrix.
41 *
42 * AINV (input) COMPLEX*16 array, dimension (N*(N+1)/2)
43 * The (Hermitian) inverse of the matrix A, stored as a packed
44 * triangular matrix.
45 *
46 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N)
47 *
48 * LDWORK (input) INTEGER
49 * The leading dimension of the array WORK. LDWORK >= max(1,N).
50 *
51 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
52 *
53 * RCOND (output) DOUBLE PRECISION
54 * The reciprocal of the condition number of A, computed as
55 * ( 1/norm(A) ) / norm(AINV).
56 *
57 * RESID (output) DOUBLE PRECISION
58 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ZERO, ONE
64 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
65 COMPLEX*16 CZERO, CONE
66 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
67 $ CONE = ( 1.0D+0, 0.0D+0 ) )
68 * ..
69 * .. Local Scalars ..
70 INTEGER I, J, JJ
71 DOUBLE PRECISION AINVNM, ANORM, EPS
72 * ..
73 * .. External Functions ..
74 LOGICAL LSAME
75 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHP
76 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHP
77 * ..
78 * .. Intrinsic Functions ..
79 INTRINSIC DBLE, DCONJG
80 * ..
81 * .. External Subroutines ..
82 EXTERNAL ZCOPY, ZHPMV
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 = ZLANHP( '1', UPLO, N, A, RWORK )
98 AINVNM = ZLANHP( '1', UPLO, N, AINV, RWORK )
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 * UPLO = 'U':
107 * Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
108 * expand it to a full matrix, then multiply by A one column at a
109 * time, moving the result one column to the left.
110 *
111 IF( LSAME( UPLO, 'U' ) ) THEN
112 *
113 * Copy AINV
114 *
115 JJ = 1
116 DO 20 J = 1, N - 1
117 CALL ZCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
118 DO 10 I = 1, J - 1
119 WORK( J, I+1 ) = DCONJG( AINV( JJ+I-1 ) )
120 10 CONTINUE
121 JJ = JJ + J
122 20 CONTINUE
123 JJ = ( ( N-1 )*N ) / 2 + 1
124 DO 30 I = 1, N - 1
125 WORK( N, I+1 ) = DCONJG( AINV( JJ+I-1 ) )
126 30 CONTINUE
127 *
128 * Multiply by A
129 *
130 DO 40 J = 1, N - 1
131 CALL ZHPMV( 'Upper', N, -CONE, A, WORK( 1, J+1 ), 1, CZERO,
132 $ WORK( 1, J ), 1 )
133 40 CONTINUE
134 CALL ZHPMV( 'Upper', N, -CONE, A, AINV( JJ ), 1, CZERO,
135 $ WORK( 1, N ), 1 )
136 *
137 * UPLO = 'L':
138 * Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
139 * and multiply by A, moving each column to the right.
140 *
141 ELSE
142 *
143 * Copy AINV
144 *
145 DO 50 I = 1, N - 1
146 WORK( 1, I ) = DCONJG( AINV( I+1 ) )
147 50 CONTINUE
148 JJ = N + 1
149 DO 70 J = 2, N
150 CALL ZCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
151 DO 60 I = 1, N - J
152 WORK( J, J+I-1 ) = DCONJG( AINV( JJ+I ) )
153 60 CONTINUE
154 JJ = JJ + N - J + 1
155 70 CONTINUE
156 *
157 * Multiply by A
158 *
159 DO 80 J = N, 2, -1
160 CALL ZHPMV( 'Lower', N, -CONE, A, WORK( 1, J-1 ), 1, CZERO,
161 $ WORK( 1, J ), 1 )
162 80 CONTINUE
163 CALL ZHPMV( 'Lower', N, -CONE, A, AINV( 1 ), 1, CZERO,
164 $ WORK( 1, 1 ), 1 )
165 *
166 END IF
167 *
168 * Add the identity matrix to WORK .
169 *
170 DO 90 I = 1, N
171 WORK( I, I ) = WORK( I, I ) + CONE
172 90 CONTINUE
173 *
174 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
175 *
176 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
177 *
178 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
179 *
180 RETURN
181 *
182 * End of ZPPT03
183 *
184 END
2 $ RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDWORK, N
11 DOUBLE PRECISION RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( * ), AINV( * ), WORK( LDWORK, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZPPT03 computes the residual for a Hermitian packed matrix times its
22 * inverse:
23 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the upper or lower triangular part of the
31 * Hermitian matrix A is stored:
32 * = 'U': Upper triangular
33 * = 'L': Lower triangular
34 *
35 * N (input) INTEGER
36 * The number of rows and columns of the matrix A. N >= 0.
37 *
38 * A (input) COMPLEX*16 array, dimension (N*(N+1)/2)
39 * The original Hermitian matrix A, stored as a packed
40 * triangular matrix.
41 *
42 * AINV (input) COMPLEX*16 array, dimension (N*(N+1)/2)
43 * The (Hermitian) inverse of the matrix A, stored as a packed
44 * triangular matrix.
45 *
46 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N)
47 *
48 * LDWORK (input) INTEGER
49 * The leading dimension of the array WORK. LDWORK >= max(1,N).
50 *
51 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
52 *
53 * RCOND (output) DOUBLE PRECISION
54 * The reciprocal of the condition number of A, computed as
55 * ( 1/norm(A) ) / norm(AINV).
56 *
57 * RESID (output) DOUBLE PRECISION
58 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ZERO, ONE
64 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
65 COMPLEX*16 CZERO, CONE
66 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
67 $ CONE = ( 1.0D+0, 0.0D+0 ) )
68 * ..
69 * .. Local Scalars ..
70 INTEGER I, J, JJ
71 DOUBLE PRECISION AINVNM, ANORM, EPS
72 * ..
73 * .. External Functions ..
74 LOGICAL LSAME
75 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHP
76 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHP
77 * ..
78 * .. Intrinsic Functions ..
79 INTRINSIC DBLE, DCONJG
80 * ..
81 * .. External Subroutines ..
82 EXTERNAL ZCOPY, ZHPMV
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 = ZLANHP( '1', UPLO, N, A, RWORK )
98 AINVNM = ZLANHP( '1', UPLO, N, AINV, RWORK )
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 * UPLO = 'U':
107 * Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
108 * expand it to a full matrix, then multiply by A one column at a
109 * time, moving the result one column to the left.
110 *
111 IF( LSAME( UPLO, 'U' ) ) THEN
112 *
113 * Copy AINV
114 *
115 JJ = 1
116 DO 20 J = 1, N - 1
117 CALL ZCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
118 DO 10 I = 1, J - 1
119 WORK( J, I+1 ) = DCONJG( AINV( JJ+I-1 ) )
120 10 CONTINUE
121 JJ = JJ + J
122 20 CONTINUE
123 JJ = ( ( N-1 )*N ) / 2 + 1
124 DO 30 I = 1, N - 1
125 WORK( N, I+1 ) = DCONJG( AINV( JJ+I-1 ) )
126 30 CONTINUE
127 *
128 * Multiply by A
129 *
130 DO 40 J = 1, N - 1
131 CALL ZHPMV( 'Upper', N, -CONE, A, WORK( 1, J+1 ), 1, CZERO,
132 $ WORK( 1, J ), 1 )
133 40 CONTINUE
134 CALL ZHPMV( 'Upper', N, -CONE, A, AINV( JJ ), 1, CZERO,
135 $ WORK( 1, N ), 1 )
136 *
137 * UPLO = 'L':
138 * Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
139 * and multiply by A, moving each column to the right.
140 *
141 ELSE
142 *
143 * Copy AINV
144 *
145 DO 50 I = 1, N - 1
146 WORK( 1, I ) = DCONJG( AINV( I+1 ) )
147 50 CONTINUE
148 JJ = N + 1
149 DO 70 J = 2, N
150 CALL ZCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
151 DO 60 I = 1, N - J
152 WORK( J, J+I-1 ) = DCONJG( AINV( JJ+I ) )
153 60 CONTINUE
154 JJ = JJ + N - J + 1
155 70 CONTINUE
156 *
157 * Multiply by A
158 *
159 DO 80 J = N, 2, -1
160 CALL ZHPMV( 'Lower', N, -CONE, A, WORK( 1, J-1 ), 1, CZERO,
161 $ WORK( 1, J ), 1 )
162 80 CONTINUE
163 CALL ZHPMV( 'Lower', N, -CONE, A, AINV( 1 ), 1, CZERO,
164 $ WORK( 1, 1 ), 1 )
165 *
166 END IF
167 *
168 * Add the identity matrix to WORK .
169 *
170 DO 90 I = 1, N
171 WORK( I, I ) = WORK( I, I ) + CONE
172 90 CONTINUE
173 *
174 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
175 *
176 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
177 *
178 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
179 *
180 RETURN
181 *
182 * End of ZPPT03
183 *
184 END