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