1 /* -- translated by f2c (version 20100827).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "f2c.h"
14
15 /* Table of constant values */
16
17 static doublecomplex c_b1 = {1.,0.};
18 static blasint c__1 = 1;
19
20 /** ZHETRF_ROOK_REC2 computes a partial factorization of a complex Hermitian indefinite matrix using the boun ded Bunch-Kaufman ("rook") diagonal pivoting method
21 *
22 * This routine is a minor modification of LAPACK's zlahef_rook.
23 * It serves as an unblocked kernel in the recursive algorithms.
24 * The blocked BLAS Level 3 updates were removed and moved to the
25 * recursive algorithm.
26 * */
RELAPACK_zhetrf_rook_rec2(char * uplo,blasint * n,int * nb,blasint * kb,doublecomplex * a,blasint * lda,blasint * ipiv,doublecomplex * w,blasint * ldw,blasint * info,ftnlen uplo_len)27 /* Subroutine */ void RELAPACK_zhetrf_rook_rec2(char *uplo, blasint *n,
28 int *nb, blasint *kb, doublecomplex *a, blasint *lda, blasint *
29 ipiv, doublecomplex *w, blasint *ldw, blasint *info, ftnlen uplo_len)
30 {
31 /* System generated locals */
32 blasint a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4;
33 double d__1, d__2;
34 doublecomplex z__1, z__2, z__3, z__4, z__5;
35
36 /* Builtin functions */
37 double sqrt(double), d_imag(doublecomplex *);
38 void d_cnjg(doublecomplex *, doublecomplex *), z_div(doublecomplex *,
39 doublecomplex *, doublecomplex *);
40
41 /* Local variables */
42 static blasint j, k, p;
43 static double t, r1;
44 static doublecomplex d11, d21, d22;
45 static blasint ii, jj, kk, kp, kw, jp1, jp2, kkw;
46 static logical done;
47 static blasint imax, jmax;
48 static double alpha;
49 extern logical lsame_(char *, char *, ftnlen, ftnlen);
50 static double dtemp, sfmin;
51 static blasint itemp, kstep;
52 extern /* Subroutine */ blasint zgemv_(char *, blasint *, blasint *,
53 doublecomplex *, doublecomplex *, blasint *, doublecomplex *,
54 blasint *, doublecomplex *, doublecomplex *, blasint *, ftnlen),
55 zcopy_(int *, doublecomplex *, blasint *, doublecomplex *,
56 blasint *), zswap_(int *, doublecomplex *, blasint *,
57 doublecomplex *, blasint *);
58 extern double dlamch_(char *, ftnlen);
59 static double absakk;
60 extern /* Subroutine */ blasint zdscal_(int *, double *,
61 doublecomplex *, blasint *);
62 static double colmax;
63 extern /* Subroutine */ blasint zlacgv_(int *, doublecomplex *, blasint *)
64 ;
65 extern blasint izamax_(int *, doublecomplex *, blasint *);
66 static double rowmax;
67
68 /* Parameter adjustments */
69 a_dim1 = *lda;
70 a_offset = 1 + a_dim1;
71 a -= a_offset;
72 --ipiv;
73 w_dim1 = *ldw;
74 w_offset = 1 + w_dim1;
75 w -= w_offset;
76
77 /* Function Body */
78 *info = 0;
79 alpha = (sqrt(17.) + 1.) / 8.;
80 sfmin = dlamch_("S", (ftnlen)1);
81 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
82 k = *n;
83 L10:
84 kw = *nb + k - *n;
85 if ((k <= *n - *nb + 1 && *nb < *n) || k < 1) {
86 goto L30;
87 }
88 kstep = 1;
89 p = k;
90 if (k > 1) {
91 i__1 = k - 1;
92 zcopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &w[kw * w_dim1 + 1], &
93 c__1);
94 }
95 i__1 = k + kw * w_dim1;
96 i__2 = k + k * a_dim1;
97 d__1 = a[i__2].r;
98 w[i__1].r = d__1, w[i__1].i = 0.;
99 if (k < *n) {
100 i__1 = *n - k;
101 z__1.r = -1., z__1.i = -0.;
102 zgemv_("No transpose", &k, &i__1, &z__1, &a[(k + 1) * a_dim1 + 1],
103 lda, &w[k + (kw + 1) * w_dim1], ldw, &c_b1, &w[kw *
104 w_dim1 + 1], &c__1, (ftnlen)12);
105 i__1 = k + kw * w_dim1;
106 i__2 = k + kw * w_dim1;
107 d__1 = w[i__2].r;
108 w[i__1].r = d__1, w[i__1].i = 0.;
109 }
110 i__1 = k + kw * w_dim1;
111 absakk = (d__1 = w[i__1].r, abs(d__1));
112 if (k > 1) {
113 i__1 = k - 1;
114 imax = izamax_(&i__1, &w[kw * w_dim1 + 1], &c__1);
115 i__1 = imax + kw * w_dim1;
116 colmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[imax +
117 kw * w_dim1]), abs(d__2));
118 } else {
119 colmax = 0.;
120 }
121 if (max(absakk,colmax) == 0.) {
122 if (*info == 0) {
123 *info = k;
124 }
125 kp = k;
126 i__1 = k + k * a_dim1;
127 i__2 = k + kw * w_dim1;
128 d__1 = w[i__2].r;
129 a[i__1].r = d__1, a[i__1].i = 0.;
130 if (k > 1) {
131 i__1 = k - 1;
132 zcopy_(&i__1, &w[kw * w_dim1 + 1], &c__1, &a[k * a_dim1 + 1],
133 &c__1);
134 }
135 } else {
136 if (! (absakk < alpha * colmax)) {
137 kp = k;
138 } else {
139 done = FALSE_;
140 L12:
141 if (imax > 1) {
142 i__1 = imax - 1;
143 zcopy_(&i__1, &a[imax * a_dim1 + 1], &c__1, &w[(kw - 1) *
144 w_dim1 + 1], &c__1);
145 }
146 i__1 = imax + (kw - 1) * w_dim1;
147 i__2 = imax + imax * a_dim1;
148 d__1 = a[i__2].r;
149 w[i__1].r = d__1, w[i__1].i = 0.;
150 i__1 = k - imax;
151 zcopy_(&i__1, &a[imax + (imax + 1) * a_dim1], lda, &w[imax +
152 1 + (kw - 1) * w_dim1], &c__1);
153 i__1 = k - imax;
154 zlacgv_(&i__1, &w[imax + 1 + (kw - 1) * w_dim1], &c__1);
155 if (k < *n) {
156 i__1 = *n - k;
157 z__1.r = -1., z__1.i = -0.;
158 zgemv_("No transpose", &k, &i__1, &z__1, &a[(k + 1) *
159 a_dim1 + 1], lda, &w[imax + (kw + 1) * w_dim1],
160 ldw, &c_b1, &w[(kw - 1) * w_dim1 + 1], &c__1, (
161 ftnlen)12);
162 i__1 = imax + (kw - 1) * w_dim1;
163 i__2 = imax + (kw - 1) * w_dim1;
164 d__1 = w[i__2].r;
165 w[i__1].r = d__1, w[i__1].i = 0.;
166 }
167 if (imax != k) {
168 i__1 = k - imax;
169 jmax = imax + izamax_(&i__1, &w[imax + 1 + (kw - 1) *
170 w_dim1], &c__1);
171 i__1 = jmax + (kw - 1) * w_dim1;
172 rowmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&
173 w[jmax + (kw - 1) * w_dim1]), abs(d__2));
174 } else {
175 rowmax = 0.;
176 }
177 if (imax > 1) {
178 i__1 = imax - 1;
179 itemp = izamax_(&i__1, &w[(kw - 1) * w_dim1 + 1], &c__1);
180 i__1 = itemp + (kw - 1) * w_dim1;
181 dtemp = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[
182 itemp + (kw - 1) * w_dim1]), abs(d__2));
183 if (dtemp > rowmax) {
184 rowmax = dtemp;
185 jmax = itemp;
186 }
187 }
188 i__1 = imax + (kw - 1) * w_dim1;
189 if (! ((d__1 = w[i__1].r, abs(d__1)) < alpha * rowmax)) {
190 kp = imax;
191 zcopy_(&k, &w[(kw - 1) * w_dim1 + 1], &c__1, &w[kw *
192 w_dim1 + 1], &c__1);
193 done = TRUE_;
194 } else if (p == jmax || rowmax <= colmax) {
195 kp = imax;
196 kstep = 2;
197 done = TRUE_;
198 } else {
199 p = imax;
200 colmax = rowmax;
201 imax = jmax;
202 zcopy_(&k, &w[(kw - 1) * w_dim1 + 1], &c__1, &w[kw *
203 w_dim1 + 1], &c__1);
204 }
205 if (! done) {
206 goto L12;
207 }
208 }
209 kk = k - kstep + 1;
210 kkw = *nb + kk - *n;
211 if (kstep == 2 && p != k) {
212 i__1 = p + p * a_dim1;
213 i__2 = k + k * a_dim1;
214 d__1 = a[i__2].r;
215 a[i__1].r = d__1, a[i__1].i = 0.;
216 i__1 = k - 1 - p;
217 zcopy_(&i__1, &a[p + 1 + k * a_dim1], &c__1, &a[p + (p + 1) *
218 a_dim1], lda);
219 i__1 = k - 1 - p;
220 zlacgv_(&i__1, &a[p + (p + 1) * a_dim1], lda);
221 if (p > 1) {
222 i__1 = p - 1;
223 zcopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &a[p * a_dim1 +
224 1], &c__1);
225 }
226 if (k < *n) {
227 i__1 = *n - k;
228 zswap_(&i__1, &a[k + (k + 1) * a_dim1], lda, &a[p + (k +
229 1) * a_dim1], lda);
230 }
231 i__1 = *n - kk + 1;
232 zswap_(&i__1, &w[k + kkw * w_dim1], ldw, &w[p + kkw * w_dim1],
233 ldw);
234 }
235 if (kp != kk) {
236 i__1 = kp + kp * a_dim1;
237 i__2 = kk + kk * a_dim1;
238 d__1 = a[i__2].r;
239 a[i__1].r = d__1, a[i__1].i = 0.;
240 i__1 = kk - 1 - kp;
241 zcopy_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + (kp +
242 1) * a_dim1], lda);
243 i__1 = kk - 1 - kp;
244 zlacgv_(&i__1, &a[kp + (kp + 1) * a_dim1], lda);
245 if (kp > 1) {
246 i__1 = kp - 1;
247 zcopy_(&i__1, &a[kk * a_dim1 + 1], &c__1, &a[kp * a_dim1
248 + 1], &c__1);
249 }
250 if (k < *n) {
251 i__1 = *n - k;
252 zswap_(&i__1, &a[kk + (k + 1) * a_dim1], lda, &a[kp + (k
253 + 1) * a_dim1], lda);
254 }
255 i__1 = *n - kk + 1;
256 zswap_(&i__1, &w[kk + kkw * w_dim1], ldw, &w[kp + kkw *
257 w_dim1], ldw);
258 }
259 if (kstep == 1) {
260 zcopy_(&k, &w[kw * w_dim1 + 1], &c__1, &a[k * a_dim1 + 1], &
261 c__1);
262 if (k > 1) {
263 i__1 = k + k * a_dim1;
264 t = a[i__1].r;
265 if (abs(t) >= sfmin) {
266 r1 = 1. / t;
267 i__1 = k - 1;
268 zdscal_(&i__1, &r1, &a[k * a_dim1 + 1], &c__1);
269 } else {
270 i__1 = k - 1;
271 for (ii = 1; ii <= i__1; ++ii) {
272 i__2 = ii + k * a_dim1;
273 i__3 = ii + k * a_dim1;
274 z__1.r = a[i__3].r / t, z__1.i = a[i__3].i / t;
275 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
276 /* L14: */
277 }
278 }
279 i__1 = k - 1;
280 zlacgv_(&i__1, &w[kw * w_dim1 + 1], &c__1);
281 }
282 } else {
283 if (k > 2) {
284 i__1 = k - 1 + kw * w_dim1;
285 d21.r = w[i__1].r, d21.i = w[i__1].i;
286 d_cnjg(&z__2, &d21);
287 z_div(&z__1, &w[k + kw * w_dim1], &z__2);
288 d11.r = z__1.r, d11.i = z__1.i;
289 z_div(&z__1, &w[k - 1 + (kw - 1) * w_dim1], &d21);
290 d22.r = z__1.r, d22.i = z__1.i;
291 z__1.r = d11.r * d22.r - d11.i * d22.i, z__1.i = d11.r *
292 d22.i + d11.i * d22.r;
293 t = 1. / (z__1.r - 1.);
294 i__1 = k - 2;
295 for (j = 1; j <= i__1; ++j) {
296 i__2 = j + (k - 1) * a_dim1;
297 i__3 = j + (kw - 1) * w_dim1;
298 z__4.r = d11.r * w[i__3].r - d11.i * w[i__3].i,
299 z__4.i = d11.r * w[i__3].i + d11.i * w[i__3]
300 .r;
301 i__4 = j + kw * w_dim1;
302 z__3.r = z__4.r - w[i__4].r, z__3.i = z__4.i - w[i__4]
303 .i;
304 z_div(&z__2, &z__3, &d21);
305 z__1.r = t * z__2.r, z__1.i = t * z__2.i;
306 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
307 i__2 = j + k * a_dim1;
308 i__3 = j + kw * w_dim1;
309 z__4.r = d22.r * w[i__3].r - d22.i * w[i__3].i,
310 z__4.i = d22.r * w[i__3].i + d22.i * w[i__3]
311 .r;
312 i__4 = j + (kw - 1) * w_dim1;
313 z__3.r = z__4.r - w[i__4].r, z__3.i = z__4.i - w[i__4]
314 .i;
315 d_cnjg(&z__5, &d21);
316 z_div(&z__2, &z__3, &z__5);
317 z__1.r = t * z__2.r, z__1.i = t * z__2.i;
318 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
319 /* L20: */
320 }
321 }
322 i__1 = k - 1 + (k - 1) * a_dim1;
323 i__2 = k - 1 + (kw - 1) * w_dim1;
324 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
325 i__1 = k - 1 + k * a_dim1;
326 i__2 = k - 1 + kw * w_dim1;
327 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
328 i__1 = k + k * a_dim1;
329 i__2 = k + kw * w_dim1;
330 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
331 i__1 = k - 1;
332 zlacgv_(&i__1, &w[kw * w_dim1 + 1], &c__1);
333 i__1 = k - 2;
334 zlacgv_(&i__1, &w[(kw - 1) * w_dim1 + 1], &c__1);
335 }
336 }
337 if (kstep == 1) {
338 ipiv[k] = kp;
339 } else {
340 ipiv[k] = -p;
341 ipiv[k - 1] = -kp;
342 }
343 k -= kstep;
344 goto L10;
345 L30:
346 j = k + 1;
347 L60:
348 kstep = 1;
349 jp1 = 1;
350 jj = j;
351 jp2 = ipiv[j];
352 if (jp2 < 0) {
353 jp2 = -jp2;
354 ++j;
355 jp1 = -ipiv[j];
356 kstep = 2;
357 }
358 ++j;
359 if (jp2 != jj && j <= *n) {
360 i__1 = *n - j + 1;
361 zswap_(&i__1, &a[jp2 + j * a_dim1], lda, &a[jj + j * a_dim1], lda)
362 ;
363 }
364 ++jj;
365 if (kstep == 2 && jp1 != jj && j <= *n) {
366 i__1 = *n - j + 1;
367 zswap_(&i__1, &a[jp1 + j * a_dim1], lda, &a[jj + j * a_dim1], lda)
368 ;
369 }
370 if (j < *n) {
371 goto L60;
372 }
373 *kb = *n - k;
374 } else {
375 k = 1;
376 L70:
377 if ((k >= *nb && *nb < *n) || k > *n) {
378 goto L90;
379 }
380 kstep = 1;
381 p = k;
382 i__1 = k + k * w_dim1;
383 i__2 = k + k * a_dim1;
384 d__1 = a[i__2].r;
385 w[i__1].r = d__1, w[i__1].i = 0.;
386 if (k < *n) {
387 i__1 = *n - k;
388 zcopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &w[k + 1 + k *
389 w_dim1], &c__1);
390 }
391 if (k > 1) {
392 i__1 = *n - k + 1;
393 i__2 = k - 1;
394 z__1.r = -1., z__1.i = -0.;
395 zgemv_("No transpose", &i__1, &i__2, &z__1, &a[k + a_dim1], lda, &
396 w[k + w_dim1], ldw, &c_b1, &w[k + k * w_dim1], &c__1, (
397 ftnlen)12);
398 i__1 = k + k * w_dim1;
399 i__2 = k + k * w_dim1;
400 d__1 = w[i__2].r;
401 w[i__1].r = d__1, w[i__1].i = 0.;
402 }
403 i__1 = k + k * w_dim1;
404 absakk = (d__1 = w[i__1].r, abs(d__1));
405 if (k < *n) {
406 i__1 = *n - k;
407 imax = k + izamax_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
408 i__1 = imax + k * w_dim1;
409 colmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[imax +
410 k * w_dim1]), abs(d__2));
411 } else {
412 colmax = 0.;
413 }
414 if (max(absakk,colmax) == 0.) {
415 if (*info == 0) {
416 *info = k;
417 }
418 kp = k;
419 i__1 = k + k * a_dim1;
420 i__2 = k + k * w_dim1;
421 d__1 = w[i__2].r;
422 a[i__1].r = d__1, a[i__1].i = 0.;
423 if (k < *n) {
424 i__1 = *n - k;
425 zcopy_(&i__1, &w[k + 1 + k * w_dim1], &c__1, &a[k + 1 + k *
426 a_dim1], &c__1);
427 }
428 } else {
429 if (! (absakk < alpha * colmax)) {
430 kp = k;
431 } else {
432 done = FALSE_;
433 L72:
434 i__1 = imax - k;
435 zcopy_(&i__1, &a[imax + k * a_dim1], lda, &w[k + (k + 1) *
436 w_dim1], &c__1);
437 i__1 = imax - k;
438 zlacgv_(&i__1, &w[k + (k + 1) * w_dim1], &c__1);
439 i__1 = imax + (k + 1) * w_dim1;
440 i__2 = imax + imax * a_dim1;
441 d__1 = a[i__2].r;
442 w[i__1].r = d__1, w[i__1].i = 0.;
443 if (imax < *n) {
444 i__1 = *n - imax;
445 zcopy_(&i__1, &a[imax + 1 + imax * a_dim1], &c__1, &w[
446 imax + 1 + (k + 1) * w_dim1], &c__1);
447 }
448 if (k > 1) {
449 i__1 = *n - k + 1;
450 i__2 = k - 1;
451 z__1.r = -1., z__1.i = -0.;
452 zgemv_("No transpose", &i__1, &i__2, &z__1, &a[k + a_dim1]
453 , lda, &w[imax + w_dim1], ldw, &c_b1, &w[k + (k +
454 1) * w_dim1], &c__1, (ftnlen)12);
455 i__1 = imax + (k + 1) * w_dim1;
456 i__2 = imax + (k + 1) * w_dim1;
457 d__1 = w[i__2].r;
458 w[i__1].r = d__1, w[i__1].i = 0.;
459 }
460 if (imax != k) {
461 i__1 = imax - k;
462 jmax = k - 1 + izamax_(&i__1, &w[k + (k + 1) * w_dim1], &
463 c__1);
464 i__1 = jmax + (k + 1) * w_dim1;
465 rowmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&
466 w[jmax + (k + 1) * w_dim1]), abs(d__2));
467 } else {
468 rowmax = 0.;
469 }
470 if (imax < *n) {
471 i__1 = *n - imax;
472 itemp = imax + izamax_(&i__1, &w[imax + 1 + (k + 1) *
473 w_dim1], &c__1);
474 i__1 = itemp + (k + 1) * w_dim1;
475 dtemp = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[
476 itemp + (k + 1) * w_dim1]), abs(d__2));
477 if (dtemp > rowmax) {
478 rowmax = dtemp;
479 jmax = itemp;
480 }
481 }
482 i__1 = imax + (k + 1) * w_dim1;
483 if (! ((d__1 = w[i__1].r, abs(d__1)) < alpha * rowmax)) {
484 kp = imax;
485 i__1 = *n - k + 1;
486 zcopy_(&i__1, &w[k + (k + 1) * w_dim1], &c__1, &w[k + k *
487 w_dim1], &c__1);
488 done = TRUE_;
489 } else if (p == jmax || rowmax <= colmax) {
490 kp = imax;
491 kstep = 2;
492 done = TRUE_;
493 } else {
494 p = imax;
495 colmax = rowmax;
496 imax = jmax;
497 i__1 = *n - k + 1;
498 zcopy_(&i__1, &w[k + (k + 1) * w_dim1], &c__1, &w[k + k *
499 w_dim1], &c__1);
500 }
501 if (! done) {
502 goto L72;
503 }
504 }
505 kk = k + kstep - 1;
506 if (kstep == 2 && p != k) {
507 i__1 = p + p * a_dim1;
508 i__2 = k + k * a_dim1;
509 d__1 = a[i__2].r;
510 a[i__1].r = d__1, a[i__1].i = 0.;
511 i__1 = p - k - 1;
512 zcopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &a[p + (k + 1) *
513 a_dim1], lda);
514 i__1 = p - k - 1;
515 zlacgv_(&i__1, &a[p + (k + 1) * a_dim1], lda);
516 if (p < *n) {
517 i__1 = *n - p;
518 zcopy_(&i__1, &a[p + 1 + k * a_dim1], &c__1, &a[p + 1 + p
519 * a_dim1], &c__1);
520 }
521 if (k > 1) {
522 i__1 = k - 1;
523 zswap_(&i__1, &a[k + a_dim1], lda, &a[p + a_dim1], lda);
524 }
525 zswap_(&kk, &w[k + w_dim1], ldw, &w[p + w_dim1], ldw);
526 }
527 if (kp != kk) {
528 i__1 = kp + kp * a_dim1;
529 i__2 = kk + kk * a_dim1;
530 d__1 = a[i__2].r;
531 a[i__1].r = d__1, a[i__1].i = 0.;
532 i__1 = kp - kk - 1;
533 zcopy_(&i__1, &a[kk + 1 + kk * a_dim1], &c__1, &a[kp + (kk +
534 1) * a_dim1], lda);
535 i__1 = kp - kk - 1;
536 zlacgv_(&i__1, &a[kp + (kk + 1) * a_dim1], lda);
537 if (kp < *n) {
538 i__1 = *n - kp;
539 zcopy_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + 1
540 + kp * a_dim1], &c__1);
541 }
542 if (k > 1) {
543 i__1 = k - 1;
544 zswap_(&i__1, &a[kk + a_dim1], lda, &a[kp + a_dim1], lda);
545 }
546 zswap_(&kk, &w[kk + w_dim1], ldw, &w[kp + w_dim1], ldw);
547 }
548 if (kstep == 1) {
549 i__1 = *n - k + 1;
550 zcopy_(&i__1, &w[k + k * w_dim1], &c__1, &a[k + k * a_dim1], &
551 c__1);
552 if (k < *n) {
553 i__1 = k + k * a_dim1;
554 t = a[i__1].r;
555 if (abs(t) >= sfmin) {
556 r1 = 1. / t;
557 i__1 = *n - k;
558 zdscal_(&i__1, &r1, &a[k + 1 + k * a_dim1], &c__1);
559 } else {
560 i__1 = *n;
561 for (ii = k + 1; ii <= i__1; ++ii) {
562 i__2 = ii + k * a_dim1;
563 i__3 = ii + k * a_dim1;
564 z__1.r = a[i__3].r / t, z__1.i = a[i__3].i / t;
565 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
566 /* L74: */
567 }
568 }
569 i__1 = *n - k;
570 zlacgv_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
571 }
572 } else {
573 if (k < *n - 1) {
574 i__1 = k + 1 + k * w_dim1;
575 d21.r = w[i__1].r, d21.i = w[i__1].i;
576 z_div(&z__1, &w[k + 1 + (k + 1) * w_dim1], &d21);
577 d11.r = z__1.r, d11.i = z__1.i;
578 d_cnjg(&z__2, &d21);
579 z_div(&z__1, &w[k + k * w_dim1], &z__2);
580 d22.r = z__1.r, d22.i = z__1.i;
581 z__1.r = d11.r * d22.r - d11.i * d22.i, z__1.i = d11.r *
582 d22.i + d11.i * d22.r;
583 t = 1. / (z__1.r - 1.);
584 i__1 = *n;
585 for (j = k + 2; j <= i__1; ++j) {
586 i__2 = j + k * a_dim1;
587 i__3 = j + k * w_dim1;
588 z__4.r = d11.r * w[i__3].r - d11.i * w[i__3].i,
589 z__4.i = d11.r * w[i__3].i + d11.i * w[i__3]
590 .r;
591 i__4 = j + (k + 1) * w_dim1;
592 z__3.r = z__4.r - w[i__4].r, z__3.i = z__4.i - w[i__4]
593 .i;
594 d_cnjg(&z__5, &d21);
595 z_div(&z__2, &z__3, &z__5);
596 z__1.r = t * z__2.r, z__1.i = t * z__2.i;
597 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
598 i__2 = j + (k + 1) * a_dim1;
599 i__3 = j + (k + 1) * w_dim1;
600 z__4.r = d22.r * w[i__3].r - d22.i * w[i__3].i,
601 z__4.i = d22.r * w[i__3].i + d22.i * w[i__3]
602 .r;
603 i__4 = j + k * w_dim1;
604 z__3.r = z__4.r - w[i__4].r, z__3.i = z__4.i - w[i__4]
605 .i;
606 z_div(&z__2, &z__3, &d21);
607 z__1.r = t * z__2.r, z__1.i = t * z__2.i;
608 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
609 /* L80: */
610 }
611 }
612 i__1 = k + k * a_dim1;
613 i__2 = k + k * w_dim1;
614 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
615 i__1 = k + 1 + k * a_dim1;
616 i__2 = k + 1 + k * w_dim1;
617 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
618 i__1 = k + 1 + (k + 1) * a_dim1;
619 i__2 = k + 1 + (k + 1) * w_dim1;
620 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
621 i__1 = *n - k;
622 zlacgv_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
623 i__1 = *n - k - 1;
624 zlacgv_(&i__1, &w[k + 2 + (k + 1) * w_dim1], &c__1);
625 }
626 }
627 if (kstep == 1) {
628 ipiv[k] = kp;
629 } else {
630 ipiv[k] = -p;
631 ipiv[k + 1] = -kp;
632 }
633 k += kstep;
634 goto L70;
635 L90:
636 j = k - 1;
637 L120:
638 kstep = 1;
639 jp1 = 1;
640 jj = j;
641 jp2 = ipiv[j];
642 if (jp2 < 0) {
643 jp2 = -jp2;
644 --j;
645 jp1 = -ipiv[j];
646 kstep = 2;
647 }
648 --j;
649 if (jp2 != jj && j >= 1) {
650 zswap_(&j, &a[jp2 + a_dim1], lda, &a[jj + a_dim1], lda);
651 }
652 --jj;
653 if (kstep == 2 && jp1 != jj && j >= 1) {
654 zswap_(&j, &a[jp1 + a_dim1], lda, &a[jj + a_dim1], lda);
655 }
656 if (j > 1) {
657 goto L120;
658 }
659 *kb = k - 1;
660 }
661 return;
662 }
663