1*> \brief \b DGET02
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*  Definition:
9*  ===========
10*
11*       SUBROUTINE DGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB,
12*                          RWORK, RESID )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          TRANS
16*       INTEGER            LDA, LDB, LDX, M, N, NRHS
17*       DOUBLE PRECISION   RESID
18*       ..
19*       .. Array Arguments ..
20*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), RWORK( * ),
21*      $                   X( LDX, * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> DGET02 computes the residual for a solution of a system of linear
31*> equations op(A)*X = B:
32*>    RESID = norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ),
33*> where op(A) = A or A**T, depending on TRANS, and EPS is the
34*> machine epsilon.
35*> The norm used is the 1-norm.
36*> \endverbatim
37*
38*  Arguments:
39*  ==========
40*
41*> \param[in] TRANS
42*> \verbatim
43*>          TRANS is CHARACTER*1
44*>          Specifies the form of the system of equations:
45*>          = 'N':  A    * X = B  (No transpose)
46*>          = 'T':  A**T * X = B  (Transpose)
47*>          = 'C':  A**H * X = B  (Conjugate transpose = Transpose)
48*> \endverbatim
49*>
50*> \param[in] M
51*> \verbatim
52*>          M is INTEGER
53*>          The number of rows of the matrix A.  M >= 0.
54*> \endverbatim
55*>
56*> \param[in] N
57*> \verbatim
58*>          N is INTEGER
59*>          The number of columns of the matrix A.  N >= 0.
60*> \endverbatim
61*>
62*> \param[in] NRHS
63*> \verbatim
64*>          NRHS is INTEGER
65*>          The number of columns of B, the matrix of right hand sides.
66*>          NRHS >= 0.
67*> \endverbatim
68*>
69*> \param[in] A
70*> \verbatim
71*>          A is DOUBLE PRECISION array, dimension (LDA,N)
72*>          The original M x N matrix A.
73*> \endverbatim
74*>
75*> \param[in] LDA
76*> \verbatim
77*>          LDA is INTEGER
78*>          The leading dimension of the array A.  LDA >= max(1,M).
79*> \endverbatim
80*>
81*> \param[in] X
82*> \verbatim
83*>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
84*>          The computed solution vectors for the system of linear
85*>          equations.
86*> \endverbatim
87*>
88*> \param[in] LDX
89*> \verbatim
90*>          LDX is INTEGER
91*>          The leading dimension of the array X.  If TRANS = 'N',
92*>          LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
93*> \endverbatim
94*>
95*> \param[in,out] B
96*> \verbatim
97*>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
98*>          On entry, the right hand side vectors for the system of
99*>          linear equations.
100*>          On exit, B is overwritten with the difference B - op(A)*X.
101*> \endverbatim
102*>
103*> \param[in] LDB
104*> \verbatim
105*>          LDB is INTEGER
106*>          The leading dimension of the array B.  IF TRANS = 'N',
107*>          LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
108*> \endverbatim
109*>
110*> \param[out] RWORK
111*> \verbatim
112*>          RWORK is DOUBLE PRECISION array, dimension (M)
113*> \endverbatim
114*>
115*> \param[out] RESID
116*> \verbatim
117*>          RESID is DOUBLE PRECISION
118*>          The maximum over the number of right hand sides of
119*>          norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ).
120*> \endverbatim
121*
122*  Authors:
123*  ========
124*
125*> \author Univ. of Tennessee
126*> \author Univ. of California Berkeley
127*> \author Univ. of Colorado Denver
128*> \author NAG Ltd.
129*
130*> \ingroup double_lin
131*
132*  =====================================================================
133      SUBROUTINE DGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB,
134     $                   RWORK, RESID )
135*
136*  -- LAPACK test routine --
137*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
138*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
139*
140*     .. Scalar Arguments ..
141      CHARACTER          TRANS
142      INTEGER            LDA, LDB, LDX, M, N, NRHS
143      DOUBLE PRECISION   RESID
144*     ..
145*     .. Array Arguments ..
146      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), RWORK( * ),
147     $                   X( LDX, * )
148*     ..
149*
150*  =====================================================================
151*
152*     .. Parameters ..
153      DOUBLE PRECISION   ZERO, ONE
154      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
155*     ..
156*     .. Local Scalars ..
157      INTEGER            J, N1, N2
158      DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM
159*     ..
160*     .. External Functions ..
161      LOGICAL            LSAME
162      DOUBLE PRECISION   DASUM, DLAMCH, DLANGE
163      EXTERNAL           LSAME, DASUM, DLAMCH, DLANGE
164*     ..
165*     .. External Subroutines ..
166      EXTERNAL           DGEMM
167*     ..
168*     .. Intrinsic Functions ..
169      INTRINSIC          MAX
170*     ..
171*     .. Executable Statements ..
172*
173*     Quick exit if M = 0 or N = 0 or NRHS = 0
174*
175      IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
176         RESID = ZERO
177         RETURN
178      END IF
179*
180      IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
181         N1 = N
182         N2 = M
183      ELSE
184         N1 = M
185         N2 = N
186      END IF
187*
188*     Exit with RESID = 1/EPS if ANORM = 0.
189*
190      EPS = DLAMCH( 'Epsilon' )
191      IF( LSAME( TRANS, 'N' ) ) THEN
192         ANORM = DLANGE( '1', M, N, A, LDA, RWORK )
193      ELSE
194         ANORM = DLANGE( 'I', M, N, A, LDA, RWORK )
195      END IF
196      IF( ANORM.LE.ZERO ) THEN
197         RESID = ONE / EPS
198         RETURN
199      END IF
200*
201*     Compute B - op(A)*X and store in B.
202*
203      CALL DGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X,
204     $            LDX, ONE, B, LDB )
205*
206*     Compute the maximum over the number of right hand sides of
207*        norm(B - op(A)*X) / ( norm(op(A)) * norm(X) * EPS ) .
208*
209      RESID = ZERO
210      DO 10 J = 1, NRHS
211         BNORM = DASUM( N1, B( 1, J ), 1 )
212         XNORM = DASUM( N2, X( 1, J ), 1 )
213         IF( XNORM.LE.ZERO ) THEN
214            RESID = ONE / EPS
215         ELSE
216            RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
217         END IF
218   10 CONTINUE
219*
220      RETURN
221*
222*     End of DGET02
223*
224      END
225