1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
SUBROUTINE ZLA_LIN_BERR ( N, NZ, NRHS, RES, AYB, BERR )
*
* -- LAPACK routine (version 3.2.2) --
* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
* -- Jason Riedy of Univ. of California Berkeley. --
* -- June 2010 --
*
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley and NAG Ltd. --
*
IMPLICIT NONE
* ..
* .. Scalar Arguments ..
INTEGER N, NZ, NRHS
* ..
* .. Array Arguments ..
DOUBLE PRECISION AYB( N, NRHS ), BERR( NRHS )
COMPLEX*16 RES( N, NRHS )
* ..
*
* Purpose
* =======
*
* ZLA_LIN_BERR computes componentwise relative backward error from
* the formula
* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
* where abs(Z) is the componentwise absolute value of the matrix
* or vector Z.
*
* Arguments
* =========
*
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
*
* NZ (input) INTEGER
* We add (NZ+1)*SLAMCH( 'Safe minimum' ) to R(i) in the numerator to
* guard against spuriously zero residuals. Default value is N.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices AYB, RES, and BERR. NRHS >= 0.
*
* RES (input) DOUBLE PRECISION array, dimension (N,NRHS)
* The residual matrix, i.e., the matrix R in the relative backward
* error formula above.
*
* AYB (input) DOUBLE PRECISION array, dimension (N, NRHS)
* The denominator in the relative backward error formula above, i.e.,
* the matrix abs(op(A_s))*abs(Y) + abs(B_s). The matrices A, Y, and B
* are from iterative refinement (see zla_gerfsx_extended.f).
*
* BERR (output) COMPLEX*16 array, dimension (NRHS)
* The componentwise relative backward error from the formula above.
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION TMP
INTEGER I, J
COMPLEX*16 CDUM
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, REAL, DIMAG, MAX
* ..
* .. External Functions ..
EXTERNAL DLAMCH
DOUBLE PRECISION DLAMCH
DOUBLE PRECISION SAFE1
* ..
* .. Statement Functions ..
COMPLEX*16 CABS1
* ..
* .. Statement Function Definitions ..
CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
* ..
* .. Executable Statements ..
*
* Adding SAFE1 to the numerator guards against spuriously zero
* residuals. A similar safeguard is in the CLA_yyAMV routine used
* to compute AYB.
*
SAFE1 = DLAMCH( 'Safe minimum' )
SAFE1 = (NZ+1)*SAFE1
DO J = 1, NRHS
BERR(J) = 0.0D+0
DO I = 1, N
IF (AYB(I,J) .NE. 0.0D+0) THEN
TMP = (SAFE1 + CABS1(RES(I,J)))/AYB(I,J)
BERR(J) = MAX( BERR(J), TMP )
END IF
*
* If AYB is exactly 0.0 (and if computed by CLA_yyAMV), then we know
* the true residual also must be exactly 0.0.
*
END DO
END DO
END
|
