LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *> \date November 2011
182 *
183 *> \ingroup complexGEcomputational
184 *
185 * =====================================================================
186  SUBROUTINE cgerfs( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
187  $ x, ldx, ferr, berr, work, rwork, info )
188 *
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..--
192 * November 2011
193 *
194 * .. Scalar Arguments ..
195  CHARACTER trans
196  INTEGER info, lda, ldaf, ldb, ldx, n, nrhs
197 * ..
198 * .. Array Arguments ..
199  INTEGER ipiv( * )
200  REAL berr( * ), ferr( * ), rwork( * )
201  COMPLEX a( lda, * ), af( ldaf, * ), b( ldb, * ),
202  $ work( * ), x( ldx, * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  INTEGER itmax
209  parameter( itmax = 5 )
210  REAL zero
211  parameter( zero = 0.0e+0 )
212  COMPLEX one
213  parameter( one = ( 1.0e+0, 0.0e+0 ) )
214  REAL two
215  parameter( two = 2.0e+0 )
216  REAL three
217  parameter( three = 3.0e+0 )
218 * ..
219 * .. Local Scalars ..
220  LOGICAL notran
221  CHARACTER transn, transt
222  INTEGER count, i, j, k, kase, nz
223  REAL eps, lstres, s, safe1, safe2, safmin, xk
224  COMPLEX zdum
225 * ..
226 * .. Local Arrays ..
227  INTEGER isave( 3 )
228 * ..
229 * .. External Functions ..
230  LOGICAL lsame
231  REAL slamch
232  EXTERNAL lsame, slamch
233 * ..
234 * .. External Subroutines ..
235  EXTERNAL caxpy, ccopy, cgemv, cgetrs, clacn2, xerbla
236 * ..
237 * .. Intrinsic Functions ..
238  INTRINSIC abs, aimag, max, real
239 * ..
240 * .. Statement Functions ..
241  REAL cabs1
242 * ..
243 * .. Statement Function definitions ..
244  cabs1( zdum ) = abs( REAL( ZDUM ) ) + abs( aimag( zdum ) )
245 * ..
246 * .. Executable Statements ..
247 *
248 * Test the input parameters.
249 *
250  info = 0
251  notran = lsame( trans, 'N' )
252  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
253  $ lsame( trans, 'C' ) ) THEN
254  info = -1
255  ELSE IF( n.LT.0 ) THEN
256  info = -2
257  ELSE IF( nrhs.LT.0 ) THEN
258  info = -3
259  ELSE IF( lda.LT.max( 1, n ) ) THEN
260  info = -5
261  ELSE IF( ldaf.LT.max( 1, n ) ) THEN
262  info = -7
263  ELSE IF( ldb.LT.max( 1, n ) ) THEN
264  info = -10
265  ELSE IF( ldx.LT.max( 1, n ) ) THEN
266  info = -12
267  END IF
268  IF( info.NE.0 ) THEN
269  CALL xerbla( 'CGERFS', -info )
270  RETURN
271  END IF
272 *
273 * Quick return if possible
274 *
275  IF( n.EQ.0 .OR. nrhs.EQ.0 ) THEN
276  DO 10 j = 1, nrhs
277  ferr( j ) = zero
278  berr( j ) = zero
279  10 CONTINUE
280  RETURN
281  END IF
282 *
283  IF( notran ) THEN
284  transn = 'N'
285  transt = 'C'
286  ELSE
287  transn = 'C'
288  transt = 'N'
289  END IF
290 *
291 * NZ = maximum number of nonzero elements in each row of A, plus 1
292 *
293  nz = n + 1
294  eps = slamch( 'Epsilon' )
295  safmin = slamch( 'Safe minimum' )
296  safe1 = nz*safmin
297  safe2 = safe1 / eps
298 *
299 * Do for each right hand side
300 *
301  DO 140 j = 1, nrhs
302 *
303  count = 1
304  lstres = three
305  20 CONTINUE
306 *
307 * Loop until stopping criterion is satisfied.
308 *
309 * Compute residual R = B - op(A) * X,
310 * where op(A) = A, A**T, or A**H, depending on TRANS.
311 *
312  CALL ccopy( n, b( 1, j ), 1, work, 1 )
313  CALL cgemv( trans, n, n, -one, a, lda, x( 1, j ), 1, one, work,
314  $ 1 )
315 *
316 * Compute componentwise relative backward error from formula
317 *
318 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
319 *
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.
324 *
325  DO 30 i = 1, n
326  rwork( i ) = cabs1( b( i, j ) )
327  30 CONTINUE
328 *
329 * Compute abs(op(A))*abs(X) + abs(B).
330 *
331  IF( notran ) THEN
332  DO 50 k = 1, n
333  xk = cabs1( x( k, j ) )
334  DO 40 i = 1, n
335  rwork( i ) = rwork( i ) + cabs1( a( i, k ) )*xk
336  40 CONTINUE
337  50 CONTINUE
338  ELSE
339  DO 70 k = 1, n
340  s = zero
341  DO 60 i = 1, n
342  s = s + cabs1( a( i, k ) )*cabs1( x( i, j ) )
343  60 CONTINUE
344  rwork( k ) = rwork( k ) + s
345  70 CONTINUE
346  END IF
347  s = zero
348  DO 80 i = 1, n
349  IF( rwork( i ).GT.safe2 ) THEN
350  s = max( s, cabs1( work( i ) ) / rwork( i ) )
351  ELSE
352  s = max( s, ( cabs1( work( i ) )+safe1 ) /
353  $ ( rwork( i )+safe1 ) )
354  END IF
355  80 CONTINUE
356  berr( j ) = s
357 *
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.
363 *
364  IF( berr( j ).GT.eps .AND. two*berr( j ).LE.lstres .AND.
365  $ count.LE.itmax ) THEN
366 *
367 * Update solution and try again.
368 *
369  CALL cgetrs( trans, n, 1, af, ldaf, ipiv, work, n, info )
370  CALL caxpy( n, one, work, 1, x( 1, j ), 1 )
371  lstres = berr( j )
372  count = count + 1
373  go to 20
374  END IF
375 *
376 * Bound error from formula
377 *
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)
381 *
382 * where
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
386 * vector Z
387 * NZ is the maximum number of nonzeros in any row of A, plus 1
388 * EPS is machine epsilon
389 *
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.
393 *
394 * Use CLACN2 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) )))
397 *
398  DO 90 i = 1, n
399  IF( rwork( i ).GT.safe2 ) THEN
400  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i )
401  ELSE
402  rwork( i ) = cabs1( work( i ) ) + nz*eps*rwork( i ) +
403  $ safe1
404  END IF
405  90 CONTINUE
406 *
407  kase = 0
408  100 CONTINUE
409  CALL clacn2( n, work( n+1 ), work, ferr( j ), kase, isave )
410  IF( kase.NE.0 ) THEN
411  IF( kase.EQ.1 ) THEN
412 *
413 * Multiply by diag(W)*inv(op(A)**H).
414 *
415  CALL cgetrs( transt, n, 1, af, ldaf, ipiv, work, n,
416  $ info )
417  DO 110 i = 1, n
418  work( i ) = rwork( i )*work( i )
419  110 CONTINUE
420  ELSE
421 *
422 * Multiply by inv(op(A))*diag(W).
423 *
424  DO 120 i = 1, n
425  work( i ) = rwork( i )*work( i )
426  120 CONTINUE
427  CALL cgetrs( transn, n, 1, af, ldaf, ipiv, work, n,
428  $ info )
429  END IF
430  go to 100
431  END IF
432 *
433 * Normalize error.
434 *
435  lstres = zero
436  DO 130 i = 1, n
437  lstres = max( lstres, cabs1( x( i, j ) ) )
438  130 CONTINUE
439  IF( lstres.NE.zero )
440  $ ferr( j ) = ferr( j ) / lstres
441 *
442  140 CONTINUE
443 *
444  RETURN
445 *
446 * End of CGERFS
447 *
448  END