LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clatdf.f
Go to the documentation of this file.
1*> \brief \b CLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution to the reciprocal Dif-estimate.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CLATDF + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatdf.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatdf.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatdf.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
20* JPIV )
21*
22* .. Scalar Arguments ..
23* INTEGER IJOB, LDZ, N
24* REAL RDSCAL, RDSUM
25* ..
26* .. Array Arguments ..
27* INTEGER IPIV( * ), JPIV( * )
28* COMPLEX RHS( * ), Z( LDZ, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> CLATDF computes the contribution to the reciprocal Dif-estimate
38*> by solving for x in Z * x = b, where b is chosen such that the norm
39*> of x is as large as possible. It is assumed that LU decomposition
40*> of Z has been computed by CGETC2. On entry RHS = f holds the
41*> contribution from earlier solved sub-systems, and on return RHS = x.
42*>
43*> The factorization of Z returned by CGETC2 has the form
44*> Z = P * L * U * Q, where P and Q are permutation matrices. L is lower
45*> triangular with unit diagonal elements and U is upper triangular.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] IJOB
52*> \verbatim
53*> IJOB is INTEGER
54*> IJOB = 2: First compute an approximative null-vector e
55*> of Z using CGECON, e is normalized and solve for
56*> Zx = +-e - f with the sign giving the greater value of
57*> 2-norm(x). About 5 times as expensive as Default.
58*> IJOB .ne. 2: Local look ahead strategy where
59*> all entries of the r.h.s. b is chosen as either +1 or
60*> -1. Default.
61*> \endverbatim
62*>
63*> \param[in] N
64*> \verbatim
65*> N is INTEGER
66*> The number of columns of the matrix Z.
67*> \endverbatim
68*>
69*> \param[in] Z
70*> \verbatim
71*> Z is COMPLEX array, dimension (LDZ, N)
72*> On entry, the LU part of the factorization of the n-by-n
73*> matrix Z computed by CGETC2: Z = P * L * U * Q
74*> \endverbatim
75*>
76*> \param[in] LDZ
77*> \verbatim
78*> LDZ is INTEGER
79*> The leading dimension of the array Z. LDA >= max(1, N).
80*> \endverbatim
81*>
82*> \param[in,out] RHS
83*> \verbatim
84*> RHS is COMPLEX array, dimension (N).
85*> On entry, RHS contains contributions from other subsystems.
86*> On exit, RHS contains the solution of the subsystem with
87*> entries according to the value of IJOB (see above).
88*> \endverbatim
89*>
90*> \param[in,out] RDSUM
91*> \verbatim
92*> RDSUM is REAL
93*> On entry, the sum of squares of computed contributions to
94*> the Dif-estimate under computation by CTGSYL, where the
95*> scaling factor RDSCAL (see below) has been factored out.
96*> On exit, the corresponding sum of squares updated with the
97*> contributions from the current sub-system.
98*> If TRANS = 'T' RDSUM is not touched.
99*> NOTE: RDSUM only makes sense when CTGSY2 is called by CTGSYL.
100*> \endverbatim
101*>
102*> \param[in,out] RDSCAL
103*> \verbatim
104*> RDSCAL is REAL
105*> On entry, scaling factor used to prevent overflow in RDSUM.
106*> On exit, RDSCAL is updated w.r.t. the current contributions
107*> in RDSUM.
108*> If TRANS = 'T', RDSCAL is not touched.
109*> NOTE: RDSCAL only makes sense when CTGSY2 is called by
110*> CTGSYL.
111*> \endverbatim
112*>
113*> \param[in] IPIV
114*> \verbatim
115*> IPIV is INTEGER array, dimension (N).
116*> The pivot indices; for 1 <= i <= N, row i of the
117*> matrix has been interchanged with row IPIV(i).
118*> \endverbatim
119*>
120*> \param[in] JPIV
121*> \verbatim
122*> JPIV is INTEGER array, dimension (N).
123*> The pivot indices; for 1 <= j <= N, column j of the
124*> matrix has been interchanged with column JPIV(j).
125*> \endverbatim
126*
127* Authors:
128* ========
129*
130*> \author Univ. of Tennessee
131*> \author Univ. of California Berkeley
132*> \author Univ. of Colorado Denver
133*> \author NAG Ltd.
134*
135*> \ingroup latdf
136*
137*> \par Further Details:
138* =====================
139*>
140*> This routine is a further developed implementation of algorithm
141*> BSOLVE in [1] using complete pivoting in the LU factorization.
142*
143*> \par Contributors:
144* ==================
145*>
146*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
147*> Umea University, S-901 87 Umea, Sweden.
148*
149*> \par References:
150* ================
151*>
152*> [1] Bo Kagstrom and Lars Westin,
153*> Generalized Schur Methods with Condition Estimators for
154*> Solving the Generalized Sylvester Equation, IEEE Transactions
155*> on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
156*>
157*> [2] Peter Poromaa,
158*> On Efficient and Robust Estimators for the Separation
159*> between two Regular Matrix Pairs with Applications in
160*> Condition Estimation. Report UMINF-95.05, Department of
161*> Computing Science, Umea University, S-901 87 Umea, Sweden,
162*> 1995.
163*
164* =====================================================================
165 SUBROUTINE clatdf( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
166 $ JPIV )
167*
168* -- LAPACK auxiliary routine --
169* -- LAPACK is a software package provided by Univ. of Tennessee, --
170* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171*
172* .. Scalar Arguments ..
173 INTEGER IJOB, LDZ, N
174 REAL RDSCAL, RDSUM
175* ..
176* .. Array Arguments ..
177 INTEGER IPIV( * ), JPIV( * )
178 COMPLEX RHS( * ), Z( LDZ, * )
179* ..
180*
181* =====================================================================
182*
183* .. Parameters ..
184 INTEGER MAXDIM
185 parameter( maxdim = 2 )
186 REAL ZERO, ONE
187 parameter( zero = 0.0e+0, one = 1.0e+0 )
188 COMPLEX CONE
189 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
190* ..
191* .. Local Scalars ..
192 INTEGER I, INFO, J, K
193 REAL RTEMP, SCALE, SMINU, SPLUS
194 COMPLEX BM, BP, PMONE, TEMP
195* ..
196* .. Local Arrays ..
197 REAL RWORK( MAXDIM )
198 COMPLEX WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
199* ..
200* .. External Subroutines ..
201 EXTERNAL caxpy, ccopy, cgecon, cgesc2, classq,
202 $ claswp,
203 $ cscal
204* ..
205* .. External Functions ..
206 REAL SCASUM
207 COMPLEX CDOTC
208 EXTERNAL scasum, cdotc
209* ..
210* .. Intrinsic Functions ..
211 INTRINSIC abs, real, sqrt
212* ..
213* .. Executable Statements ..
214*
215 IF( ijob.NE.2 ) THEN
216*
217* Apply permutations IPIV to RHS
218*
219 CALL claswp( 1, rhs, ldz, 1, n-1, ipiv, 1 )
220*
221* Solve for L-part choosing RHS either to +1 or -1.
222*
223 pmone = -cone
224 DO 10 j = 1, n - 1
225 bp = rhs( j ) + cone
226 bm = rhs( j ) - cone
227 splus = one
228*
229* Look-ahead for L- part RHS(1:N-1) = +-1
230* SPLUS and SMIN computed more efficiently than in BSOLVE[1].
231*
232 splus = splus + real( cdotc( n-j, z( j+1, j ), 1, z( j+1,
233 $ j ), 1 ) )
234 sminu = real( cdotc( n-j, z( j+1, j ), 1, rhs( j+1 ),
235 $ 1 ) )
236 splus = splus*real( rhs( j ) )
237 IF( splus.GT.sminu ) THEN
238 rhs( j ) = bp
239 ELSE IF( sminu.GT.splus ) THEN
240 rhs( j ) = bm
241 ELSE
242*
243* In this case the updating sums are equal and we can
244* choose RHS(J) +1 or -1. The first time this happens we
245* choose -1, thereafter +1. This is a simple way to get
246* good estimates of matrices like Byers well-known example
247* (see [1]). (Not done in BSOLVE.)
248*
249 rhs( j ) = rhs( j ) + pmone
250 pmone = cone
251 END IF
252*
253* Compute the remaining r.h.s.
254*
255 temp = -rhs( j )
256 CALL caxpy( n-j, temp, z( j+1, j ), 1, rhs( j+1 ), 1 )
257 10 CONTINUE
258*
259* Solve for U- part, lockahead for RHS(N) = +-1. This is not done
260* In BSOLVE and will hopefully give us a better estimate because
261* any ill-conditioning of the original matrix is transferred to U
262* and not to L. U(N, N) is an approximation to sigma_min(LU).
263*
264 CALL ccopy( n-1, rhs, 1, work, 1 )
265 work( n ) = rhs( n ) + cone
266 rhs( n ) = rhs( n ) - cone
267 splus = zero
268 sminu = zero
269 DO 30 i = n, 1, -1
270 temp = cone / z( i, i )
271 work( i ) = work( i )*temp
272 rhs( i ) = rhs( i )*temp
273 DO 20 k = i + 1, n
274 work( i ) = work( i ) - work( k )*( z( i, k )*temp )
275 rhs( i ) = rhs( i ) - rhs( k )*( z( i, k )*temp )
276 20 CONTINUE
277 splus = splus + abs( work( i ) )
278 sminu = sminu + abs( rhs( i ) )
279 30 CONTINUE
280 IF( splus.GT.sminu )
281 $ CALL ccopy( n, work, 1, rhs, 1 )
282*
283* Apply the permutations JPIV to the computed solution (RHS)
284*
285 CALL claswp( 1, rhs, ldz, 1, n-1, jpiv, -1 )
286*
287* Compute the sum of squares
288*
289 CALL classq( n, rhs, 1, rdscal, rdsum )
290 RETURN
291 END IF
292*
293* ENTRY IJOB = 2
294*
295* Compute approximate nullvector XM of Z
296*
297 CALL cgecon( 'I', n, z, ldz, one, rtemp, work, rwork, info )
298 CALL ccopy( n, work( n+1 ), 1, xm, 1 )
299*
300* Compute RHS
301*
302 CALL claswp( 1, xm, ldz, 1, n-1, ipiv, -1 )
303 temp = cone / sqrt( cdotc( n, xm, 1, xm, 1 ) )
304 CALL cscal( n, temp, xm, 1 )
305 CALL ccopy( n, xm, 1, xp, 1 )
306 CALL caxpy( n, cone, rhs, 1, xp, 1 )
307 CALL caxpy( n, -cone, xm, 1, rhs, 1 )
308 CALL cgesc2( n, z, ldz, rhs, ipiv, jpiv, scale )
309 CALL cgesc2( n, z, ldz, xp, ipiv, jpiv, scale )
310 IF( scasum( n, xp, 1 ).GT.scasum( n, rhs, 1 ) )
311 $ CALL ccopy( n, xp, 1, rhs, 1 )
312*
313* Compute the sum of squares
314*
315 CALL classq( n, rhs, 1, rdscal, rdsum )
316 RETURN
317*
318* End of CLATDF
319*
320 END
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgecon(norm, n, a, lda, anorm, rcond, work, rwork, info)
CGECON
Definition cgecon.f:130
subroutine cgesc2(n, a, lda, rhs, ipiv, jpiv, scale)
CGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition cgesc2.f:113
subroutine classq(n, x, incx, scale, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:122
subroutine claswp(n, a, lda, k1, k2, ipiv, incx)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition claswp.f:113
subroutine clatdf(ijob, n, z, ldz, rhs, rdsum, rdscal, ipiv, jpiv)
CLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition clatdf.f:167
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78