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