LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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