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