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