1 SUBROUTINE ZTPT02( UPLO, TRANS, DIAG, N, NRHS, AP, X, LDX, B, LDB,
2 $ WORK, RWORK, 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 DIAG, TRANS, UPLO
10 INTEGER LDB, LDX, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZTPT02 computes the residual for the computed solution to a
22 * triangular system of linear equations A*x = b, A**T *x = b, or
23 * A**H *x = b, when the triangular matrix A is stored in packed format.
24 * Here A**T denotes the transpose of A, A**H denotes the conjugate
25 * transpose of A, and x and b are N by NRHS matrices. The test ratio
26 * is the maximum over the number of right hand sides of
27 * the maximum over the number of right hand sides of
28 * norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
29 * where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the matrix A is upper or lower triangular.
36 * = 'U': Upper triangular
37 * = 'L': Lower triangular
38 *
39 * TRANS (input) CHARACTER*1
40 * Specifies the operation applied to A.
41 * = 'N': A *x = b (No transpose)
42 * = 'T': A**T *x = b (Transpose)
43 * = 'C': A**H *x = b (Conjugate transpose)
44 *
45 * DIAG (input) CHARACTER*1
46 * Specifies whether or not the matrix A is unit triangular.
47 * = 'N': Non-unit triangular
48 * = 'U': Unit triangular
49 *
50 * N (input) INTEGER
51 * The order of the matrix A. N >= 0.
52 *
53 * NRHS (input) INTEGER
54 * The number of right hand sides, i.e., the number of columns
55 * of the matrices X and B. NRHS >= 0.
56 *
57 * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
58 * The upper or lower triangular matrix A, packed columnwise in
59 * a linear array. The j-th column of A is stored in the array
60 * AP as follows:
61 * if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
62 * if UPLO = 'L',
63 * AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
64 *
65 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
66 * The computed solution vectors for the system of linear
67 * equations.
68 *
69 * LDX (input) INTEGER
70 * The leading dimension of the array X. LDX >= max(1,N).
71 *
72 * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
73 * The right hand side vectors for the system of linear
74 * equations.
75 *
76 * LDB (input) INTEGER
77 * The leading dimension of the array B. LDB >= max(1,N).
78 *
79 * WORK (workspace) COMPLEX*16 array, dimension (N)
80 *
81 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
82 *
83 * RESID (output) DOUBLE PRECISION
84 * The maximum over the number of right hand sides of
85 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
86 *
87 * =====================================================================
88 *
89 * .. Parameters ..
90 DOUBLE PRECISION ZERO, ONE
91 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
92 * ..
93 * .. Local Scalars ..
94 INTEGER J
95 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
96 * ..
97 * .. External Functions ..
98 LOGICAL LSAME
99 DOUBLE PRECISION DLAMCH, DZASUM, ZLANTP
100 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANTP
101 * ..
102 * .. External Subroutines ..
103 EXTERNAL ZAXPY, ZCOPY, ZTPMV
104 * ..
105 * .. Intrinsic Functions ..
106 INTRINSIC DCMPLX, MAX
107 * ..
108 * .. Executable Statements ..
109 *
110 * Quick exit if N = 0 or NRHS = 0
111 *
112 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
113 RESID = ZERO
114 RETURN
115 END IF
116 *
117 * Compute the 1-norm of A or A**H.
118 *
119 IF( LSAME( TRANS, 'N' ) ) THEN
120 ANORM = ZLANTP( '1', UPLO, DIAG, N, AP, RWORK )
121 ELSE
122 ANORM = ZLANTP( 'I', UPLO, DIAG, N, AP, RWORK )
123 END IF
124 *
125 * Exit with RESID = 1/EPS if ANORM = 0.
126 *
127 EPS = DLAMCH( 'Epsilon' )
128 IF( ANORM.LE.ZERO ) THEN
129 RESID = ONE / EPS
130 RETURN
131 END IF
132 *
133 * Compute the maximum over the number of right hand sides of
134 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
135 *
136 RESID = ZERO
137 DO 10 J = 1, NRHS
138 CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
139 CALL ZTPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 )
140 CALL ZAXPY( N, DCMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
141 BNORM = DZASUM( N, WORK, 1 )
142 XNORM = DZASUM( N, X( 1, J ), 1 )
143 IF( XNORM.LE.ZERO ) THEN
144 RESID = ONE / EPS
145 ELSE
146 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
147 END IF
148 10 CONTINUE
149 *
150 RETURN
151 *
152 * End of ZTPT02
153 *
154 END
2 $ WORK, RWORK, 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 DIAG, TRANS, UPLO
10 INTEGER LDB, LDX, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 AP( * ), B( LDB, * ), WORK( * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZTPT02 computes the residual for the computed solution to a
22 * triangular system of linear equations A*x = b, A**T *x = b, or
23 * A**H *x = b, when the triangular matrix A is stored in packed format.
24 * Here A**T denotes the transpose of A, A**H denotes the conjugate
25 * transpose of A, and x and b are N by NRHS matrices. The test ratio
26 * is the maximum over the number of right hand sides of
27 * the maximum over the number of right hand sides of
28 * norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
29 * where op(A) denotes A, A**T, or A**H, and EPS is the machine epsilon.
30 *
31 * Arguments
32 * =========
33 *
34 * UPLO (input) CHARACTER*1
35 * Specifies whether the matrix A is upper or lower triangular.
36 * = 'U': Upper triangular
37 * = 'L': Lower triangular
38 *
39 * TRANS (input) CHARACTER*1
40 * Specifies the operation applied to A.
41 * = 'N': A *x = b (No transpose)
42 * = 'T': A**T *x = b (Transpose)
43 * = 'C': A**H *x = b (Conjugate transpose)
44 *
45 * DIAG (input) CHARACTER*1
46 * Specifies whether or not the matrix A is unit triangular.
47 * = 'N': Non-unit triangular
48 * = 'U': Unit triangular
49 *
50 * N (input) INTEGER
51 * The order of the matrix A. N >= 0.
52 *
53 * NRHS (input) INTEGER
54 * The number of right hand sides, i.e., the number of columns
55 * of the matrices X and B. NRHS >= 0.
56 *
57 * AP (input) COMPLEX*16 array, dimension (N*(N+1)/2)
58 * The upper or lower triangular matrix A, packed columnwise in
59 * a linear array. The j-th column of A is stored in the array
60 * AP as follows:
61 * if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
62 * if UPLO = 'L',
63 * AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
64 *
65 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
66 * The computed solution vectors for the system of linear
67 * equations.
68 *
69 * LDX (input) INTEGER
70 * The leading dimension of the array X. LDX >= max(1,N).
71 *
72 * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
73 * The right hand side vectors for the system of linear
74 * equations.
75 *
76 * LDB (input) INTEGER
77 * The leading dimension of the array B. LDB >= max(1,N).
78 *
79 * WORK (workspace) COMPLEX*16 array, dimension (N)
80 *
81 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
82 *
83 * RESID (output) DOUBLE PRECISION
84 * The maximum over the number of right hand sides of
85 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
86 *
87 * =====================================================================
88 *
89 * .. Parameters ..
90 DOUBLE PRECISION ZERO, ONE
91 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
92 * ..
93 * .. Local Scalars ..
94 INTEGER J
95 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
96 * ..
97 * .. External Functions ..
98 LOGICAL LSAME
99 DOUBLE PRECISION DLAMCH, DZASUM, ZLANTP
100 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANTP
101 * ..
102 * .. External Subroutines ..
103 EXTERNAL ZAXPY, ZCOPY, ZTPMV
104 * ..
105 * .. Intrinsic Functions ..
106 INTRINSIC DCMPLX, MAX
107 * ..
108 * .. Executable Statements ..
109 *
110 * Quick exit if N = 0 or NRHS = 0
111 *
112 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
113 RESID = ZERO
114 RETURN
115 END IF
116 *
117 * Compute the 1-norm of A or A**H.
118 *
119 IF( LSAME( TRANS, 'N' ) ) THEN
120 ANORM = ZLANTP( '1', UPLO, DIAG, N, AP, RWORK )
121 ELSE
122 ANORM = ZLANTP( 'I', UPLO, DIAG, N, AP, RWORK )
123 END IF
124 *
125 * Exit with RESID = 1/EPS if ANORM = 0.
126 *
127 EPS = DLAMCH( 'Epsilon' )
128 IF( ANORM.LE.ZERO ) THEN
129 RESID = ONE / EPS
130 RETURN
131 END IF
132 *
133 * Compute the maximum over the number of right hand sides of
134 * norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
135 *
136 RESID = ZERO
137 DO 10 J = 1, NRHS
138 CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
139 CALL ZTPMV( UPLO, TRANS, DIAG, N, AP, WORK, 1 )
140 CALL ZAXPY( N, DCMPLX( -ONE ), B( 1, J ), 1, WORK, 1 )
141 BNORM = DZASUM( N, WORK, 1 )
142 XNORM = DZASUM( N, X( 1, J ), 1 )
143 IF( XNORM.LE.ZERO ) THEN
144 RESID = ONE / EPS
145 ELSE
146 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
147 END IF
148 10 CONTINUE
149 *
150 RETURN
151 *
152 * End of ZTPT02
153 *
154 END