1 /* ./src_f77/zlatdf.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static doublecomplex c_b1 = {1.,0.};
11 static integer c__1 = 1;
12 static integer c_n1 = -1;
13 static doublereal c_b24 = 1.;
14
zlatdf_(integer * ijob,integer * n,doublecomplex * z__,integer * ldz,doublecomplex * rhs,doublereal * rdsum,doublereal * rdscal,integer * ipiv,integer * jpiv)15 /* Subroutine */ int zlatdf_(integer *ijob, integer *n, doublecomplex *z__,
16 integer *ldz, doublecomplex *rhs, doublereal *rdsum, doublereal *
17 rdscal, integer *ipiv, integer *jpiv)
18 {
19 /* System generated locals */
20 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
21 doublecomplex z__1, z__2, z__3;
22
23 /* Builtin functions */
24 void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
25 double z_abs(doublecomplex *);
26 void z_sqrt(doublecomplex *, doublecomplex *);
27
28 /* Local variables */
29 static integer i__, j, k;
30 static doublecomplex bm, bp, xm[2], xp[2];
31 static integer info;
32 static doublecomplex temp, work[8];
33 static doublereal scale;
34 extern /* Subroutine */ int zscal_(integer *, doublecomplex *,
35 doublecomplex *, integer *);
36 static doublecomplex pmone;
37 extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *,
38 doublecomplex *, integer *, doublecomplex *, integer *);
39 static doublereal rtemp, sminu, rwork[2];
40 extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *,
41 doublecomplex *, integer *);
42 static doublereal splus;
43 extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *,
44 doublecomplex *, integer *, doublecomplex *, integer *), zgesc2_(
45 integer *, doublecomplex *, integer *, doublecomplex *, integer *,
46 integer *, doublereal *), zgecon_(char *, integer *,
47 doublecomplex *, integer *, doublereal *, doublereal *,
48 doublecomplex *, doublereal *, integer *, ftnlen);
49 extern doublereal dzasum_(integer *, doublecomplex *, integer *);
50 extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *,
51 doublereal *, doublereal *), zlaswp_(integer *, doublecomplex *,
52 integer *, integer *, integer *, integer *, integer *);
53
54
55 /* -- LAPACK auxiliary routine (version 3.0) -- */
56 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
57 /* Courant Institute, Argonne National Lab, and Rice University */
58 /* June 30, 1999 */
59
60 /* .. Scalar Arguments .. */
61 /* .. */
62 /* .. Array Arguments .. */
63 /* .. */
64
65 /* Purpose */
66 /* ======= */
67
68 /* ZLATDF computes the contribution to the reciprocal Dif-estimate */
69 /* by solving for x in Z * x = b, where b is chosen such that the norm */
70 /* of x is as large as possible. It is assumed that LU decomposition */
71 /* of Z has been computed by ZGETC2. On entry RHS = f holds the */
72 /* contribution from earlier solved sub-systems, and on return RHS = x. */
73
74 /* The factorization of Z returned by ZGETC2 has the form */
75 /* Z = P * L * U * Q, where P and Q are permutation matrices. L is lower */
76 /* triangular with unit diagonal elements and U is upper triangular. */
77
78 /* Arguments */
79 /* ========= */
80
81 /* IJOB (input) INTEGER */
82 /* IJOB = 2: First compute an approximative null-vector e */
83 /* of Z using ZGECON, e is normalized and solve for */
84 /* Zx = +-e - f with the sign giving the greater value of */
85 /* 2-norm(x). About 5 times as expensive as Default. */
86 /* IJOB .ne. 2: Local look ahead strategy where */
87 /* all entries of the r.h.s. b is choosen as either +1 or */
88 /* -1. Default. */
89
90 /* N (input) INTEGER */
91 /* The number of columns of the matrix Z. */
92
93 /* Z (input) DOUBLE PRECISION array, dimension (LDZ, N) */
94 /* On entry, the LU part of the factorization of the n-by-n */
95 /* matrix Z computed by ZGETC2: Z = P * L * U * Q */
96
97 /* LDZ (input) INTEGER */
98 /* The leading dimension of the array Z. LDA >= max(1, N). */
99
100 /* RHS (input/output) DOUBLE PRECISION array, dimension (N). */
101 /* On entry, RHS contains contributions from other subsystems. */
102 /* On exit, RHS contains the solution of the subsystem with */
103 /* entries according to the value of IJOB (see above). */
104
105 /* RDSUM (input/output) DOUBLE PRECISION */
106 /* On entry, the sum of squares of computed contributions to */
107 /* the Dif-estimate under computation by ZTGSYL, where the */
108 /* scaling factor RDSCAL (see below) has been factored out. */
109 /* On exit, the corresponding sum of squares updated with the */
110 /* contributions from the current sub-system. */
111 /* If TRANS = 'T' RDSUM is not touched. */
112 /* NOTE: RDSUM only makes sense when ZTGSY2 is called by CTGSYL. */
113
114 /* RDSCAL (input/output) DOUBLE PRECISION */
115 /* On entry, scaling factor used to prevent overflow in RDSUM. */
116 /* On exit, RDSCAL is updated w.r.t. the current contributions */
117 /* in RDSUM. */
118 /* If TRANS = 'T', RDSCAL is not touched. */
119 /* NOTE: RDSCAL only makes sense when ZTGSY2 is called by */
120 /* ZTGSYL. */
121
122 /* IPIV (input) INTEGER array, dimension (N). */
123 /* The pivot indices; for 1 <= i <= N, row i of the */
124 /* matrix has been interchanged with row IPIV(i). */
125
126 /* JPIV (input) INTEGER array, dimension (N). */
127 /* The pivot indices; for 1 <= j <= N, column j of the */
128 /* matrix has been interchanged with column JPIV(j). */
129
130 /* Further Details */
131 /* =============== */
132
133 /* Based on contributions by */
134 /* Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
135 /* Umea University, S-901 87 Umea, Sweden. */
136
137 /* This routine is a further developed implementation of algorithm */
138 /* BSOLVE in [1] using complete pivoting in the LU factorization. */
139
140 /* [1] Bo Kagstrom and Lars Westin, */
141 /* Generalized Schur Methods with Condition Estimators for */
142 /* Solving the Generalized Sylvester Equation, IEEE Transactions */
143 /* on Automatic Control, Vol. 34, No. 7, July 1989, pp 745-751. */
144
145 /* [2] Peter Poromaa, */
146 /* On Efficient and Robust Estimators for the Separation */
147 /* between two Regular Matrix Pairs with Applications in */
148 /* Condition Estimation. Report UMINF-95.05, Department of */
149 /* Computing Science, Umea University, S-901 87 Umea, Sweden, */
150 /* 1995. */
151
152 /* ===================================================================== */
153
154 /* .. Parameters .. */
155 /* .. */
156 /* .. Local Scalars .. */
157 /* .. */
158 /* .. Local Arrays .. */
159 /* .. */
160 /* .. External Subroutines .. */
161 /* .. */
162 /* .. External Functions .. */
163 /* .. */
164 /* .. Intrinsic Functions .. */
165 /* .. */
166 /* .. Executable Statements .. */
167
168 /* Parameter adjustments */
169 z_dim1 = *ldz;
170 z_offset = 1 + z_dim1;
171 z__ -= z_offset;
172 --rhs;
173 --ipiv;
174 --jpiv;
175
176 /* Function Body */
177 if (*ijob != 2) {
178
179 /* Apply permutations IPIV to RHS */
180
181 i__1 = *n - 1;
182 zlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &ipiv[1], &c__1);
183
184 /* Solve for L-part choosing RHS either to +1 or -1. */
185
186 z__1.r = -1., z__1.i = -0.;
187 pmone.r = z__1.r, pmone.i = z__1.i;
188 i__1 = *n - 1;
189 for (j = 1; j <= i__1; ++j) {
190 i__2 = j;
191 z__1.r = rhs[i__2].r + 1., z__1.i = rhs[i__2].i + 0.;
192 bp.r = z__1.r, bp.i = z__1.i;
193 i__2 = j;
194 z__1.r = rhs[i__2].r - 1., z__1.i = rhs[i__2].i - 0.;
195 bm.r = z__1.r, bm.i = z__1.i;
196 splus = 1.;
197
198 /* Lockahead for L- part RHS(1:N-1) = +-1 */
199 /* SPLUS and SMIN computed more efficiently than in BSOLVE[1]. */
200
201 i__2 = *n - j;
202 zdotc_(&z__1, &i__2, &z__[j + 1 + j * z_dim1], &c__1, &z__[j + 1
203 + j * z_dim1], &c__1);
204 splus += z__1.r;
205 i__2 = *n - j;
206 zdotc_(&z__1, &i__2, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1],
207 &c__1);
208 sminu = z__1.r;
209 i__2 = j;
210 splus *= rhs[i__2].r;
211 if (splus > sminu) {
212 i__2 = j;
213 rhs[i__2].r = bp.r, rhs[i__2].i = bp.i;
214 } else if (sminu > splus) {
215 i__2 = j;
216 rhs[i__2].r = bm.r, rhs[i__2].i = bm.i;
217 } else {
218
219 /* In this case the updating sums are equal and we can */
220 /* choose RHS(J) +1 or -1. The first time this happens we */
221 /* choose -1, thereafter +1. This is a simple way to get */
222 /* good estimates of matrices like Byers well-known example */
223 /* (see [1]). (Not done in BSOLVE.) */
224
225 i__2 = j;
226 i__3 = j;
227 z__1.r = rhs[i__3].r + pmone.r, z__1.i = rhs[i__3].i +
228 pmone.i;
229 rhs[i__2].r = z__1.r, rhs[i__2].i = z__1.i;
230 pmone.r = 1., pmone.i = 0.;
231 }
232
233 /* Compute the remaining r.h.s. */
234
235 i__2 = j;
236 z__1.r = -rhs[i__2].r, z__1.i = -rhs[i__2].i;
237 temp.r = z__1.r, temp.i = z__1.i;
238 i__2 = *n - j;
239 zaxpy_(&i__2, &temp, &z__[j + 1 + j * z_dim1], &c__1, &rhs[j + 1],
240 &c__1);
241 /* L10: */
242 }
243
244 /* Solve for U- part, lockahead for RHS(N) = +-1. This is not done */
245 /* In BSOLVE and will hopefully give us a better estimate because */
246 /* any ill-conditioning of the original matrix is transfered to U */
247 /* and not to L. U(N, N) is an approximation to sigma_min(LU). */
248
249 i__1 = *n - 1;
250 zcopy_(&i__1, &rhs[1], &c__1, work, &c__1);
251 i__1 = *n - 1;
252 i__2 = *n;
253 z__1.r = rhs[i__2].r + 1., z__1.i = rhs[i__2].i + 0.;
254 work[i__1].r = z__1.r, work[i__1].i = z__1.i;
255 i__1 = *n;
256 i__2 = *n;
257 z__1.r = rhs[i__2].r - 1., z__1.i = rhs[i__2].i - 0.;
258 rhs[i__1].r = z__1.r, rhs[i__1].i = z__1.i;
259 splus = 0.;
260 sminu = 0.;
261 for (i__ = *n; i__ >= 1; --i__) {
262 z_div(&z__1, &c_b1, &z__[i__ + i__ * z_dim1]);
263 temp.r = z__1.r, temp.i = z__1.i;
264 i__1 = i__ - 1;
265 i__2 = i__ - 1;
266 z__1.r = work[i__2].r * temp.r - work[i__2].i * temp.i, z__1.i =
267 work[i__2].r * temp.i + work[i__2].i * temp.r;
268 work[i__1].r = z__1.r, work[i__1].i = z__1.i;
269 i__1 = i__;
270 i__2 = i__;
271 z__1.r = rhs[i__2].r * temp.r - rhs[i__2].i * temp.i, z__1.i =
272 rhs[i__2].r * temp.i + rhs[i__2].i * temp.r;
273 rhs[i__1].r = z__1.r, rhs[i__1].i = z__1.i;
274 i__1 = *n;
275 for (k = i__ + 1; k <= i__1; ++k) {
276 i__2 = i__ - 1;
277 i__3 = i__ - 1;
278 i__4 = k - 1;
279 i__5 = i__ + k * z_dim1;
280 z__3.r = z__[i__5].r * temp.r - z__[i__5].i * temp.i, z__3.i =
281 z__[i__5].r * temp.i + z__[i__5].i * temp.r;
282 z__2.r = work[i__4].r * z__3.r - work[i__4].i * z__3.i,
283 z__2.i = work[i__4].r * z__3.i + work[i__4].i *
284 z__3.r;
285 z__1.r = work[i__3].r - z__2.r, z__1.i = work[i__3].i -
286 z__2.i;
287 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
288 i__2 = i__;
289 i__3 = i__;
290 i__4 = k;
291 i__5 = i__ + k * z_dim1;
292 z__3.r = z__[i__5].r * temp.r - z__[i__5].i * temp.i, z__3.i =
293 z__[i__5].r * temp.i + z__[i__5].i * temp.r;
294 z__2.r = rhs[i__4].r * z__3.r - rhs[i__4].i * z__3.i, z__2.i =
295 rhs[i__4].r * z__3.i + rhs[i__4].i * z__3.r;
296 z__1.r = rhs[i__3].r - z__2.r, z__1.i = rhs[i__3].i - z__2.i;
297 rhs[i__2].r = z__1.r, rhs[i__2].i = z__1.i;
298 /* L20: */
299 }
300 splus += z_abs(&work[i__ - 1]);
301 sminu += z_abs(&rhs[i__]);
302 /* L30: */
303 }
304 if (splus > sminu) {
305 zcopy_(n, work, &c__1, &rhs[1], &c__1);
306 }
307
308 /* Apply the permutations JPIV to the computed solution (RHS) */
309
310 i__1 = *n - 1;
311 zlaswp_(&c__1, &rhs[1], ldz, &c__1, &i__1, &jpiv[1], &c_n1);
312
313 /* Compute the sum of squares */
314
315 zlassq_(n, &rhs[1], &c__1, rdscal, rdsum);
316 return 0;
317 }
318
319 /* ENTRY IJOB = 2 */
320
321 /* Compute approximate nullvector XM of Z */
322
323 zgecon_("I", n, &z__[z_offset], ldz, &c_b24, &rtemp, work, rwork, &info, (
324 ftnlen)1);
325 zcopy_(n, &work[*n], &c__1, xm, &c__1);
326
327 /* Compute RHS */
328
329 i__1 = *n - 1;
330 zlaswp_(&c__1, xm, ldz, &c__1, &i__1, &ipiv[1], &c_n1);
331 zdotc_(&z__3, n, xm, &c__1, xm, &c__1);
332 z_sqrt(&z__2, &z__3);
333 z_div(&z__1, &c_b1, &z__2);
334 temp.r = z__1.r, temp.i = z__1.i;
335 zscal_(n, &temp, xm, &c__1);
336 zcopy_(n, xm, &c__1, xp, &c__1);
337 zaxpy_(n, &c_b1, &rhs[1], &c__1, xp, &c__1);
338 z__1.r = -1., z__1.i = -0.;
339 zaxpy_(n, &z__1, xm, &c__1, &rhs[1], &c__1);
340 zgesc2_(n, &z__[z_offset], ldz, &rhs[1], &ipiv[1], &jpiv[1], &scale);
341 zgesc2_(n, &z__[z_offset], ldz, xp, &ipiv[1], &jpiv[1], &scale);
342 if (dzasum_(n, xp, &c__1) > dzasum_(n, &rhs[1], &c__1)) {
343 zcopy_(n, xp, &c__1, &rhs[1], &c__1);
344 }
345
346 /* Compute the sum of squares */
347
348 zlassq_(n, &rhs[1], &c__1, rdscal, rdsum);
349 return 0;
350
351 /* End of ZLATDF */
352
353 } /* zlatdf_ */
354
355