3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download ZGERFS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgerfs.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgerfs.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgerfs.f">
21 * SUBROUTINE ZGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
22 * X, LDX, FERR, BERR, WORK, RWORK, INFO )
24 * .. Scalar Arguments ..
26 * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
28 * .. Array Arguments ..
30 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
31 * COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
32 * $ WORK( * ), X( LDX, * )
41 *> ZGERFS improves the computed solution to a system of linear
42 *> equations and provides error bounds and backward error estimates for
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)
61 *> The order of the matrix A. N >= 0.
67 *> The number of right hand sides, i.e., the number of columns
68 *> of the matrices B and X. NRHS >= 0.
73 *> A is COMPLEX*16 array, dimension (LDA,N)
74 *> The original N-by-N matrix A.
80 *> The leading dimension of the array A. LDA >= max(1,N).
85 *> AF is COMPLEX*16 array, dimension (LDAF,N)
86 *> The factors L and U from the factorization A = P*L*U
87 *> as computed by ZGETRF.
93 *> The leading dimension of the array AF. LDAF >= max(1,N).
98 *> IPIV is INTEGER array, dimension (N)
99 *> The pivot indices from ZGETRF; for 1<=i<=N, row i of the
100 *> matrix was interchanged with row IPIV(i).
105 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
106 *> The right hand side matrix B.
112 *> The leading dimension of the array B. LDB >= max(1,N).
117 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
118 *> On entry, the solution matrix X, as computed by ZGETRS.
119 *> On exit, the improved solution matrix X.
125 *> The leading dimension of the array X. LDX >= max(1,N).
130 *> FERR is DOUBLE PRECISION 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.
143 *> BERR is DOUBLE PRECISION 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).
151 *> WORK is COMPLEX*16 array, dimension (2*N)
156 *> RWORK is DOUBLE PRECISION array, dimension (N)
162 *> = 0: successful exit
163 *> < 0: if INFO = -i, the i-th argument had an illegal value
166 *> \par Internal Parameters:
167 * =========================
170 *> ITMAX is the maximum number of steps of iterative refinement.
176 *> \author Univ. of Tennessee
177 *> \author Univ. of California Berkeley
178 *> \author Univ. of Colorado Denver
181 *> \date November 2011
183 *> \ingroup complex16GEcomputational
185 * =====================================================================
186 SUBROUTINE ZGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
187 $ X, LDX, FERR, BERR, WORK, RWORK, INFO )
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..--
194 * .. Scalar Arguments ..
196 INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
198 * .. Array Arguments ..
200 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
201 COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
202 $ WORK( * ), X( LDX, * )
205 * =====================================================================
209 PARAMETER ( ITMAX = 5 )
210 DOUBLE PRECISION ZERO
211 PARAMETER ( ZERO = 0.0D+0 )
213 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
215 PARAMETER ( TWO = 2.0D+0 )
216 DOUBLE PRECISION THREE
217 PARAMETER ( THREE = 3.0D+0 )
219 * .. Local Scalars ..
221 CHARACTER TRANSN, TRANST
222 INTEGER COUNT, I, J, K, KASE, NZ
223 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
229 * .. External Functions ..
231 DOUBLE PRECISION DLAMCH
232 EXTERNAL LSAME, DLAMCH
234 * .. External Subroutines ..
235 EXTERNAL XERBLA, ZAXPY, ZCOPY, ZGEMV, ZGETRS, ZLACN2
237 * .. Intrinsic Functions ..
238 INTRINSIC ABS, DBLE, DIMAG, MAX
240 * .. Statement Functions ..
241 DOUBLE PRECISION CABS1
243 * .. Statement Function definitions ..
244 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
246 * .. Executable Statements ..
248 * Test the input parameters.
251 NOTRAN = LSAME( TRANS, 'N' )
252 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
253 $ LSAME( TRANS, 'C' ) ) THEN
255 ELSE IF( N.LT.0 ) THEN
257 ELSE IF( NRHS.LT.0 ) THEN
259 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
261 ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
263 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
265 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
269 CALL XERBLA( 'ZGERFS', -INFO )
273 * Quick return if possible
275 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
291 * NZ = maximum number of nonzero elements in each row of A, plus 1
294 EPS = DLAMCH( 'Epsilon' )
295 SAFMIN = DLAMCH( 'Safe minimum' )
299 * Do for each right hand side
307 * Loop until stopping criterion is satisfied.
309 * Compute residual R = B - op(A) * X,
310 * where op(A) = A, A**T, or A**H, depending on TRANS.
312 CALL ZCOPY( N, B( 1, J ), 1, WORK, 1 )
313 CALL ZGEMV( TRANS, N, N, -ONE, A, LDA, X( 1, J ), 1, ONE, WORK,
316 * Compute componentwise relative backward error from formula
318 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
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.
326 RWORK( I ) = CABS1( B( I, J ) )
329 * Compute abs(op(A))*abs(X) + abs(B).
333 XK = CABS1( X( K, J ) )
335 RWORK( I ) = RWORK( I ) + CABS1( A( I, K ) )*XK
342 S = S + CABS1( A( I, K ) )*CABS1( X( I, J ) )
344 RWORK( K ) = RWORK( K ) + S
349 IF( RWORK( I ).GT.SAFE2 ) THEN
350 S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
352 S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
353 $ ( RWORK( I )+SAFE1 ) )
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.
364 IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
365 $ COUNT.LE.ITMAX ) THEN
367 * Update solution and try again.
369 CALL ZGETRS( TRANS, N, 1, AF, LDAF, IPIV, WORK, N, INFO )
370 CALL ZAXPY( N, ONE, WORK, 1, X( 1, J ), 1 )
376 * Bound error from formula
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)
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
387 * NZ is the maximum number of nonzeros in any row of A, plus 1
388 * EPS is machine epsilon
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.
394 * Use ZLACN2 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) )))
399 IF( RWORK( I ).GT.SAFE2 ) THEN
400 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
402 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
409 CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
413 * Multiply by diag(W)*inv(op(A)**H).
415 CALL ZGETRS( TRANST, N, 1, AF, LDAF, IPIV, WORK, N,
418 WORK( I ) = RWORK( I )*WORK( I )
422 * Multiply by inv(op(A))*diag(W).
425 WORK( I ) = RWORK( I )*WORK( I )
427 CALL ZGETRS( TRANSN, N, 1, AF, LDAF, IPIV, WORK, N,
437 LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
440 $ FERR( J ) = FERR( J ) / LSTRES