1*> \brief \b DTRT02
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 DTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B,
12*                          LDB, WORK, RESID )
13*
14*       .. Scalar Arguments ..
15*       CHARACTER          DIAG, TRANS, UPLO
16*       INTEGER            LDA, LDB, LDX, N, NRHS
17*       DOUBLE PRECISION   RESID
18*       ..
19*       .. Array Arguments ..
20*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * ),
21*      $                   X( LDX, * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> DTRT02 computes the residual for the computed solution to a
31*> triangular system of linear equations  A*x = b  or  A'*x = b.
32*> Here A is a triangular matrix, A' is the transpose of A, and x and b
33*> are N by NRHS matrices.  The test ratio is the maximum over the
34*> number of right hand sides of
35*>    norm(b - op(A)*x) / ( norm(op(A)) * norm(x) * EPS ),
36*> where op(A) denotes A or A' and EPS is the machine epsilon.
37*> \endverbatim
38*
39*  Arguments:
40*  ==========
41*
42*> \param[in] UPLO
43*> \verbatim
44*>          UPLO is CHARACTER*1
45*>          Specifies whether the matrix A is upper or lower triangular.
46*>          = 'U':  Upper triangular
47*>          = 'L':  Lower triangular
48*> \endverbatim
49*>
50*> \param[in] TRANS
51*> \verbatim
52*>          TRANS is CHARACTER*1
53*>          Specifies the operation applied to A.
54*>          = 'N':  A *x = b  (No transpose)
55*>          = 'T':  A'*x = b  (Transpose)
56*>          = 'C':  A'*x = b  (Conjugate transpose = Transpose)
57*> \endverbatim
58*>
59*> \param[in] DIAG
60*> \verbatim
61*>          DIAG is CHARACTER*1
62*>          Specifies whether or not the matrix A is unit triangular.
63*>          = 'N':  Non-unit triangular
64*>          = 'U':  Unit triangular
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*>          N is INTEGER
70*>          The order of the matrix A.  N >= 0.
71*> \endverbatim
72*>
73*> \param[in] NRHS
74*> \verbatim
75*>          NRHS is INTEGER
76*>          The number of right hand sides, i.e., the number of columns
77*>          of the matrices X and B.  NRHS >= 0.
78*> \endverbatim
79*>
80*> \param[in] A
81*> \verbatim
82*>          A is DOUBLE PRECISION array, dimension (LDA,N)
83*>          The triangular matrix A.  If UPLO = 'U', the leading n by n
84*>          upper triangular part of the array A contains the upper
85*>          triangular matrix, and the strictly lower triangular part of
86*>          A is not referenced.  If UPLO = 'L', the leading n by n lower
87*>          triangular part of the array A contains the lower triangular
88*>          matrix, and the strictly upper triangular part of A is not
89*>          referenced.  If DIAG = 'U', the diagonal elements of A are
90*>          also not referenced and are assumed to be 1.
91*> \endverbatim
92*>
93*> \param[in] LDA
94*> \verbatim
95*>          LDA is INTEGER
96*>          The leading dimension of the array A.  LDA >= max(1,N).
97*> \endverbatim
98*>
99*> \param[in] X
100*> \verbatim
101*>          X is DOUBLE PRECISION array, dimension (LDX,NRHS)
102*>          The computed solution vectors for the system of linear
103*>          equations.
104*> \endverbatim
105*>
106*> \param[in] LDX
107*> \verbatim
108*>          LDX is INTEGER
109*>          The leading dimension of the array X.  LDX >= max(1,N).
110*> \endverbatim
111*>
112*> \param[in] B
113*> \verbatim
114*>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
115*>          The right hand side vectors for the system of linear
116*>          equations.
117*> \endverbatim
118*>
119*> \param[in] LDB
120*> \verbatim
121*>          LDB is INTEGER
122*>          The leading dimension of the array B.  LDB >= max(1,N).
123*> \endverbatim
124*>
125*> \param[out] WORK
126*> \verbatim
127*>          WORK is DOUBLE PRECISION array, dimension (N)
128*> \endverbatim
129*>
130*> \param[out] RESID
131*> \verbatim
132*>          RESID is DOUBLE PRECISION
133*>          The maximum over the number of right hand sides of
134*>          norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS ).
135*> \endverbatim
136*
137*  Authors:
138*  ========
139*
140*> \author Univ. of Tennessee
141*> \author Univ. of California Berkeley
142*> \author Univ. of Colorado Denver
143*> \author NAG Ltd.
144*
145*> \date November 2011
146*
147*> \ingroup double_lin
148*
149*  =====================================================================
150      SUBROUTINE DTRT02( UPLO, TRANS, DIAG, N, NRHS, A, LDA, X, LDX, B,
151     $                   LDB, WORK, RESID )
152*
153*  -- LAPACK test routine (version 3.4.0) --
154*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
155*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
156*     November 2011
157*
158*     .. Scalar Arguments ..
159      CHARACTER          DIAG, TRANS, UPLO
160      INTEGER            LDA, LDB, LDX, N, NRHS
161      DOUBLE PRECISION   RESID
162*     ..
163*     .. Array Arguments ..
164      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * ),
165     $                   X( LDX, * )
166*     ..
167*
168*  =====================================================================
169*
170*     .. Parameters ..
171      DOUBLE PRECISION   ZERO, ONE
172      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
173*     ..
174*     .. Local Scalars ..
175      INTEGER            J
176      DOUBLE PRECISION   ANORM, BNORM, EPS, XNORM
177*     ..
178*     .. External Functions ..
179      LOGICAL            LSAME
180      DOUBLE PRECISION   DASUM, DLAMCH, DLANTR
181      EXTERNAL           LSAME, DASUM, DLAMCH, DLANTR
182*     ..
183*     .. External Subroutines ..
184      EXTERNAL           DAXPY, DCOPY, DTRMV
185*     ..
186*     .. Intrinsic Functions ..
187      INTRINSIC          MAX
188*     ..
189*     .. Executable Statements ..
190*
191*     Quick exit if N = 0 or NRHS = 0
192*
193      IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
194         RESID = ZERO
195         RETURN
196      END IF
197*
198*     Compute the 1-norm of A or A'.
199*
200      IF( LSAME( TRANS, 'N' ) ) THEN
201         ANORM = DLANTR( '1', UPLO, DIAG, N, N, A, LDA, WORK )
202      ELSE
203         ANORM = DLANTR( 'I', UPLO, DIAG, N, N, A, LDA, WORK )
204      END IF
205*
206*     Exit with RESID = 1/EPS if ANORM = 0.
207*
208      EPS = DLAMCH( 'Epsilon' )
209      IF( ANORM.LE.ZERO ) THEN
210         RESID = ONE / EPS
211         RETURN
212      END IF
213*
214*     Compute the maximum over the number of right hand sides of
215*        norm(op(A)*x - b) / ( norm(op(A)) * norm(x) * EPS )
216*
217      RESID = ZERO
218      DO 10 J = 1, NRHS
219         CALL DCOPY( N, X( 1, J ), 1, WORK, 1 )
220         CALL DTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
221         CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
222         BNORM = DASUM( N, WORK, 1 )
223         XNORM = DASUM( N, X( 1, J ), 1 )
224         IF( XNORM.LE.ZERO ) THEN
225            RESID = ONE / EPS
226         ELSE
227            RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
228         END IF
229   10 CONTINUE
230*
231      RETURN
232*
233*     End of DTRT02
234*
235      END
236