LAPACK  3.6.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
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 *> \date June 2016
138 *
139 *> \ingroup complexOTHERauxiliary
140 *
141 *> \par Further Details:
142 * =====================
143 *>
144 *> This routine is a further developed implementation of algorithm
145 *> BSOLVE in [1] using complete pivoting in the LU factorization.
146 *
147 *> \par Contributors:
148 * ==================
149 *>
150 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
151 *> Umea University, S-901 87 Umea, Sweden.
152 *
153 *> \par References:
154 * ================
155 *>
156 *> [1] Bo Kagstrom and Lars Westin,
157 *> Generalized Schur Methods with Condition Estimators for
158 *> Solving the Generalized Sylvester Equation, IEEE Transactions
159 *> on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751.
160 *>
161 *> [2] Peter Poromaa,
162 *> On Efficient and Robust Estimators for the Separation
163 *> between two Regular Matrix Pairs with Applications in
164 *> Condition Estimation. Report UMINF-95.05, Department of
165 *> Computing Science, Umea University, S-901 87 Umea, Sweden,
166 *> 1995.
167 *
168 * =====================================================================
169  SUBROUTINE clatdf( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
170  \$ jpiv )
171 *
172 * -- LAPACK auxiliary routine (version 3.6.1) --
173 * -- LAPACK is a software package provided by Univ. of Tennessee, --
174 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
175 * June 2016
176 *
177 * .. Scalar Arguments ..
178  INTEGER IJOB, LDZ, N
179  REAL RDSCAL, RDSUM
180 * ..
181 * .. Array Arguments ..
182  INTEGER IPIV( * ), JPIV( * )
183  COMPLEX RHS( * ), Z( ldz, * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189  INTEGER MAXDIM
190  parameter ( maxdim = 2 )
191  REAL ZERO, ONE
192  parameter ( zero = 0.0e+0, one = 1.0e+0 )
193  COMPLEX CONE
194  parameter ( cone = ( 1.0e+0, 0.0e+0 ) )
195 * ..
196 * .. Local Scalars ..
197  INTEGER I, INFO, J, K
198  REAL RTEMP, SCALE, SMINU, SPLUS
199  COMPLEX BM, BP, PMONE, TEMP
200 * ..
201 * .. Local Arrays ..
202  REAL RWORK( maxdim )
203  COMPLEX WORK( 4*maxdim ), XM( maxdim ), XP( maxdim )
204 * ..
205 * .. External Subroutines ..
206  EXTERNAL caxpy, ccopy, cgecon, cgesc2, classq, claswp,
207  \$ cscal
208 * ..
209 * .. External Functions ..
210  REAL SCASUM
211  COMPLEX CDOTC
212  EXTERNAL scasum, cdotc
213 * ..
214 * .. Intrinsic Functions ..
215  INTRINSIC abs, REAL, SQRT
216 * ..
217 * .. Executable Statements ..
218 *
219  IF( ijob.NE.2 ) THEN
220 *
221 * Apply permutations IPIV to RHS
222 *
223  CALL claswp( 1, rhs, ldz, 1, n-1, ipiv, 1 )
224 *
225 * Solve for L-part choosing RHS either to +1 or -1.
226 *
227  pmone = -cone
228  DO 10 j = 1, n - 1
229  bp = rhs( j ) + cone
230  bm = rhs( j ) - cone
231  splus = one
232 *
233 * Lockahead for L- part RHS(1:N-1) = +-1
234 * SPLUS and SMIN computed more efficiently than in BSOLVE[1].
235 *
236  splus = splus + REAL( CDOTC( N-J, Z( J+1, J ), 1, Z( J+1, \$ J ), 1 ) )
237  sminu = REAL( CDOTC( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 ) )
238  splus = splus*REAL( RHS( J ) )
239  IF( splus.GT.sminu ) THEN
240  rhs( j ) = bp
241  ELSE IF( sminu.GT.splus ) THEN
242  rhs( j ) = bm
243  ELSE
244 *
245 * In this case the updating sums are equal and we can
246 * choose RHS(J) +1 or -1. The first time this happens we
247 * choose -1, thereafter +1. This is a simple way to get
248 * good estimates of matrices like Byers well-known example
249 * (see [1]). (Not done in BSOLVE.)
250 *
251  rhs( j ) = rhs( j ) + pmone
252  pmone = cone
253  END IF
254 *
255 * Compute the remaining r.h.s.
256 *
257  temp = -rhs( j )
258  CALL caxpy( n-j, temp, z( j+1, j ), 1, rhs( j+1 ), 1 )
259  10 CONTINUE
260 *
261 * Solve for U- part, lockahead for RHS(N) = +-1. This is not done
262 * In BSOLVE and will hopefully give us a better estimate because
263 * any ill-conditioning of the original matrix is transfered to U
264 * and not to L. U(N, N) is an approximation to sigma_min(LU).
265 *
266  CALL ccopy( n-1, rhs, 1, work, 1 )
267  work( n ) = rhs( n ) + cone
268  rhs( n ) = rhs( n ) - cone
269  splus = zero
270  sminu = zero
271  DO 30 i = n, 1, -1
272  temp = cone / z( i, i )
273  work( i ) = work( i )*temp
274  rhs( i ) = rhs( i )*temp
275  DO 20 k = i + 1, n
276  work( i ) = work( i ) - work( k )*( z( i, k )*temp )
277  rhs( i ) = rhs( i ) - rhs( k )*( z( i, k )*temp )
278  20 CONTINUE
279  splus = splus + abs( work( i ) )
280  sminu = sminu + abs( rhs( i ) )
281  30 CONTINUE
282  IF( splus.GT.sminu )
283  \$ CALL ccopy( n, work, 1, rhs, 1 )
284 *
285 * Apply the permutations JPIV to the computed solution (RHS)
286 *
287  CALL claswp( 1, rhs, ldz, 1, n-1, jpiv, -1 )
288 *
289 * Compute the sum of squares
290 *
291  CALL classq( n, rhs, 1, rdscal, rdsum )
292  RETURN
293  END IF
294 *
295 * ENTRY IJOB = 2
296 *
297 * Compute approximate nullvector XM of Z
298 *
299  CALL cgecon( 'I', n, z, ldz, one, rtemp, work, rwork, info )
300  CALL ccopy( n, work( n+1 ), 1, xm, 1 )
301 *
302 * Compute RHS
303 *
304  CALL claswp( 1, xm, ldz, 1, n-1, ipiv, -1 )
305  temp = cone / sqrt( cdotc( n, xm, 1, xm, 1 ) )
306  CALL cscal( n, temp, xm, 1 )
307  CALL ccopy( n, xm, 1, xp, 1 )
308  CALL caxpy( n, cone, rhs, 1, xp, 1 )
309  CALL caxpy( n, -cone, xm, 1, rhs, 1 )
310  CALL cgesc2( n, z, ldz, rhs, ipiv, jpiv, scale )
311  CALL cgesc2( n, z, ldz, xp, ipiv, jpiv, scale )
312  IF( scasum( n, xp, 1 ).GT.scasum( n, rhs, 1 ) )
313  \$ CALL ccopy( n, xp, 1, rhs, 1 )
314 *
315 * Compute the sum of squares
316 *
317  CALL classq( n, rhs, 1, rdscal, rdsum )
318  RETURN
319 *
320 * End of CLATDF
321 *
322  END
323
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:171
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f:108
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
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:117
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine cgecon(NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK, INFO)
CGECON
Definition: cgecon.f:126
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
subroutine claswp(N, A, LDA, K1, K2, IPIV, INCX)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: claswp.f:116