LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dget01.f
Go to the documentation of this file.
1*> \brief \b DGET01
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 DGET01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
12* RESID )
13*
14* .. Scalar Arguments ..
15* INTEGER LDA, LDAFAC, M, N
16* DOUBLE PRECISION RESID
17* ..
18* .. Array Arguments ..
19* INTEGER IPIV( * )
20* DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DGET01 reconstructs a matrix A from its L*U factorization and
30*> computes the residual
31*> norm(L*U - A) / ( N * norm(A) * EPS ),
32*> where EPS is the machine epsilon.
33*> \endverbatim
34*
35* Arguments:
36* ==========
37*
38*> \param[in] M
39*> \verbatim
40*> M is INTEGER
41*> The number of rows of the matrix A. M >= 0.
42*> \endverbatim
43*>
44*> \param[in] N
45*> \verbatim
46*> N is INTEGER
47*> The number of columns of the matrix A. N >= 0.
48*> \endverbatim
49*>
50*> \param[in] A
51*> \verbatim
52*> A is DOUBLE PRECISION array, dimension (LDA,N)
53*> The original M x N matrix A.
54*> \endverbatim
55*>
56*> \param[in] LDA
57*> \verbatim
58*> LDA is INTEGER
59*> The leading dimension of the array A. LDA >= max(1,M).
60*> \endverbatim
61*>
62*> \param[in,out] AFAC
63*> \verbatim
64*> AFAC is DOUBLE PRECISION array, dimension (LDAFAC,N)
65*> The factored form of the matrix A. AFAC contains the factors
66*> L and U from the L*U factorization as computed by DGETRF.
67*> Overwritten with the reconstructed matrix, and then with the
68*> difference L*U - A.
69*> \endverbatim
70*>
71*> \param[in] LDAFAC
72*> \verbatim
73*> LDAFAC is INTEGER
74*> The leading dimension of the array AFAC. LDAFAC >= max(1,M).
75*> \endverbatim
76*>
77*> \param[in] IPIV
78*> \verbatim
79*> IPIV is INTEGER array, dimension (N)
80*> The pivot indices from DGETRF.
81*> \endverbatim
82*>
83*> \param[out] RWORK
84*> \verbatim
85*> RWORK is DOUBLE PRECISION array, dimension (M)
86*> \endverbatim
87*>
88*> \param[out] RESID
89*> \verbatim
90*> RESID is DOUBLE PRECISION
91*> norm(L*U - A) / ( N * norm(A) * EPS )
92*> \endverbatim
93*
94* Authors:
95* ========
96*
97*> \author Univ. of Tennessee
98*> \author Univ. of California Berkeley
99*> \author Univ. of Colorado Denver
100*> \author NAG Ltd.
101*
102*> \ingroup double_lin
103*
104* =====================================================================
105 SUBROUTINE dget01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
106 $ RESID )
107*
108* -- LAPACK test routine --
109* -- LAPACK is a software package provided by Univ. of Tennessee, --
110* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
111*
112* .. Scalar Arguments ..
113 INTEGER LDA, LDAFAC, M, N
114 DOUBLE PRECISION RESID
115* ..
116* .. Array Arguments ..
117 INTEGER IPIV( * )
118 DOUBLE PRECISION A( LDA, * ), AFAC( LDAFAC, * ), RWORK( * )
119* ..
120*
121* =====================================================================
122*
123*
124* .. Parameters ..
125 DOUBLE PRECISION ZERO, ONE
126 parameter( zero = 0.0d+0, one = 1.0d+0 )
127* ..
128* .. Local Scalars ..
129 INTEGER I, J, K
130 DOUBLE PRECISION ANORM, EPS, T
131* ..
132* .. External Functions ..
133 DOUBLE PRECISION DDOT, DLAMCH, DLANGE
134 EXTERNAL ddot, dlamch, dlange
135* ..
136* .. External Subroutines ..
137 EXTERNAL dgemv, dlaswp, dscal, dtrmv
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC dble, min
141* ..
142* .. Executable Statements ..
143*
144* Quick exit if M = 0 or N = 0.
145*
146 IF( m.LE.0 .OR. n.LE.0 ) THEN
147 resid = zero
148 RETURN
149 END IF
150*
151* Determine EPS and the norm of A.
152*
153 eps = dlamch( 'Epsilon' )
154 anorm = dlange( '1', m, n, a, lda, rwork )
155*
156* Compute the product L*U and overwrite AFAC with the result.
157* A column at a time of the product is obtained, starting with
158* column N.
159*
160 DO 10 k = n, 1, -1
161 IF( k.GT.m ) THEN
162 CALL dtrmv( 'Lower', 'No transpose', 'Unit', m, afac,
163 $ ldafac, afac( 1, k ), 1 )
164 ELSE
165*
166* Compute elements (K+1:M,K)
167*
168 t = afac( k, k )
169 IF( k+1.LE.m ) THEN
170 CALL dscal( m-k, t, afac( k+1, k ), 1 )
171 CALL dgemv( 'No transpose', m-k, k-1, one,
172 $ afac( k+1, 1 ), ldafac, afac( 1, k ), 1, one,
173 $ afac( k+1, k ), 1 )
174 END IF
175*
176* Compute the (K,K) element
177*
178 afac( k, k ) = t + ddot( k-1, afac( k, 1 ), ldafac,
179 $ afac( 1, k ), 1 )
180*
181* Compute elements (1:K-1,K)
182*
183 CALL dtrmv( 'Lower', 'No transpose', 'Unit', k-1, afac,
184 $ ldafac, afac( 1, k ), 1 )
185 END IF
186 10 CONTINUE
187 CALL dlaswp( n, afac, ldafac, 1, min( m, n ), ipiv, -1 )
188*
189* Compute the difference L*U - A and store in AFAC.
190*
191 DO 30 j = 1, n
192 DO 20 i = 1, m
193 afac( i, j ) = afac( i, j ) - a( i, j )
194 20 CONTINUE
195 30 CONTINUE
196*
197* Compute norm( L*U - A ) / ( N * norm(A) * EPS )
198*
199 resid = dlange( '1', m, n, afac, ldafac, rwork )
200*
201 IF( anorm.LE.zero ) THEN
202 IF( resid.NE.zero )
203 $ resid = one / eps
204 ELSE
205 resid = ( ( resid / dble( n ) ) / anorm ) / eps
206 END IF
207*
208 RETURN
209*
210* End of DGET01
211*
212 END
subroutine dget01(m, n, a, lda, afac, ldafac, ipiv, rwork, resid)
DGET01
Definition dget01.f:107
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dlaswp(n, a, lda, k1, k2, ipiv, incx)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition dlaswp.f:115
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dtrmv(uplo, trans, diag, n, a, lda, x, incx)
DTRMV
Definition dtrmv.f:147