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