1 SUBROUTINE ZSYT02( UPLO, N, NRHS, A, LDA, X, LDX, B, LDB, RWORK,
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 LDA, LDB, LDX, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZSYT02 computes the residual for a solution to a complex symmetric
22 * system of linear equations A*x = b:
23 *
24 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
25 *
26 * where EPS is the machine epsilon.
27 *
28 * Arguments
29 * =========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the upper or lower triangular part of the
33 * symmetric matrix A is stored:
34 * = 'U': Upper triangular
35 * = 'L': Lower triangular
36 *
37 * N (input) INTEGER
38 * The number of rows and columns of the matrix A. N >= 0.
39 *
40 * NRHS (input) INTEGER
41 * The number of columns of B, the matrix of right hand sides.
42 * NRHS >= 0.
43 *
44 * A (input) COMPLEX*16 array, dimension (LDA,N)
45 * The original complex symmetric matrix A.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A. LDA >= max(1,N)
49 *
50 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
51 * The computed solution vectors for the system of linear
52 * equations.
53 *
54 * LDX (input) INTEGER
55 * The leading dimension of the array X. LDX >= max(1,N).
56 *
57 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
58 * On entry, the right hand side vectors for the system of
59 * linear equations.
60 * On exit, B is overwritten with the difference B - A*X.
61 *
62 * LDB (input) INTEGER
63 * The leading dimension of the array B. LDB >= max(1,N).
64 *
65 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
66 *
67 * RESID (output) DOUBLE PRECISION
68 * The maximum over the number of right hand sides of
69 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
70 *
71 * =====================================================================
72 *
73 * .. Parameters ..
74 DOUBLE PRECISION ZERO, ONE
75 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
76 COMPLEX*16 CONE
77 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
78 * ..
79 * .. Local Scalars ..
80 INTEGER J
81 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
82 * ..
83 * .. External Functions ..
84 DOUBLE PRECISION DLAMCH, DZASUM, ZLANSY
85 EXTERNAL DLAMCH, DZASUM, ZLANSY
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL ZSYMM
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC MAX
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0 or NRHS = 0
96 *
97 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 * Exit with RESID = 1/EPS if ANORM = 0.
103 *
104 EPS = DLAMCH( 'Epsilon' )
105 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
106 IF( ANORM.LE.ZERO ) THEN
107 RESID = ONE / EPS
108 RETURN
109 END IF
110 *
111 * Compute B - A*X (or B - A'*X ) and store in B .
112 *
113 CALL ZSYMM( 'Left', UPLO, N, NRHS, -CONE, A, LDA, X, LDX, CONE, B,
114 $ LDB )
115 *
116 * Compute the maximum over the number of right hand sides of
117 * norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) .
118 *
119 RESID = ZERO
120 DO 10 J = 1, NRHS
121 BNORM = DZASUM( N, B( 1, J ), 1 )
122 XNORM = DZASUM( N, X( 1, J ), 1 )
123 IF( XNORM.LE.ZERO ) THEN
124 RESID = ONE / EPS
125 ELSE
126 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
127 END IF
128 10 CONTINUE
129 *
130 RETURN
131 *
132 * End of ZSYT02
133 *
134 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 LDA, LDB, LDX, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZSYT02 computes the residual for a solution to a complex symmetric
22 * system of linear equations A*x = b:
23 *
24 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
25 *
26 * where EPS is the machine epsilon.
27 *
28 * Arguments
29 * =========
30 *
31 * UPLO (input) CHARACTER*1
32 * Specifies whether the upper or lower triangular part of the
33 * symmetric matrix A is stored:
34 * = 'U': Upper triangular
35 * = 'L': Lower triangular
36 *
37 * N (input) INTEGER
38 * The number of rows and columns of the matrix A. N >= 0.
39 *
40 * NRHS (input) INTEGER
41 * The number of columns of B, the matrix of right hand sides.
42 * NRHS >= 0.
43 *
44 * A (input) COMPLEX*16 array, dimension (LDA,N)
45 * The original complex symmetric matrix A.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A. LDA >= max(1,N)
49 *
50 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
51 * The computed solution vectors for the system of linear
52 * equations.
53 *
54 * LDX (input) INTEGER
55 * The leading dimension of the array X. LDX >= max(1,N).
56 *
57 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
58 * On entry, the right hand side vectors for the system of
59 * linear equations.
60 * On exit, B is overwritten with the difference B - A*X.
61 *
62 * LDB (input) INTEGER
63 * The leading dimension of the array B. LDB >= max(1,N).
64 *
65 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
66 *
67 * RESID (output) DOUBLE PRECISION
68 * The maximum over the number of right hand sides of
69 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
70 *
71 * =====================================================================
72 *
73 * .. Parameters ..
74 DOUBLE PRECISION ZERO, ONE
75 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
76 COMPLEX*16 CONE
77 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
78 * ..
79 * .. Local Scalars ..
80 INTEGER J
81 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
82 * ..
83 * .. External Functions ..
84 DOUBLE PRECISION DLAMCH, DZASUM, ZLANSY
85 EXTERNAL DLAMCH, DZASUM, ZLANSY
86 * ..
87 * .. External Subroutines ..
88 EXTERNAL ZSYMM
89 * ..
90 * .. Intrinsic Functions ..
91 INTRINSIC MAX
92 * ..
93 * .. Executable Statements ..
94 *
95 * Quick exit if N = 0 or NRHS = 0
96 *
97 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
98 RESID = ZERO
99 RETURN
100 END IF
101 *
102 * Exit with RESID = 1/EPS if ANORM = 0.
103 *
104 EPS = DLAMCH( 'Epsilon' )
105 ANORM = ZLANSY( '1', UPLO, N, A, LDA, RWORK )
106 IF( ANORM.LE.ZERO ) THEN
107 RESID = ONE / EPS
108 RETURN
109 END IF
110 *
111 * Compute B - A*X (or B - A'*X ) and store in B .
112 *
113 CALL ZSYMM( 'Left', UPLO, N, NRHS, -CONE, A, LDA, X, LDX, CONE, B,
114 $ LDB )
115 *
116 * Compute the maximum over the number of right hand sides of
117 * norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) .
118 *
119 RESID = ZERO
120 DO 10 J = 1, NRHS
121 BNORM = DZASUM( N, B( 1, J ), 1 )
122 XNORM = DZASUM( N, X( 1, J ), 1 )
123 IF( XNORM.LE.ZERO ) THEN
124 RESID = ONE / EPS
125 ELSE
126 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
127 END IF
128 10 CONTINUE
129 *
130 RETURN
131 *
132 * End of ZSYT02
133 *
134 END