LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
subroutine ctrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
CTRMV
Definition: ctrmv.f:147
subroutine cget01(M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK, RESID)
CGET01
Definition: cget01.f:108
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