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