LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
zlatdf.f
Go to the documentation of this file.
1 *> \brief \b ZLATDF 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/zlatdf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlatdf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlatdf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLATDF( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
22 * JPIV )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER IJOB, LDZ, N
26 * DOUBLE PRECISION RDSCAL, RDSUM
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IPIV( * ), JPIV( * )
30 * COMPLEX*16 RHS( * ), Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> ZLATDF 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 ZGETC2. 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 ZGETC2 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 ZGECON, 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 choosen 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 DOUBLE PRECISION array, dimension (LDZ, N)
74 *> On entry, the LU part of the factorization of the n-by-n
75 *> matrix Z computed by ZGETC2: 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 DOUBLE PRECISION 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 DOUBLE PRECISION
95 *> On entry, the sum of squares of computed contributions to
96 *> the Dif-estimate under computation by ZTGSYL, 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 ZTGSY2 is called by CTGSYL.
102 *> \endverbatim
103 *>
104 *> \param[in,out] RDSCAL
105 *> \verbatim
106 *> RDSCAL is DOUBLE PRECISION
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 ZTGSY2 is called by
112 *> ZTGSYL.
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 September 2012
138 *
139 *> \ingroup complex16OTHERauxiliary
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 *>\n
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 zlatdf( IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV,
170  \$ jpiv )
171 *
172 * -- LAPACK auxiliary routine (version 3.4.2) --
173 * -- LAPACK is a software package provided by Univ. of Tennessee, --
174 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
175 * September 2012
176 *
177 * .. Scalar Arguments ..
178  INTEGER ijob, ldz, n
179  DOUBLE PRECISION rdscal, rdsum
180 * ..
181 * .. Array Arguments ..
182  INTEGER ipiv( * ), jpiv( * )
183  COMPLEX*16 rhs( * ), z( ldz, * )
184 * ..
185 *
186 * =====================================================================
187 *
188 * .. Parameters ..
189  INTEGER maxdim
190  parameter( maxdim = 2 )
191  DOUBLE PRECISION zero, one
192  parameter( zero = 0.0d+0, one = 1.0d+0 )
193  COMPLEX*16 cone
194  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
195 * ..
196 * .. Local Scalars ..
197  INTEGER i, info, j, k
198  DOUBLE PRECISION rtemp, scale, sminu, splus
199  COMPLEX*16 bm, bp, pmone, temp
200 * ..
201 * .. Local Arrays ..
202  DOUBLE PRECISION rwork( maxdim )
203  COMPLEX*16 work( 4*maxdim ), xm( maxdim ), xp( maxdim )
204 * ..
205 * .. External Subroutines ..
206  EXTERNAL zaxpy, zcopy, zgecon, zgesc2, zlassq, zlaswp,
207  \$ zscal
208 * ..
209 * .. External Functions ..
210  DOUBLE PRECISION dzasum
211  COMPLEX*16 zdotc
212  EXTERNAL dzasum, zdotc
213 * ..
214 * .. Intrinsic Functions ..
215  INTRINSIC abs, dble, sqrt
216 * ..
217 * .. Executable Statements ..
218 *
219  IF( ijob.NE.2 ) THEN
220 *
221 * Apply permutations IPIV to RHS
222 *
223  CALL zlaswp( 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 + dble( zdotc( n-j, z( j+1, j ), 1, z( j+1,
237  \$ j ), 1 ) )
238  sminu = dble( zdotc( n-j, z( j+1, j ), 1, rhs( j+1 ), 1 ) )
239  splus = splus*dble( rhs( j ) )
240  IF( splus.GT.sminu ) THEN
241  rhs( j ) = bp
242  ELSE IF( sminu.GT.splus ) THEN
243  rhs( j ) = bm
244  ELSE
245 *
246 * In this case the updating sums are equal and we can
247 * choose RHS(J) +1 or -1. The first time this happens we
248 * choose -1, thereafter +1. This is a simple way to get
249 * good estimates of matrices like Byers well-known example
250 * (see [1]). (Not done in BSOLVE.)
251 *
252  rhs( j ) = rhs( j ) + pmone
253  pmone = cone
254  END IF
255 *
256 * Compute the remaining r.h.s.
257 *
258  temp = -rhs( j )
259  CALL zaxpy( n-j, temp, z( j+1, j ), 1, rhs( j+1 ), 1 )
260  10 continue
261 *
262 * Solve for U- part, lockahead for RHS(N) = +-1. This is not done
263 * In BSOLVE and will hopefully give us a better estimate because
264 * any ill-conditioning of the original matrix is transfered to U
265 * and not to L. U(N, N) is an approximation to sigma_min(LU).
266 *
267  CALL zcopy( n-1, rhs, 1, work, 1 )
268  work( n ) = rhs( n ) + cone
269  rhs( n ) = rhs( n ) - cone
270  splus = zero
271  sminu = zero
272  DO 30 i = n, 1, -1
273  temp = cone / z( i, i )
274  work( i ) = work( i )*temp
275  rhs( i ) = rhs( i )*temp
276  DO 20 k = i + 1, n
277  work( i ) = work( i ) - work( k )*( z( i, k )*temp )
278  rhs( i ) = rhs( i ) - rhs( k )*( z( i, k )*temp )
279  20 continue
280  splus = splus + abs( work( i ) )
281  sminu = sminu + abs( rhs( i ) )
282  30 continue
283  IF( splus.GT.sminu )
284  \$ CALL zcopy( n, work, 1, rhs, 1 )
285 *
286 * Apply the permutations JPIV to the computed solution (RHS)
287 *
288  CALL zlaswp( 1, rhs, ldz, 1, n-1, jpiv, -1 )
289 *
290 * Compute the sum of squares
291 *
292  CALL zlassq( n, rhs, 1, rdscal, rdsum )
293  return
294  END IF
295 *
296 * ENTRY IJOB = 2
297 *
298 * Compute approximate nullvector XM of Z
299 *
300  CALL zgecon( 'I', n, z, ldz, one, rtemp, work, rwork, info )
301  CALL zcopy( n, work( n+1 ), 1, xm, 1 )
302 *
303 * Compute RHS
304 *
305  CALL zlaswp( 1, xm, ldz, 1, n-1, ipiv, -1 )
306  temp = cone / sqrt( zdotc( n, xm, 1, xm, 1 ) )
307  CALL zscal( n, temp, xm, 1 )
308  CALL zcopy( n, xm, 1, xp, 1 )
309  CALL zaxpy( n, cone, rhs, 1, xp, 1 )
310  CALL zaxpy( n, -cone, xm, 1, rhs, 1 )
311  CALL zgesc2( n, z, ldz, rhs, ipiv, jpiv, scale )
312  CALL zgesc2( n, z, ldz, xp, ipiv, jpiv, scale )
313  IF( dzasum( n, xp, 1 ).GT.dzasum( n, rhs, 1 ) )
314  \$ CALL zcopy( n, xp, 1, rhs, 1 )
315 *
316 * Compute the sum of squares
317 *
318  CALL zlassq( n, rhs, 1, rdscal, rdsum )
319  return
320 *
321 * End of ZLATDF
322 *
323  END