1 SUBROUTINE ZSYT03( 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 * ZSYT03 computes the residual for a complex symmetric matrix times
23 * its 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 * complex symmetric 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 complex symmetric 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 symmetric
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 * RCOND = 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 *
73 * .. Parameters ..
74 DOUBLE PRECISION ZERO, ONE
75 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
76 COMPLEX*16 CZERO, CONE
77 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
78 $ CONE = ( 1.0D+0, 0.0D+0 ) )
79 * ..
80 * .. Local Scalars ..
81 INTEGER I, J
82 DOUBLE PRECISION AINVNM, ANORM, EPS
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
87 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANSY
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL ZSYMM
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC DBLE
94 * ..
95 * .. Executable Statements ..
96 *
97 * Quick exit if N = 0
98 *
99 IF( N.LE.0 ) THEN
100 RCOND = ONE
101 RESID = ZERO
102 RETURN
103 END IF
104 *
105 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
106 *
107 EPS = DLAMCH( 'Epsilon' )
108 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
109 AINVNM = ZLANSY( '1', UPLO, N, AINV, LDAINV, RWORK )
110 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
111 RCOND = ZERO
112 RESID = ONE / EPS
113 RETURN
114 END IF
115 RCOND = ( ONE / ANORM ) / AINVNM
116 *
117 * Expand AINV into a full matrix and call ZSYMM to multiply
118 * AINV on the left by A (store the result in WORK).
119 *
120 IF( LSAME( UPLO, 'U' ) ) THEN
121 DO 20 J = 1, N
122 DO 10 I = 1, J - 1
123 AINV( J, I ) = AINV( I, J )
124 10 CONTINUE
125 20 CONTINUE
126 ELSE
127 DO 40 J = 1, N
128 DO 30 I = J + 1, N
129 AINV( J, I ) = AINV( I, J )
130 30 CONTINUE
131 40 CONTINUE
132 END IF
133 CALL ZSYMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
134 $ CZERO, WORK, LDWORK )
135 *
136 * Add the identity matrix to WORK .
137 *
138 DO 50 I = 1, N
139 WORK( I, I ) = WORK( I, I ) + CONE
140 50 CONTINUE
141 *
142 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
143 *
144 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
145 *
146 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
147 *
148 RETURN
149 *
150 * End of ZSYT03
151 *
152 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 * ZSYT03 computes the residual for a complex symmetric matrix times
23 * its 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 * complex symmetric 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 complex symmetric 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 symmetric
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 * RCOND = 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 *
73 * .. Parameters ..
74 DOUBLE PRECISION ZERO, ONE
75 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
76 COMPLEX*16 CZERO, CONE
77 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
78 $ CONE = ( 1.0D+0, 0.0D+0 ) )
79 * ..
80 * .. Local Scalars ..
81 INTEGER I, J
82 DOUBLE PRECISION AINVNM, ANORM, EPS
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
87 EXTERNAL LSAME, DLAMCH, ZLANGE, ZLANSY
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL ZSYMM
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC DBLE
94 * ..
95 * .. Executable Statements ..
96 *
97 * Quick exit if N = 0
98 *
99 IF( N.LE.0 ) THEN
100 RCOND = ONE
101 RESID = ZERO
102 RETURN
103 END IF
104 *
105 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
106 *
107 EPS = DLAMCH( 'Epsilon' )
108 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
109 AINVNM = ZLANSY( '1', UPLO, N, AINV, LDAINV, RWORK )
110 IF( ANORM.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
111 RCOND = ZERO
112 RESID = ONE / EPS
113 RETURN
114 END IF
115 RCOND = ( ONE / ANORM ) / AINVNM
116 *
117 * Expand AINV into a full matrix and call ZSYMM to multiply
118 * AINV on the left by A (store the result in WORK).
119 *
120 IF( LSAME( UPLO, 'U' ) ) THEN
121 DO 20 J = 1, N
122 DO 10 I = 1, J - 1
123 AINV( J, I ) = AINV( I, J )
124 10 CONTINUE
125 20 CONTINUE
126 ELSE
127 DO 40 J = 1, N
128 DO 30 I = J + 1, N
129 AINV( J, I ) = AINV( I, J )
130 30 CONTINUE
131 40 CONTINUE
132 END IF
133 CALL ZSYMM( 'Left', UPLO, N, N, -CONE, A, LDA, AINV, LDAINV,
134 $ CZERO, WORK, LDWORK )
135 *
136 * Add the identity matrix to WORK .
137 *
138 DO 50 I = 1, N
139 WORK( I, I ) = WORK( I, I ) + CONE
140 50 CONTINUE
141 *
142 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
143 *
144 RESID = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
145 *
146 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
147 *
148 RETURN
149 *
150 * End of ZSYT03
151 *
152 END