1*> \brief \b CTRRFS
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTRRFS + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrrfs.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrrfs.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrrfs.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
22*                          LDX, FERR, BERR, WORK, RWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          DIAG, TRANS, UPLO
26*       INTEGER            INFO, LDA, LDB, LDX, N, NRHS
27*       ..
28*       .. Array Arguments ..
29*       REAL               BERR( * ), FERR( * ), RWORK( * )
30*       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * ),
31*      $                   X( LDX, * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> CTRRFS provides error bounds and backward error estimates for the
41*> solution to a system of linear equations with a triangular
42*> coefficient matrix.
43*>
44*> The solution matrix X must be computed by CTRTRS or some other
45*> means before entering this routine.  CTRRFS does not do iterative
46*> refinement because doing so cannot improve the backward error.
47*> \endverbatim
48*
49*  Arguments:
50*  ==========
51*
52*> \param[in] UPLO
53*> \verbatim
54*>          UPLO is CHARACTER*1
55*>          = 'U':  A is upper triangular;
56*>          = 'L':  A is lower triangular.
57*> \endverbatim
58*>
59*> \param[in] TRANS
60*> \verbatim
61*>          TRANS is CHARACTER*1
62*>          Specifies the form of the system of equations:
63*>          = 'N':  A * X = B     (No transpose)
64*>          = 'T':  A**T * X = B  (Transpose)
65*>          = 'C':  A**H * X = B  (Conjugate transpose)
66*> \endverbatim
67*>
68*> \param[in] DIAG
69*> \verbatim
70*>          DIAG is CHARACTER*1
71*>          = 'N':  A is non-unit triangular;
72*>          = 'U':  A is unit triangular.
73*> \endverbatim
74*>
75*> \param[in] N
76*> \verbatim
77*>          N is INTEGER
78*>          The order of the matrix A.  N >= 0.
79*> \endverbatim
80*>
81*> \param[in] NRHS
82*> \verbatim
83*>          NRHS is INTEGER
84*>          The number of right hand sides, i.e., the number of columns
85*>          of the matrices B and X.  NRHS >= 0.
86*> \endverbatim
87*>
88*> \param[in] A
89*> \verbatim
90*>          A is COMPLEX array, dimension (LDA,N)
91*>          The triangular matrix A.  If UPLO = 'U', the leading N-by-N
92*>          upper triangular part of the array A contains the upper
93*>          triangular matrix, and the strictly lower triangular part of
94*>          A is not referenced.  If UPLO = 'L', the leading N-by-N lower
95*>          triangular part of the array A contains the lower triangular
96*>          matrix, and the strictly upper triangular part of A is not
97*>          referenced.  If DIAG = 'U', the diagonal elements of A are
98*>          also not referenced and are assumed to be 1.
99*> \endverbatim
100*>
101*> \param[in] LDA
102*> \verbatim
103*>          LDA is INTEGER
104*>          The leading dimension of the array A.  LDA >= max(1,N).
105*> \endverbatim
106*>
107*> \param[in] B
108*> \verbatim
109*>          B is COMPLEX array, dimension (LDB,NRHS)
110*>          The right hand side matrix B.
111*> \endverbatim
112*>
113*> \param[in] LDB
114*> \verbatim
115*>          LDB is INTEGER
116*>          The leading dimension of the array B.  LDB >= max(1,N).
117*> \endverbatim
118*>
119*> \param[in] X
120*> \verbatim
121*>          X is COMPLEX array, dimension (LDX,NRHS)
122*>          The solution matrix X.
123*> \endverbatim
124*>
125*> \param[in] LDX
126*> \verbatim
127*>          LDX is INTEGER
128*>          The leading dimension of the array X.  LDX >= max(1,N).
129*> \endverbatim
130*>
131*> \param[out] FERR
132*> \verbatim
133*>          FERR is REAL array, dimension (NRHS)
134*>          The estimated forward error bound for each solution vector
135*>          X(j) (the j-th column of the solution matrix X).
136*>          If XTRUE is the true solution corresponding to X(j), FERR(j)
137*>          is an estimated upper bound for the magnitude of the largest
138*>          element in (X(j) - XTRUE) divided by the magnitude of the
139*>          largest element in X(j).  The estimate is as reliable as
140*>          the estimate for RCOND, and is almost always a slight
141*>          overestimate of the true error.
142*> \endverbatim
143*>
144*> \param[out] BERR
145*> \verbatim
146*>          BERR is REAL array, dimension (NRHS)
147*>          The componentwise relative backward error of each solution
148*>          vector X(j) (i.e., the smallest relative change in
149*>          any element of A or B that makes X(j) an exact solution).
150*> \endverbatim
151*>
152*> \param[out] WORK
153*> \verbatim
154*>          WORK is COMPLEX array, dimension (2*N)
155*> \endverbatim
156*>
157*> \param[out] RWORK
158*> \verbatim
159*>          RWORK is REAL array, dimension (N)
160*> \endverbatim
161*>
162*> \param[out] INFO
163*> \verbatim
164*>          INFO is INTEGER
165*>          = 0:  successful exit
166*>          < 0:  if INFO = -i, the i-th argument had an illegal value
167*> \endverbatim
168*
169*  Authors:
170*  ========
171*
172*> \author Univ. of Tennessee
173*> \author Univ. of California Berkeley
174*> \author Univ. of Colorado Denver
175*> \author NAG Ltd.
176*
177*> \ingroup complexOTHERcomputational
178*
179*  =====================================================================
180      SUBROUTINE CTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
181     $                   LDX, FERR, BERR, WORK, RWORK, INFO )
182*
183*  -- LAPACK computational routine --
184*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
185*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187*     .. Scalar Arguments ..
188      CHARACTER          DIAG, TRANS, UPLO
189      INTEGER            INFO, LDA, LDB, LDX, N, NRHS
190*     ..
191*     .. Array Arguments ..
192      REAL               BERR( * ), FERR( * ), RWORK( * )
193      COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * ),
194     $                   X( LDX, * )
195*     ..
196*
197*  =====================================================================
198*
199*     .. Parameters ..
200      REAL               ZERO
201      PARAMETER          ( ZERO = 0.0E+0 )
202      COMPLEX            ONE
203      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
204*     ..
205*     .. Local Scalars ..
206      LOGICAL            NOTRAN, NOUNIT, UPPER
207      CHARACTER          TRANSN, TRANST
208      INTEGER            I, J, K, KASE, NZ
209      REAL               EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
210      COMPLEX            ZDUM
211*     ..
212*     .. Local Arrays ..
213      INTEGER            ISAVE( 3 )
214*     ..
215*     .. External Subroutines ..
216      EXTERNAL           CAXPY, CCOPY, CLACN2, CTRMV, CTRSV, XERBLA
217*     ..
218*     .. Intrinsic Functions ..
219      INTRINSIC          ABS, AIMAG, MAX, REAL
220*     ..
221*     .. External Functions ..
222      LOGICAL            LSAME
223      REAL               SLAMCH
224      EXTERNAL           LSAME, SLAMCH
225*     ..
226*     .. Statement Functions ..
227      REAL               CABS1
228*     ..
229*     .. Statement Function definitions ..
230      CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
231*     ..
232*     .. Executable Statements ..
233*
234*     Test the input parameters.
235*
236      INFO = 0
237      UPPER = LSAME( UPLO, 'U' )
238      NOTRAN = LSAME( TRANS, 'N' )
239      NOUNIT = LSAME( DIAG, 'N' )
240*
241      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
242         INFO = -1
243      ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
244     $         LSAME( TRANS, 'C' ) ) THEN
245         INFO = -2
246      ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
247         INFO = -3
248      ELSE IF( N.LT.0 ) THEN
249         INFO = -4
250      ELSE IF( NRHS.LT.0 ) THEN
251         INFO = -5
252      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
253         INFO = -7
254      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
255         INFO = -9
256      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
257         INFO = -11
258      END IF
259      IF( INFO.NE.0 ) THEN
260         CALL XERBLA( 'CTRRFS', -INFO )
261         RETURN
262      END IF
263*
264*     Quick return if possible
265*
266      IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
267         DO 10 J = 1, NRHS
268            FERR( J ) = ZERO
269            BERR( J ) = ZERO
270   10    CONTINUE
271         RETURN
272      END IF
273*
274      IF( NOTRAN ) THEN
275         TRANSN = 'N'
276         TRANST = 'C'
277      ELSE
278         TRANSN = 'C'
279         TRANST = 'N'
280      END IF
281*
282*     NZ = maximum number of nonzero elements in each row of A, plus 1
283*
284      NZ = N + 1
285      EPS = SLAMCH( 'Epsilon' )
286      SAFMIN = SLAMCH( 'Safe minimum' )
287      SAFE1 = NZ*SAFMIN
288      SAFE2 = SAFE1 / EPS
289*
290*     Do for each right hand side
291*
292      DO 250 J = 1, NRHS
293*
294*        Compute residual R = B - op(A) * X,
295*        where op(A) = A, A**T, or A**H, depending on TRANS.
296*
297         CALL CCOPY( N, X( 1, J ), 1, WORK, 1 )
298         CALL CTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK, 1 )
299         CALL CAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
300*
301*        Compute componentwise relative backward error from formula
302*
303*        max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
304*
305*        where abs(Z) is the componentwise absolute value of the matrix
306*        or vector Z.  If the i-th component of the denominator is less
307*        than SAFE2, then SAFE1 is added to the i-th components of the
308*        numerator and denominator before dividing.
309*
310         DO 20 I = 1, N
311            RWORK( I ) = CABS1( B( I, J ) )
312   20    CONTINUE
313*
314         IF( NOTRAN ) THEN
315*
316*           Compute abs(A)*abs(X) + abs(B).
317*
318            IF( UPPER ) THEN
319               IF( NOUNIT ) THEN
320                  DO 40 K = 1, N
321                     XK = CABS1( X( K, J ) )
322                     DO 30 I = 1, K
323                        RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
324   30                CONTINUE
325   40             CONTINUE
326               ELSE
327                  DO 60 K = 1, N
328                     XK = CABS1( X( K, J ) )
329                     DO 50 I = 1, K - 1
330                        RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
331   50                CONTINUE
332                     RWORK( K ) = RWORK( K ) + XK
333   60             CONTINUE
334               END IF
335            ELSE
336               IF( NOUNIT ) THEN
337                  DO 80 K = 1, N
338                     XK = CABS1( X( K, J ) )
339                     DO 70 I = K, N
340                        RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
341   70                CONTINUE
342   80             CONTINUE
343               ELSE
344                  DO 100 K = 1, N
345                     XK = CABS1( X( K, J ) )
346                     DO 90 I = K + 1, N
347                        RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
348   90                CONTINUE
349                     RWORK( K ) = RWORK( K ) + XK
350  100             CONTINUE
351               END IF
352            END IF
353         ELSE
354*
355*           Compute abs(A**H)*abs(X) + abs(B).
356*
357            IF( UPPER ) THEN
358               IF( NOUNIT ) THEN
359                  DO 120 K = 1, N
360                     S = ZERO
361                     DO 110 I = 1, K
362                        S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
363  110                CONTINUE
364                     RWORK( K ) = RWORK( K ) + S
365  120             CONTINUE
366               ELSE
367                  DO 140 K = 1, N
368                     S = CABS1( X( K, J ) )
369                     DO 130 I = 1, K - 1
370                        S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
371  130                CONTINUE
372                     RWORK( K ) = RWORK( K ) + S
373  140             CONTINUE
374               END IF
375            ELSE
376               IF( NOUNIT ) THEN
377                  DO 160 K = 1, N
378                     S = ZERO
379                     DO 150 I = K, N
380                        S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
381  150                CONTINUE
382                     RWORK( K ) = RWORK( K ) + S
383  160             CONTINUE
384               ELSE
385                  DO 180 K = 1, N
386                     S = CABS1( X( K, J ) )
387                     DO 170 I = K + 1, N
388                        S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
389  170                CONTINUE
390                     RWORK( K ) = RWORK( K ) + S
391  180             CONTINUE
392               END IF
393            END IF
394         END IF
395         S = ZERO
396         DO 190 I = 1, N
397            IF( RWORK( I ).GT.SAFE2 ) THEN
398               S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
399            ELSE
400               S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
401     $             ( RWORK( I )+SAFE1 ) )
402            END IF
403  190    CONTINUE
404         BERR( J ) = S
405*
406*        Bound error from formula
407*
408*        norm(X - XTRUE) / norm(X) .le. FERR =
409*        norm( abs(inv(op(A)))*
410*           ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
411*
412*        where
413*          norm(Z) is the magnitude of the largest component of Z
414*          inv(op(A)) is the inverse of op(A)
415*          abs(Z) is the componentwise absolute value of the matrix or
416*             vector Z
417*          NZ is the maximum number of nonzeros in any row of A, plus 1
418*          EPS is machine epsilon
419*
420*        The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
421*        is incremented by SAFE1 if the i-th component of
422*        abs(op(A))*abs(X) + abs(B) is less than SAFE2.
423*
424*        Use CLACN2 to estimate the infinity-norm of the matrix
425*           inv(op(A)) * diag(W),
426*        where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
427*
428         DO 200 I = 1, N
429            IF( RWORK( I ).GT.SAFE2 ) THEN
430               RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
431            ELSE
432               RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
433     $                      SAFE1
434            END IF
435  200    CONTINUE
436*
437         KASE = 0
438  210    CONTINUE
439         CALL CLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
440         IF( KASE.NE.0 ) THEN
441            IF( KASE.EQ.1 ) THEN
442*
443*              Multiply by diag(W)*inv(op(A)**H).
444*
445               CALL CTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK, 1 )
446               DO 220 I = 1, N
447                  WORK( I ) = RWORK( I )*WORK( I )
448  220          CONTINUE
449            ELSE
450*
451*              Multiply by inv(op(A))*diag(W).
452*
453               DO 230 I = 1, N
454                  WORK( I ) = RWORK( I )*WORK( I )
455  230          CONTINUE
456               CALL CTRSV( UPLO, TRANSN, DIAG, N, A, LDA, WORK, 1 )
457            END IF
458            GO TO 210
459         END IF
460*
461*        Normalize error.
462*
463         LSTRES = ZERO
464         DO 240 I = 1, N
465            LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
466  240    CONTINUE
467         IF( LSTRES.NE.ZERO )
468     $      FERR( J ) = FERR( J ) / LSTRES
469*
470  250 CONTINUE
471*
472      RETURN
473*
474*     End of CTRRFS
475*
476      END
477