LAPACK 3.11.0
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*> \htmlonly
9*> Download CLATDF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clatdf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clatdf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clatdf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
22* JPIV )
23*
24* .. Scalar Arguments ..
25* INTEGER IJOB, LDZ, N
26* REAL RDSCAL, RDSUM
27* ..
28* .. Array Arguments ..
29* INTEGER IPIV( * ), JPIV( * )
30* COMPLEX RHS( * ), Z( LDZ, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CLATDF computes the contribution to the reciprocal Dif-estimate
40*> by solving for x in Z * x = b, where b is chosen such that the norm
41*> of x is as large as possible. It is assumed that LU decomposition
42*> of Z has been computed by CGETC2. On entry RHS = f holds the
43*> contribution from earlier solved sub-systems, and on return RHS = x.
44*>
45*> The factorization of Z returned by CGETC2 has the form
46*> Z = P * L * U * Q, where P and Q are permutation matrices. L is lower
47*> triangular with unit diagonal elements and U is upper triangular.
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] IJOB
54*> \verbatim
55*> IJOB is INTEGER
56*> IJOB = 2: First compute an approximative null-vector e
57*> of Z using CGECON, e is normalized and solve for
58*> Zx = +-e - f with the sign giving the greater value of
59*> 2-norm(x). About 5 times as expensive as Default.
60*> IJOB .ne. 2: Local look ahead strategy where
61*> all entries of the r.h.s. b is chosen as either +1 or
62*> -1. Default.
63*> \endverbatim
64*>
65*> \param[in] N
66*> \verbatim
67*> N is INTEGER
68*> The number of columns of the matrix Z.
69*> \endverbatim
70*>
71*> \param[in] Z
72*> \verbatim
73*> Z is COMPLEX array, dimension (LDZ, N)
74*> On entry, the LU part of the factorization of the n-by-n
75*> matrix Z computed by CGETC2: Z = P * L * U * Q
76*> \endverbatim
77*>
78*> \param[in] LDZ
79*> \verbatim
80*> LDZ is INTEGER
81*> The leading dimension of the array Z. LDA >= max(1, N).
82*> \endverbatim
83*>
84*> \param[in,out] RHS
85*> \verbatim
86*> RHS is COMPLEX array, dimension (N).
87*> On entry, RHS contains contributions from other subsystems.
88*> On exit, RHS contains the solution of the subsystem with
89*> entries according to the value of IJOB (see above).
90*> \endverbatim
91*>
92*> \param[in,out] RDSUM
93*> \verbatim
94*> RDSUM is REAL
95*> On entry, the sum of squares of computed contributions to
96*> the Dif-estimate under computation by CTGSYL, where the
97*> scaling factor RDSCAL (see below) has been factored out.
98*> On exit, the corresponding sum of squares updated with the
99*> contributions from the current sub-system.
100*> If TRANS = 'T' RDSUM is not touched.
101*> NOTE: RDSUM only makes sense when CTGSY2 is called by CTGSYL.
102*> \endverbatim
103*>
104*> \param[in,out] RDSCAL
105*> \verbatim
106*> RDSCAL is REAL
107*> On entry, scaling factor used to prevent overflow in RDSUM.
108*> On exit, RDSCAL is updated w.r.t. the current contributions
109*> in RDSUM.
110*> If TRANS = 'T', RDSCAL is not touched.
111*> NOTE: RDSCAL only makes sense when CTGSY2 is called by
112*> CTGSYL.
113*> \endverbatim
114*>
115*> \param[in] IPIV
116*> \verbatim
117*> IPIV is INTEGER array, dimension (N).
118*> The pivot indices; for 1 <= i <= N, row i of the
119*> matrix has been interchanged with row IPIV(i).
120*> \endverbatim
121*>
122*> \param[in] JPIV
123*> \verbatim
124*> JPIV is INTEGER array, dimension (N).
125*> The pivot indices; for 1 <= j <= N, column j of the
126*> matrix has been interchanged with column JPIV(j).
127*> \endverbatim
128*
129* Authors:
130* ========
131*
132*> \author Univ. of Tennessee
133*> \author Univ. of California Berkeley
134*> \author Univ. of Colorado Denver
135*> \author NAG Ltd.
136*
137*> \ingroup complexOTHERauxiliary
138*
139*> \par Further Details:
140* =====================
141*>
142*> This routine is a further developed implementation of algorithm
143*> BSOLVE in [1] using complete pivoting in the LU factorization.
144*
145*> \par Contributors:
146* ==================
147*>
148*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
149*> Umea University, S-901 87 Umea, Sweden.
150*
151*> \par References:
152* ================
153*>
154*> [1] Bo Kagstrom and Lars Westin,
155*> Generalized Schur Methods with Condition Estimators for
156*> Solving the Generalized Sylvester Equation, IEEE Transactions
157*> on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
158*>
159*> [2] Peter Poromaa,
160*> On Efficient and Robust Estimators for the Separation
161*> between two Regular Matrix Pairs with Applications in
162*> Condition Estimation. Report UMINF-95.05, Department of
163*> Computing Science, Umea University, S-901 87 Umea, Sweden,
164*> 1995.
165*
166* =====================================================================
167 SUBROUTINE clatdf( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
168 $ JPIV )
169*
170* -- LAPACK auxiliary routine --
171* -- LAPACK is a software package provided by Univ. of Tennessee, --
172* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173*
174* .. Scalar Arguments ..
175 INTEGER IJOB, LDZ, N
176 REAL RDSCAL, RDSUM
177* ..
178* .. Array Arguments ..
179 INTEGER IPIV( * ), JPIV( * )
180 COMPLEX RHS( * ), Z( LDZ, * )
181* ..
182*
183* =====================================================================
184*
185* .. Parameters ..
186 INTEGER MAXDIM
187 parameter( maxdim = 2 )
188 REAL ZERO, ONE
189 parameter( zero = 0.0e+0, one = 1.0e+0 )
190 COMPLEX CONE
191 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
192* ..
193* .. Local Scalars ..
194 INTEGER I, INFO, J, K
195 REAL RTEMP, SCALE, SMINU, SPLUS
196 COMPLEX BM, BP, PMONE, TEMP
197* ..
198* .. Local Arrays ..
199 REAL RWORK( MAXDIM )
200 COMPLEX WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
201* ..
202* .. External Subroutines ..
203 EXTERNAL caxpy, ccopy, cgecon, cgesc2, classq, claswp,
204 $ cscal
205* ..
206* .. External Functions ..
207 REAL SCASUM
208 COMPLEX CDOTC
209 EXTERNAL scasum, cdotc
210* ..
211* .. Intrinsic Functions ..
212 INTRINSIC abs, real, sqrt
213* ..
214* .. Executable Statements ..
215*
216 IF( ijob.NE.2 ) THEN
217*
218* Apply permutations IPIV to RHS
219*
220 CALL claswp( 1, rhs, ldz, 1, n-1, ipiv, 1 )
221*
222* Solve for L-part choosing RHS either to +1 or -1.
223*
224 pmone = -cone
225 DO 10 j = 1, n - 1
226 bp = rhs( j ) + cone
227 bm = rhs( j ) - cone
228 splus = one
229*
230* Lockahead for L- part RHS(1:N-1) = +-1
231* SPLUS and SMIN computed more efficiently than in BSOLVE[1].
232*
233 splus = splus + real( cdotc( n-j, z( j+1, j ), 1, z( j+1,
234 $ j ), 1 ) )
235 sminu = real( cdotc( n-j, z( j+1, j ), 1, rhs( j+1 ), 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 classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f90:137
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
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:115
subroutine cgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
CGECON
Definition: cgecon.f:124
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 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:169