1 SUBROUTINE SQRT16( 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 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 REAL A( LDA, * ), B( LDB, * ), RWORK( * ),
15 $ X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * SQRT16 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'*x = b, where A' is the transpose of A
33 * = 'C': A'*x = b, where A' is the 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) REAL 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) REAL 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) REAL 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) REAL array, dimension (M)
69 *
70 * RESID (output) REAL
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 REAL ZERO, ONE
78 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
79 * ..
80 * .. Local Scalars ..
81 INTEGER J, N1, N2
82 REAL ANORM, BNORM, EPS, XNORM
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 REAL SASUM, SLAMCH, SLANGE
87 EXTERNAL LSAME, SASUM, SLAMCH, SLANGE
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL SGEMM
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC MAX
94 * ..
95 * .. Executable Statements ..
96 *
97 * Quick exit if M = 0 or N = 0 or NRHS = 0
98 *
99 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
100 RESID = ZERO
101 RETURN
102 END IF
103 *
104 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
105 ANORM = SLANGE( 'I', M, N, A, LDA, RWORK )
106 N1 = N
107 N2 = M
108 ELSE
109 ANORM = SLANGE( '1', M, N, A, LDA, RWORK )
110 N1 = M
111 N2 = N
112 END IF
113 *
114 EPS = SLAMCH( 'Epsilon' )
115 *
116 * Compute B - A*X (or B - A'*X ) and store in B.
117 *
118 CALL SGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X,
119 $ LDX, ONE, B, LDB )
120 *
121 * Compute the maximum over the number of right hand sides of
122 * norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ) .
123 *
124 RESID = ZERO
125 DO 10 J = 1, NRHS
126 BNORM = SASUM( N1, B( 1, J ), 1 )
127 XNORM = SASUM( N2, X( 1, J ), 1 )
128 IF( ANORM.EQ.ZERO .AND. BNORM.EQ.ZERO ) THEN
129 RESID = ZERO
130 ELSE IF( ANORM.LE.ZERO .OR. XNORM.LE.ZERO ) THEN
131 RESID = ONE / EPS
132 ELSE
133 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) /
134 $ ( MAX( M, N )*EPS ) )
135 END IF
136 10 CONTINUE
137 *
138 RETURN
139 *
140 * End of SQRT16
141 *
142 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 REAL RESID
12 * ..
13 * .. Array Arguments ..
14 REAL A( LDA, * ), B( LDB, * ), RWORK( * ),
15 $ X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * SQRT16 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'*x = b, where A' is the transpose of A
33 * = 'C': A'*x = b, where A' is the 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) REAL 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) REAL 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) REAL 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) REAL array, dimension (M)
69 *
70 * RESID (output) REAL
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 REAL ZERO, ONE
78 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
79 * ..
80 * .. Local Scalars ..
81 INTEGER J, N1, N2
82 REAL ANORM, BNORM, EPS, XNORM
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 REAL SASUM, SLAMCH, SLANGE
87 EXTERNAL LSAME, SASUM, SLAMCH, SLANGE
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL SGEMM
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC MAX
94 * ..
95 * .. Executable Statements ..
96 *
97 * Quick exit if M = 0 or N = 0 or NRHS = 0
98 *
99 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
100 RESID = ZERO
101 RETURN
102 END IF
103 *
104 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
105 ANORM = SLANGE( 'I', M, N, A, LDA, RWORK )
106 N1 = N
107 N2 = M
108 ELSE
109 ANORM = SLANGE( '1', M, N, A, LDA, RWORK )
110 N1 = M
111 N2 = N
112 END IF
113 *
114 EPS = SLAMCH( 'Epsilon' )
115 *
116 * Compute B - A*X (or B - A'*X ) and store in B.
117 *
118 CALL SGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X,
119 $ LDX, ONE, B, LDB )
120 *
121 * Compute the maximum over the number of right hand sides of
122 * norm(B - A*X) / ( max(m,n) * norm(A) * norm(X) * EPS ) .
123 *
124 RESID = ZERO
125 DO 10 J = 1, NRHS
126 BNORM = SASUM( N1, B( 1, J ), 1 )
127 XNORM = SASUM( N2, X( 1, J ), 1 )
128 IF( ANORM.EQ.ZERO .AND. BNORM.EQ.ZERO ) THEN
129 RESID = ZERO
130 ELSE IF( ANORM.LE.ZERO .OR. XNORM.LE.ZERO ) THEN
131 RESID = ONE / EPS
132 ELSE
133 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) /
134 $ ( MAX( M, N )*EPS ) )
135 END IF
136 10 CONTINUE
137 *
138 RETURN
139 *
140 * End of SQRT16
141 *
142 END