1 SUBROUTINE ZGET08( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB,
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * June 2010
7 *
8 * .. Scalar Arguments ..
9 CHARACTER TRANS
10 INTEGER LDA, LDB, LDX, M, 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 * ZGET08 computes the residual for a solution of a system of linear
22 * equations A*x = b or A'*x = b:
23 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANS (input) CHARACTER*1
30 * Specifies the form of the system of equations:
31 * = 'N': A *x = b
32 * = 'T': A^T*x = b, where A^T is the transpose of A
33 * = 'C': A^H*x = b, where A^H is the conjugate transpose of A
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of 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 (LDA,N)
46 * The original M x N matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
52 * The computed solution vectors for the system of linear
53 * equations.
54 *
55 * LDX (input) INTEGER
56 * The leading dimension of the array X. If TRANS = 'N',
57 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
58 *
59 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
60 * On entry, the right hand side vectors for the system of
61 * linear equations.
62 * On exit, B is overwritten with the difference B - A*X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. IF TRANS = 'N',
66 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * The maximum over the number of right hand sides of
72 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
79 COMPLEX*16 CONE
80 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
81 * ..
82 * .. Local Scalars ..
83 INTEGER J, N1, N2
84 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
85 COMPLEX*16 ZDUM
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 INTEGER IZAMAX
90 DOUBLE PRECISION DLAMCH, ZLANGE
91 EXTERNAL LSAME, IZAMAX, DLAMCH, ZLANGE
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZGEMM
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC ABS, DBLE, DIMAG, MAX
98 * ..
99 * .. Statement Functions ..
100 DOUBLE PRECISION CABS1
101 * ..
102 * .. Statement Function definitions ..
103 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
104 * ..
105 * .. Executable Statements ..
106 *
107 * Quick exit if M = 0 or N = 0 or NRHS = 0
108 *
109 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
110 RESID = ZERO
111 RETURN
112 END IF
113 *
114 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
115 N1 = N
116 N2 = M
117 ELSE
118 N1 = M
119 N2 = N
120 END IF
121 *
122 * Exit with RESID = 1/EPS if ANORM = 0.
123 *
124 EPS = DLAMCH( 'Epsilon' )
125 ANORM = ZLANGE( 'I', N1, N2, A, LDA, RWORK )
126 IF( ANORM.LE.ZERO ) THEN
127 RESID = ONE / EPS
128 RETURN
129 END IF
130 *
131 * Compute B - A*X (or B - A'*X ) and store in B.
132 *
133 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X,
134 $ LDX, CONE, B, LDB )
135 *
136 * Compute the maximum over the number of right hand sides of
137 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
138 *
139 RESID = ZERO
140 DO 10 J = 1, NRHS
141 BNORM = CABS1( B( IZAMAX( N1, B( 1, J ), 1 ), J ) )
142 XNORM = CABS1( X( IZAMAX( N2, X( 1, J ), 1 ), J ) )
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 ZGET02
153 *
154 END
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * June 2010
7 *
8 * .. Scalar Arguments ..
9 CHARACTER TRANS
10 INTEGER LDA, LDB, LDX, M, 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 * ZGET08 computes the residual for a solution of a system of linear
22 * equations A*x = b or A'*x = b:
23 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANS (input) CHARACTER*1
30 * Specifies the form of the system of equations:
31 * = 'N': A *x = b
32 * = 'T': A^T*x = b, where A^T is the transpose of A
33 * = 'C': A^H*x = b, where A^H is the conjugate transpose of A
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of 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 (LDA,N)
46 * The original M x N matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
52 * The computed solution vectors for the system of linear
53 * equations.
54 *
55 * LDX (input) INTEGER
56 * The leading dimension of the array X. If TRANS = 'N',
57 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
58 *
59 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
60 * On entry, the right hand side vectors for the system of
61 * linear equations.
62 * On exit, B is overwritten with the difference B - A*X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. IF TRANS = 'N',
66 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * The maximum over the number of right hand sides of
72 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
79 COMPLEX*16 CONE
80 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
81 * ..
82 * .. Local Scalars ..
83 INTEGER J, N1, N2
84 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
85 COMPLEX*16 ZDUM
86 * ..
87 * .. External Functions ..
88 LOGICAL LSAME
89 INTEGER IZAMAX
90 DOUBLE PRECISION DLAMCH, ZLANGE
91 EXTERNAL LSAME, IZAMAX, DLAMCH, ZLANGE
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZGEMM
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC ABS, DBLE, DIMAG, MAX
98 * ..
99 * .. Statement Functions ..
100 DOUBLE PRECISION CABS1
101 * ..
102 * .. Statement Function definitions ..
103 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
104 * ..
105 * .. Executable Statements ..
106 *
107 * Quick exit if M = 0 or N = 0 or NRHS = 0
108 *
109 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
110 RESID = ZERO
111 RETURN
112 END IF
113 *
114 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
115 N1 = N
116 N2 = M
117 ELSE
118 N1 = M
119 N2 = N
120 END IF
121 *
122 * Exit with RESID = 1/EPS if ANORM = 0.
123 *
124 EPS = DLAMCH( 'Epsilon' )
125 ANORM = ZLANGE( 'I', N1, N2, A, LDA, RWORK )
126 IF( ANORM.LE.ZERO ) THEN
127 RESID = ONE / EPS
128 RETURN
129 END IF
130 *
131 * Compute B - A*X (or B - A'*X ) and store in B.
132 *
133 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X,
134 $ LDX, CONE, B, LDB )
135 *
136 * Compute the maximum over the number of right hand sides of
137 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
138 *
139 RESID = ZERO
140 DO 10 J = 1, NRHS
141 BNORM = CABS1( B( IZAMAX( N1, B( 1, J ), 1 ), J ) )
142 XNORM = CABS1( X( IZAMAX( N2, X( 1, J ), 1 ), J ) )
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 ZGET02
153 *
154 END