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