LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dget03.f
Go to the documentation of this file.
1*> \brief \b DGET03
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DGET03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK,
12* RCOND, RESID )
13*
14* .. Scalar Arguments ..
15* INTEGER LDA, LDAINV, LDWORK, N
16* DOUBLE PRECISION RCOND, RESID
17* ..
18* .. Array Arguments ..
19* DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
20* $ WORK( LDWORK, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DGET03 computes the residual for a general matrix times its inverse:
30*> norm( I - AINV*A ) / ( N * norm(A) * norm(AINV) * EPS ),
31*> where EPS is the machine epsilon.
32*> \endverbatim
33*
34* Arguments:
35* ==========
36*
37*> \param[in] N
38*> \verbatim
39*> N is INTEGER
40*> The number of rows and columns of the matrix A. N >= 0.
41*> \endverbatim
42*>
43*> \param[in] A
44*> \verbatim
45*> A is DOUBLE PRECISION array, dimension (LDA,N)
46*> The original N x N matrix A.
47*> \endverbatim
48*>
49*> \param[in] LDA
50*> \verbatim
51*> LDA is INTEGER
52*> The leading dimension of the array A. LDA >= max(1,N).
53*> \endverbatim
54*>
55*> \param[in] AINV
56*> \verbatim
57*> AINV is DOUBLE PRECISION array, dimension (LDAINV,N)
58*> The inverse of the matrix A.
59*> \endverbatim
60*>
61*> \param[in] LDAINV
62*> \verbatim
63*> LDAINV is INTEGER
64*> The leading dimension of the array AINV. LDAINV >= max(1,N).
65*> \endverbatim
66*>
67*> \param[out] WORK
68*> \verbatim
69*> WORK is DOUBLE PRECISION array, dimension (LDWORK,N)
70*> \endverbatim
71*>
72*> \param[in] LDWORK
73*> \verbatim
74*> LDWORK is INTEGER
75*> The leading dimension of the array WORK. LDWORK >= max(1,N).
76*> \endverbatim
77*>
78*> \param[out] RWORK
79*> \verbatim
80*> RWORK is DOUBLE PRECISION array, dimension (N)
81*> \endverbatim
82*>
83*> \param[out] RCOND
84*> \verbatim
85*> RCOND is DOUBLE PRECISION
86*> The reciprocal of the condition number of A, computed as
87*> ( 1/norm(A) ) / norm(AINV).
88*> \endverbatim
89*>
90*> \param[out] RESID
91*> \verbatim
92*> RESID is DOUBLE PRECISION
93*> norm(I - AINV*A) / ( N * norm(A) * norm(AINV) * EPS )
94*> \endverbatim
95*
96* Authors:
97* ========
98*
99*> \author Univ. of Tennessee
100*> \author Univ. of California Berkeley
101*> \author Univ. of Colorado Denver
102*> \author NAG Ltd.
103*
104*> \ingroup double_lin
105*
106* =====================================================================
107 SUBROUTINE dget03( N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK,
108 $ RCOND, RESID )
109*
110* -- LAPACK test routine --
111* -- LAPACK is a software package provided by Univ. of Tennessee, --
112* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
113*
114* .. Scalar Arguments ..
115 INTEGER LDA, LDAINV, LDWORK, N
116 DOUBLE PRECISION RCOND, RESID
117* ..
118* .. Array Arguments ..
119 DOUBLE PRECISION A( LDA, * ), AINV( LDAINV, * ), RWORK( * ),
120 $ work( ldwork, * )
121* ..
122*
123* =====================================================================
124*
125* .. Parameters ..
126 DOUBLE PRECISION ZERO, ONE
127 parameter( zero = 0.0d+0, one = 1.0d+0 )
128* ..
129* .. Local Scalars ..
130 INTEGER I
131 DOUBLE PRECISION AINVNM, ANORM, EPS
132* ..
133* .. External Functions ..
134 DOUBLE PRECISION DLAMCH, DLANGE
135 EXTERNAL dlamch, dlange
136* ..
137* .. External Subroutines ..
138 EXTERNAL dgemm
139* ..
140* .. Intrinsic Functions ..
141 INTRINSIC dble
142* ..
143* .. Executable Statements ..
144*
145* Quick exit if N = 0.
146*
147 IF( n.LE.0 ) THEN
148 rcond = one
149 resid = zero
150 RETURN
151 END IF
152*
153* Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
154*
155 eps = dlamch( 'Epsilon' )
156 anorm = dlange( '1', n, n, a, lda, rwork )
157 ainvnm = dlange( '1', n, n, ainv, ldainv, rwork )
158 IF( anorm.LE.zero .OR. ainvnm.LE.zero ) THEN
159 rcond = zero
160 resid = one / eps
161 RETURN
162 END IF
163 rcond = ( one / anorm ) / ainvnm
164*
165* Compute I - A * AINV
166*
167 CALL dgemm( 'No transpose', 'No transpose', n, n, n, -one, ainv,
168 $ ldainv, a, lda, zero, work, ldwork )
169 DO 10 i = 1, n
170 work( i, i ) = one + work( i, i )
171 10 CONTINUE
172*
173* Compute norm(I - AINV*A) / (N * norm(A) * norm(AINV) * EPS)
174*
175 resid = dlange( '1', n, n, work, ldwork, rwork )
176*
177 resid = ( ( resid*rcond ) / eps ) / dble( n )
178*
179 RETURN
180*
181* End of DGET03
182*
183 END
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dget03(N, A, LDA, AINV, LDAINV, WORK, LDWORK, RWORK, RCOND, RESID)
DGET03
Definition: dget03.f:109