1 SUBROUTINE ZPOT03( UPLO, N, A, LDA, AINV, LDAINV, WORK, LDWORK,
2 $ RWORK, RCOND, 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 LDA, LDAINV, LDWORK, N
11 DOUBLE PRECISION RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), AINV( LDAINV, * ),
16 $ WORK( LDWORK, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZPOT03 computes the residual for a Hermitian matrix times its
23 * inverse:
24 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
25 * where EPS is the machine epsilon.
26 *
27 * Arguments
28 * ==========
29 *
30 * UPLO (input) CHARACTER*1
31 * Specifies whether the upper or lower triangular part of the
32 * Hermitian matrix A is stored:
33 * = 'U': Upper triangular
34 * = 'L': Lower triangular
35 *
36 * N (input) INTEGER
37 * The number of rows and columns of the matrix A. N >= 0.
38 *
39 * A (input) COMPLEX*16 array, dimension (LDA,N)
40 * The original Hermitian matrix A.
41 *
42 * LDA (input) INTEGER
43 * The leading dimension of the array A. LDA >= max(1,N)
44 *
45 * AINV (input/output) COMPLEX*16 array, dimension (LDAINV,N)
46 * On entry, the inverse of the matrix A, stored as a Hermitian
47 * matrix in the same format as A.
48 * In this version, AINV is expanded into a full matrix and
49 * multiplied by A, so the opposing triangle of AINV will be
50 * changed; i.e., if the upper triangular part of AINV is
51 * stored, the lower triangular part will be used as work space.
52 *
53 * LDAINV (input) INTEGER
54 * The leading dimension of the array AINV. LDAINV >= max(1,N).
55 *
56 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N)
57 *
58 * LDWORK (input) INTEGER
59 * The leading dimension of the array WORK. LDWORK >= max(1,N).
60 *
61 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
62 *
63 * RCOND (output) DOUBLE PRECISION
64 * The reciprocal of the condition number of A, computed as
65 * ( 1/norm(A) ) / norm(AINV).
66 *
67 * RESID (output) DOUBLE PRECISION
68 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 DOUBLE PRECISION ZERO, ONE
74 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75 COMPLEX*16 CZERO, CONE
76 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
77 $ CONE = ( 1.0D+0, 0.0D+0 ) )
78 * ..
79 * .. Local Scalars ..
80 INTEGER I, J
81 DOUBLE PRECISION AINVNM, ANORM, EPS
82 * ..
83 * .. External Functions ..
84 LOGICAL LSAME
85 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
86 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHE
87 * ..
88 * .. External Subroutines ..
89 EXTERNAL ZHEMM
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC DBLE, DCONJG
93 * ..
94 * .. Executable Statements ..
95 *
96 * Quick exit if N = 0.
97 *
98 IF( N.LE.0 ) THEN
99 RCOND = ONE
100 RESID = ZERO
101 RETURN
102 END IF
103 *
104 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
105 *
106 EPS = DLAMCH( 'Epsilon' )
107 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
108 AINVNM = ZLANHE( '1', UPLO, N, AINV, LDAINV, RWORK )
109 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
110 RCOND = ZERO
111 RESID = ONE / EPS
112 RETURN
113 END IF
114 RCOND = ( ONE / ANORM ) / AINVNM
115 *
116 * Expand AINV into a full matrix and call ZHEMM to multiply
117 * AINV on the left by A.
118 *
119 IF( LSAME( UPLO, 'U' ) ) THEN
120 DO 20 J = 1, N
121 DO 10 I = 1, J - 1
122 AINV( J, I ) = DCONJG( AINV( I, J ) )
123 10 CONTINUE
124 20 CONTINUE
125 ELSE
126 DO 40 J = 1, N
127 DO 30 I = J + 1, N
128 AINV( J, I ) = DCONJG( AINV( I, J ) )
129 30 CONTINUE
130 40 CONTINUE
131 END IF
132 CALL ZHEMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
133 $ CZERO, WORK, LDWORK )
134 *
135 * Add the identity matrix to WORK .
136 *
137 DO 50 I = 1, N
138 WORK( I, I ) = WORK( I, I ) + CONE
139 50 CONTINUE
140 *
141 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
142 *
143 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
144 *
145 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
146 *
147 RETURN
148 *
149 * End of ZPOT03
150 *
151 END
2 $ RWORK, RCOND, 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 LDA, LDAINV, LDWORK, N
11 DOUBLE PRECISION RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), AINV( LDAINV, * ),
16 $ WORK( LDWORK, * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * ZPOT03 computes the residual for a Hermitian matrix times its
23 * inverse:
24 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
25 * where EPS is the machine epsilon.
26 *
27 * Arguments
28 * ==========
29 *
30 * UPLO (input) CHARACTER*1
31 * Specifies whether the upper or lower triangular part of the
32 * Hermitian matrix A is stored:
33 * = 'U': Upper triangular
34 * = 'L': Lower triangular
35 *
36 * N (input) INTEGER
37 * The number of rows and columns of the matrix A. N >= 0.
38 *
39 * A (input) COMPLEX*16 array, dimension (LDA,N)
40 * The original Hermitian matrix A.
41 *
42 * LDA (input) INTEGER
43 * The leading dimension of the array A. LDA >= max(1,N)
44 *
45 * AINV (input/output) COMPLEX*16 array, dimension (LDAINV,N)
46 * On entry, the inverse of the matrix A, stored as a Hermitian
47 * matrix in the same format as A.
48 * In this version, AINV is expanded into a full matrix and
49 * multiplied by A, so the opposing triangle of AINV will be
50 * changed; i.e., if the upper triangular part of AINV is
51 * stored, the lower triangular part will be used as work space.
52 *
53 * LDAINV (input) INTEGER
54 * The leading dimension of the array AINV. LDAINV >= max(1,N).
55 *
56 * WORK (workspace) COMPLEX*16 array, dimension (LDWORK,N)
57 *
58 * LDWORK (input) INTEGER
59 * The leading dimension of the array WORK. LDWORK >= max(1,N).
60 *
61 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
62 *
63 * RCOND (output) DOUBLE PRECISION
64 * The reciprocal of the condition number of A, computed as
65 * ( 1/norm(A) ) / norm(AINV).
66 *
67 * RESID (output) DOUBLE PRECISION
68 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
69 *
70 * =====================================================================
71 *
72 * .. Parameters ..
73 DOUBLE PRECISION ZERO, ONE
74 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
75 COMPLEX*16 CZERO, CONE
76 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
77 $ CONE = ( 1.0D+0, 0.0D+0 ) )
78 * ..
79 * .. Local Scalars ..
80 INTEGER I, J
81 DOUBLE PRECISION AINVNM, ANORM, EPS
82 * ..
83 * .. External Functions ..
84 LOGICAL LSAME
85 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
86 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANHE
87 * ..
88 * .. External Subroutines ..
89 EXTERNAL ZHEMM
90 * ..
91 * .. Intrinsic Functions ..
92 INTRINSIC DBLE, DCONJG
93 * ..
94 * .. Executable Statements ..
95 *
96 * Quick exit if N = 0.
97 *
98 IF( N.LE.0 ) THEN
99 RCOND = ONE
100 RESID = ZERO
101 RETURN
102 END IF
103 *
104 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
105 *
106 EPS = DLAMCH( 'Epsilon' )
107 ANORM = ZLANHE( '1', UPLO, N, A, LDA, RWORK )
108 AINVNM = ZLANHE( '1', UPLO, N, AINV, LDAINV, RWORK )
109 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
110 RCOND = ZERO
111 RESID = ONE / EPS
112 RETURN
113 END IF
114 RCOND = ( ONE / ANORM ) / AINVNM
115 *
116 * Expand AINV into a full matrix and call ZHEMM to multiply
117 * AINV on the left by A.
118 *
119 IF( LSAME( UPLO, 'U' ) ) THEN
120 DO 20 J = 1, N
121 DO 10 I = 1, J - 1
122 AINV( J, I ) = DCONJG( AINV( I, J ) )
123 10 CONTINUE
124 20 CONTINUE
125 ELSE
126 DO 40 J = 1, N
127 DO 30 I = J + 1, N
128 AINV( J, I ) = DCONJG( AINV( I, J ) )
129 30 CONTINUE
130 40 CONTINUE
131 END IF
132 CALL ZHEMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
133 $ CZERO, WORK, LDWORK )
134 *
135 * Add the identity matrix to WORK .
136 *
137 DO 50 I = 1, N
138 WORK( I, I ) = WORK( I, I ) + CONE
139 50 CONTINUE
140 *
141 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
142 *
143 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
144 *
145 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
146 *
147 RETURN
148 *
149 * End of ZPOT03
150 *
151 END