1 SUBROUTINE ZQRT16( 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 * November 2006
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 * ZQRT16 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) / ( max(m,n) * 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) / ( max(m,n) * 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 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
89 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANGE
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL ZGEMM
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 * Quick exit if M = 0 or N = 0 or NRHS = 0
100 *
101 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
102 RESID = ZERO
103 RETURN
104 END IF
105 *
106 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
107 ANORM = ZLANGE( 'I', M, N, A, LDA, RWORK )
108 N1 = N
109 N2 = M
110 ELSE
111 ANORM = ZLANGE( '1', M, N, A, LDA, RWORK )
112 N1 = M
113 N2 = N
114 END IF
115 *
116 EPS = DLAMCH( 'Epsilon' )
117 *
118 * Compute B - A*X (or B - A'*X ) and store in B.
119 *
120 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X,
121 $ LDX, CONE, B, LDB )
122 *
123 * Compute the maximum over the number of right hand sides of
124 * norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ) .
125 *
126 RESID = ZERO
127 DO 10 J = 1, NRHS
128 BNORM = DZASUM( N1, B( 1, J ), 1 )
129 XNORM = DZASUM( N2, X( 1, J ), 1 )
130 IF( ANORM.EQ.ZERO .AND. BNORM.EQ.ZERO ) THEN
131 RESID = ZERO
132 ELSE IF( ANORM.LE.ZERO .OR. XNORM.LE.ZERO ) THEN
133 RESID = ONE / EPS
134 ELSE
135 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) /
136 $ ( MAX( M, N )*EPS ) )
137 END IF
138 10 CONTINUE
139 *
140 RETURN
141 *
142 * End of ZQRT16
143 *
144 END
2 $ 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 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 * ZQRT16 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) / ( max(m,n) * 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) / ( max(m,n) * 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 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
89 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANGE
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL ZGEMM
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 * Quick exit if M = 0 or N = 0 or NRHS = 0
100 *
101 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
102 RESID = ZERO
103 RETURN
104 END IF
105 *
106 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
107 ANORM = ZLANGE( 'I', M, N, A, LDA, RWORK )
108 N1 = N
109 N2 = M
110 ELSE
111 ANORM = ZLANGE( '1', M, N, A, LDA, RWORK )
112 N1 = M
113 N2 = N
114 END IF
115 *
116 EPS = DLAMCH( 'Epsilon' )
117 *
118 * Compute B - A*X (or B - A'*X ) and store in B.
119 *
120 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X,
121 $ LDX, CONE, B, LDB )
122 *
123 * Compute the maximum over the number of right hand sides of
124 * norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ) .
125 *
126 RESID = ZERO
127 DO 10 J = 1, NRHS
128 BNORM = DZASUM( N1, B( 1, J ), 1 )
129 XNORM = DZASUM( N2, X( 1, J ), 1 )
130 IF( ANORM.EQ.ZERO .AND. BNORM.EQ.ZERO ) THEN
131 RESID = ZERO
132 ELSE IF( ANORM.LE.ZERO .OR. XNORM.LE.ZERO ) THEN
133 RESID = ONE / EPS
134 ELSE
135 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) /
136 $ ( MAX( M, N )*EPS ) )
137 END IF
138 10 CONTINUE
139 *
140 RETURN
141 *
142 * End of ZQRT16
143 *
144 END