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