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