LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cgerfs.f
Go to the documentation of this file.
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 *> \ingroup complexGEcomputational
182 *
183 * =====================================================================
184  SUBROUTINE cgerfs( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
185  $ X, LDX, FERR, BERR, WORK, RWORK, INFO )
186 *
187 * -- LAPACK computational routine --
188 * -- LAPACK is a software package provided by Univ. of Tennessee, --
189 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
190 *
191 * .. Scalar Arguments ..
192  CHARACTER TRANS
193  INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
194 * ..
195 * .. Array Arguments ..
196  INTEGER IPIV( * )
197  REAL BERR( * ), FERR( * ), RWORK( * )
198  COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
199  $ work( * ), x( ldx, * )
200 * ..
201 *
202 * =====================================================================
203 *
204 * .. Parameters ..
205  INTEGER ITMAX
206  parameter( itmax = 5 )
207  REAL ZERO
208  parameter( zero = 0.0e+0 )
209  COMPLEX ONE
210  parameter( one = ( 1.0e+0, 0.0e+0 ) )
211  REAL TWO
212  parameter( two = 2.0e+0 )
213  REAL THREE
214  parameter( three = 3.0e+0 )
215 * ..
216 * .. Local Scalars ..
217  LOGICAL NOTRAN
218  CHARACTER TRANSN, TRANST
219  INTEGER COUNT, I, J, K, KASE, NZ
220  REAL EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
221  COMPLEX ZDUM
222 * ..
223 * .. Local Arrays ..
224  INTEGER ISAVE( 3 )
225 * ..
226 * .. External Functions ..
227  LOGICAL LSAME
228  REAL SLAMCH
229  EXTERNAL lsame, slamch
230 * ..
231 * .. External Subroutines ..
232  EXTERNAL caxpy, ccopy, cgemv, cgetrs, clacn2, xerbla
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC abs, aimag, max, real
236 * ..
237 * .. Statement Functions ..
238  REAL CABS1
239 * ..
240 * .. Statement Function definitions ..
241  cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
242 * ..
243 * .. Executable Statements ..
244 *
245 * Test the input parameters.
246 *
247  info = 0
248  notran = lsame( trans, 'N' )
249  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
250  $ lsame( trans, 'C' ) ) THEN
251  info = -1
252  ELSE IF( n.LT.0 ) THEN
253  info = -2
254  ELSE IF( nrhs.LT.0 ) THEN
255  info = -3
256  ELSE IF( lda.LT.max( 1, n ) ) THEN
257  info = -5
258  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
259  info = -7
260  ELSE IF( ldb.LT.max( 1, n ) ) THEN
261  info = -10
262  ELSE IF( ldx.LT.max( 1, n ) ) THEN
263  info = -12
264  END IF
265  IF( info.NE.0 ) THEN
266  CALL xerbla( 'CGERFS', -info )
267  RETURN
268  END IF
269 *
270 * Quick return if possible
271 *
272  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
273  DO 10 j = 1, nrhs
274  ferr( j ) = zero
275  berr( j ) = zero
276  10 CONTINUE
277  RETURN
278  END IF
279 *
280  IF( notran ) THEN
281  transn = 'N'
282  transt = 'C'
283  ELSE
284  transn = 'C'
285  transt = 'N'
286  END IF
287 *
288 * NZ = maximum number of nonzero elements in each row of A, plus 1
289 *
290  nz = n + 1
291  eps = slamch( 'Epsilon' )
292  safmin = slamch( 'Safe minimum' )
293  safe1 = nz*safmin
294  safe2 = safe1 / eps
295 *
296 * Do for each right hand side
297 *
298  DO 140 j = 1, nrhs
299 *
300  count = 1
301  lstres = three
302  20 CONTINUE
303 *
304 * Loop until stopping criterion is satisfied.
305 *
306 * Compute residual R = B - op(A) * X,
307 * where op(A) = A, A**T, or A**H, depending on TRANS.
308 *
309  CALL ccopy( n, b( 1, j ), 1, work, 1 )
310  CALL cgemv( trans, n, n, -one, a, lda, x( 1, j ), 1, one, work,
311  $ 1 )
312 *
313 * Compute componentwise relative backward error from formula
314 *
315 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
316 *
317 * where abs(Z) is the componentwise absolute value of the matrix
318 * or vector Z. If the i-th component of the denominator is less
319 * than SAFE2, then SAFE1 is added to the i-th components of the
320 * numerator and denominator before dividing.
321 *
322  DO 30 i = 1, n
323  rwork( i ) = cabs1( b( i, j ) )
324  30 CONTINUE
325 *
326 * Compute abs(op(A))*abs(X) + abs(B).
327 *
328  IF( notran ) THEN
329  DO 50 k = 1, n
330  xk = cabs1( x( k, j ) )
331  DO 40 i = 1, n
332  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
333  40 CONTINUE
334  50 CONTINUE
335  ELSE
336  DO 70 k = 1, n
337  s = zero
338  DO 60 i = 1, n
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 cgetrs( trans, 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(op(A)))*
377 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
378 *
379 * where
380 * norm(Z) is the magnitude of the largest component of Z
381 * inv(op(A)) is the inverse of op(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(op(A))*abs(X)+abs(B))
388 * is incremented by SAFE1 if the i-th component of
389 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
390 *
391 * Use CLACN2 to estimate the infinity-norm of the matrix
392 * inv(op(A)) * diag(W),
393 * where W = abs(R) + NZ*EPS*( abs(op(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(op(A)**H).
411 *
412  CALL cgetrs( transt, n, 1, af, ldaf, ipiv, work, n,
413  $ info )
414  DO 110 i = 1, n
415  work( i ) = rwork( i )*work( i )
416  110 CONTINUE
417  ELSE
418 *
419 * Multiply by inv(op(A))*diag(W).
420 *
421  DO 120 i = 1, n
422  work( i ) = rwork( i )*work( i )
423  120 CONTINUE
424  CALL cgetrs( transn, n, 1, af, ldaf, ipiv, work, n,
425  $ info )
426  END IF
427  GO TO 100
428  END IF
429 *
430 * Normalize error.
431 *
432  lstres = zero
433  DO 130 i = 1, n
434  lstres = max( lstres, cabs1( x( i, j ) ) )
435  130 CONTINUE
436  IF( lstres.NE.zero )
437  $ ferr( j ) = ferr( j ) / lstres
438 *
439  140 CONTINUE
440 *
441  RETURN
442 *
443 * End of CGERFS
444 *
445  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine cgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
CGETRS
Definition: cgetrs.f:121
subroutine cgerfs(TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO)
CGERFS
Definition: cgerfs.f:186
subroutine clacn2(N, V, X, EST, KASE, ISAVE)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition: clacn2.f:133