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