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_REC2 computes a partial factorization of a complex Hermitian indefinite matrix using the Bunch-Kau fman diagonal pivoting method
21 *
22 * This routine is a minor modification of LAPACK's zlahef.
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_rec2(char * uplo,blasint * n,blasint * nb,blasint * kb,doublecomplex * a,blasint * lda,blasint * ipiv,doublecomplex * w,blasint * ldw,blasint * info,ftnlen uplo_len)27 /* Subroutine */ void RELAPACK_zhetrf_rec2(char *uplo, blasint *n, blasint *
28 nb, blasint *kb, doublecomplex *a, blasint *lda, blasint *ipiv,
29 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, d__3, d__4;
34 doublecomplex z__1, z__2, z__3, z__4;
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;
43 static double t, r1;
44 static doublecomplex d11, d21, d22;
45 static blasint jj, kk, jp, kp, kw, kkw, imax, jmax;
46 static double alpha;
47 extern logical lsame_(char *, char *, ftnlen, ftnlen);
48 static blasint kstep;
49 extern /* Subroutine */ blasint zgemv_(char *, blasint *, blasint *,
50 doublecomplex *, doublecomplex *, blasint *, doublecomplex *,
51 blasint *, doublecomplex *, doublecomplex *, blasint *, ftnlen),
52 zcopy_(int *, doublecomplex *, blasint *, doublecomplex *,
53 blasint *), zswap_(int *, doublecomplex *, blasint *,
54 doublecomplex *, blasint *);
55 static double absakk;
56 extern /* Subroutine */ blasint zdscal_(int *, double *,
57 doublecomplex *, blasint *);
58 static double colmax;
59 extern /* Subroutine */ blasint zlacgv_(int *, doublecomplex *, blasint *)
60 ;
61 extern blasint izamax_(int *, doublecomplex *, blasint *);
62 static double rowmax;
63
64 /* Parameter adjustments */
65 a_dim1 = *lda;
66 a_offset = 1 + a_dim1;
67 a -= a_offset;
68 --ipiv;
69 w_dim1 = *ldw;
70 w_offset = 1 + w_dim1;
71 w -= w_offset;
72
73 /* Function Body */
74 *info = 0;
75 alpha = (sqrt(17.) + 1.) / 8.;
76 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
77 k = *n;
78 L10:
79 kw = *nb + k - *n;
80 if ((k <= *n - *nb + 1 && *nb < *n) || k < 1) {
81 goto L30;
82 }
83 kstep = 1;
84 i__1 = k - 1;
85 zcopy_(&i__1, &a[k * a_dim1 + 1], &c__1, &w[kw * w_dim1 + 1], &c__1);
86 i__1 = k + kw * w_dim1;
87 i__2 = k + k * a_dim1;
88 d__1 = a[i__2].r;
89 w[i__1].r = d__1, w[i__1].i = 0.;
90 if (k < *n) {
91 i__1 = *n - k;
92 z__1.r = -1., z__1.i = -0.;
93 zgemv_("No transpose", &k, &i__1, &z__1, &a[(k + 1) * a_dim1 + 1],
94 lda, &w[k + (kw + 1) * w_dim1], ldw, &c_b1, &w[kw *
95 w_dim1 + 1], &c__1, (ftnlen)12);
96 i__1 = k + kw * w_dim1;
97 i__2 = k + kw * w_dim1;
98 d__1 = w[i__2].r;
99 w[i__1].r = d__1, w[i__1].i = 0.;
100 }
101 i__1 = k + kw * w_dim1;
102 absakk = (d__1 = w[i__1].r, abs(d__1));
103 if (k > 1) {
104 i__1 = k - 1;
105 imax = izamax_(&i__1, &w[kw * w_dim1 + 1], &c__1);
106 i__1 = imax + kw * w_dim1;
107 colmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[imax +
108 kw * w_dim1]), abs(d__2));
109 } else {
110 colmax = 0.;
111 }
112 if (max(absakk,colmax) == 0.) {
113 if (*info == 0) {
114 *info = k;
115 }
116 kp = k;
117 i__1 = k + k * a_dim1;
118 i__2 = k + k * a_dim1;
119 d__1 = a[i__2].r;
120 a[i__1].r = d__1, a[i__1].i = 0.;
121 } else {
122 if (absakk >= alpha * colmax) {
123 kp = k;
124 } else {
125 i__1 = imax - 1;
126 zcopy_(&i__1, &a[imax * a_dim1 + 1], &c__1, &w[(kw - 1) *
127 w_dim1 + 1], &c__1);
128 i__1 = imax + (kw - 1) * w_dim1;
129 i__2 = imax + imax * a_dim1;
130 d__1 = a[i__2].r;
131 w[i__1].r = d__1, w[i__1].i = 0.;
132 i__1 = k - imax;
133 zcopy_(&i__1, &a[imax + (imax + 1) * a_dim1], lda, &w[imax +
134 1 + (kw - 1) * w_dim1], &c__1);
135 i__1 = k - imax;
136 zlacgv_(&i__1, &w[imax + 1 + (kw - 1) * w_dim1], &c__1);
137 if (k < *n) {
138 i__1 = *n - k;
139 z__1.r = -1., z__1.i = -0.;
140 zgemv_("No transpose", &k, &i__1, &z__1, &a[(k + 1) *
141 a_dim1 + 1], lda, &w[imax + (kw + 1) * w_dim1],
142 ldw, &c_b1, &w[(kw - 1) * w_dim1 + 1], &c__1, (
143 ftnlen)12);
144 i__1 = imax + (kw - 1) * w_dim1;
145 i__2 = imax + (kw - 1) * w_dim1;
146 d__1 = w[i__2].r;
147 w[i__1].r = d__1, w[i__1].i = 0.;
148 }
149 i__1 = k - imax;
150 jmax = imax + izamax_(&i__1, &w[imax + 1 + (kw - 1) * w_dim1],
151 &c__1);
152 i__1 = jmax + (kw - 1) * w_dim1;
153 rowmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[
154 jmax + (kw - 1) * w_dim1]), abs(d__2));
155 if (imax > 1) {
156 i__1 = imax - 1;
157 jmax = izamax_(&i__1, &w[(kw - 1) * w_dim1 + 1], &c__1);
158 /* Computing MAX */
159 i__1 = jmax + (kw - 1) * w_dim1;
160 d__3 = rowmax, d__4 = (d__1 = w[i__1].r, abs(d__1)) + (
161 d__2 = d_imag(&w[jmax + (kw - 1) * w_dim1]), abs(
162 d__2));
163 rowmax = max(d__3,d__4);
164 }
165 if (absakk >= alpha * colmax * (colmax / rowmax)) {
166 kp = k;
167 } else /* if(complicated condition) */ {
168 i__1 = imax + (kw - 1) * w_dim1;
169 if ((d__1 = w[i__1].r, abs(d__1)) >= alpha * rowmax) {
170 kp = imax;
171 zcopy_(&k, &w[(kw - 1) * w_dim1 + 1], &c__1, &w[kw *
172 w_dim1 + 1], &c__1);
173 } else {
174 kp = imax;
175 kstep = 2;
176 }
177 }
178 }
179 kk = k - kstep + 1;
180 kkw = *nb + kk - *n;
181 if (kp != kk) {
182 i__1 = kp + kp * a_dim1;
183 i__2 = kk + kk * a_dim1;
184 d__1 = a[i__2].r;
185 a[i__1].r = d__1, a[i__1].i = 0.;
186 i__1 = kk - 1 - kp;
187 zcopy_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + (kp +
188 1) * a_dim1], lda);
189 i__1 = kk - 1 - kp;
190 zlacgv_(&i__1, &a[kp + (kp + 1) * a_dim1], lda);
191 if (kp > 1) {
192 i__1 = kp - 1;
193 zcopy_(&i__1, &a[kk * a_dim1 + 1], &c__1, &a[kp * a_dim1
194 + 1], &c__1);
195 }
196 if (k < *n) {
197 i__1 = *n - k;
198 zswap_(&i__1, &a[kk + (k + 1) * a_dim1], lda, &a[kp + (k
199 + 1) * a_dim1], lda);
200 }
201 i__1 = *n - kk + 1;
202 zswap_(&i__1, &w[kk + kkw * w_dim1], ldw, &w[kp + kkw *
203 w_dim1], ldw);
204 }
205 if (kstep == 1) {
206 zcopy_(&k, &w[kw * w_dim1 + 1], &c__1, &a[k * a_dim1 + 1], &
207 c__1);
208 if (k > 1) {
209 i__1 = k + k * a_dim1;
210 r1 = 1. / a[i__1].r;
211 i__1 = k - 1;
212 zdscal_(&i__1, &r1, &a[k * a_dim1 + 1], &c__1);
213 i__1 = k - 1;
214 zlacgv_(&i__1, &w[kw * w_dim1 + 1], &c__1);
215 }
216 } else {
217 if (k > 2) {
218 i__1 = k - 1 + kw * w_dim1;
219 d21.r = w[i__1].r, d21.i = w[i__1].i;
220 d_cnjg(&z__2, &d21);
221 z_div(&z__1, &w[k + kw * w_dim1], &z__2);
222 d11.r = z__1.r, d11.i = z__1.i;
223 z_div(&z__1, &w[k - 1 + (kw - 1) * w_dim1], &d21);
224 d22.r = z__1.r, d22.i = z__1.i;
225 z__1.r = d11.r * d22.r - d11.i * d22.i, z__1.i = d11.r *
226 d22.i + d11.i * d22.r;
227 t = 1. / (z__1.r - 1.);
228 z__2.r = t, z__2.i = 0.;
229 z_div(&z__1, &z__2, &d21);
230 d21.r = z__1.r, d21.i = z__1.i;
231 i__1 = k - 2;
232 for (j = 1; j <= i__1; ++j) {
233 i__2 = j + (k - 1) * a_dim1;
234 i__3 = j + (kw - 1) * w_dim1;
235 z__3.r = d11.r * w[i__3].r - d11.i * w[i__3].i,
236 z__3.i = d11.r * w[i__3].i + d11.i * w[i__3]
237 .r;
238 i__4 = j + kw * w_dim1;
239 z__2.r = z__3.r - w[i__4].r, z__2.i = z__3.i - w[i__4]
240 .i;
241 z__1.r = d21.r * z__2.r - d21.i * z__2.i, z__1.i =
242 d21.r * z__2.i + d21.i * z__2.r;
243 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
244 i__2 = j + k * a_dim1;
245 d_cnjg(&z__2, &d21);
246 i__3 = j + kw * w_dim1;
247 z__4.r = d22.r * w[i__3].r - d22.i * w[i__3].i,
248 z__4.i = d22.r * w[i__3].i + d22.i * w[i__3]
249 .r;
250 i__4 = j + (kw - 1) * w_dim1;
251 z__3.r = z__4.r - w[i__4].r, z__3.i = z__4.i - w[i__4]
252 .i;
253 z__1.r = z__2.r * z__3.r - z__2.i * z__3.i, z__1.i =
254 z__2.r * z__3.i + z__2.i * z__3.r;
255 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
256 /* L20: */
257 }
258 }
259 i__1 = k - 1 + (k - 1) * a_dim1;
260 i__2 = k - 1 + (kw - 1) * w_dim1;
261 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
262 i__1 = k - 1 + k * a_dim1;
263 i__2 = k - 1 + kw * w_dim1;
264 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
265 i__1 = k + k * a_dim1;
266 i__2 = k + kw * w_dim1;
267 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
268 i__1 = k - 1;
269 zlacgv_(&i__1, &w[kw * w_dim1 + 1], &c__1);
270 i__1 = k - 2;
271 zlacgv_(&i__1, &w[(kw - 1) * w_dim1 + 1], &c__1);
272 }
273 }
274 if (kstep == 1) {
275 ipiv[k] = kp;
276 } else {
277 ipiv[k] = -kp;
278 ipiv[k - 1] = -kp;
279 }
280 k -= kstep;
281 goto L10;
282 L30:
283 j = k + 1;
284 L60:
285 jj = j;
286 jp = ipiv[j];
287 if (jp < 0) {
288 jp = -jp;
289 ++j;
290 }
291 ++j;
292 if (jp != jj && j <= *n) {
293 i__1 = *n - j + 1;
294 zswap_(&i__1, &a[jp + j * a_dim1], lda, &a[jj + j * a_dim1], lda);
295 }
296 if (j < *n) {
297 goto L60;
298 }
299 *kb = *n - k;
300 } else {
301 k = 1;
302 L70:
303 if ((k >= *nb && *nb < *n) || k > *n) {
304 goto L90;
305 }
306 kstep = 1;
307 i__1 = k + k * w_dim1;
308 i__2 = k + k * a_dim1;
309 d__1 = a[i__2].r;
310 w[i__1].r = d__1, w[i__1].i = 0.;
311 if (k < *n) {
312 i__1 = *n - k;
313 zcopy_(&i__1, &a[k + 1 + k * a_dim1], &c__1, &w[k + 1 + k *
314 w_dim1], &c__1);
315 }
316 i__1 = *n - k + 1;
317 i__2 = k - 1;
318 z__1.r = -1., z__1.i = -0.;
319 zgemv_("No transpose", &i__1, &i__2, &z__1, &a[k + a_dim1], lda, &w[k
320 + w_dim1], ldw, &c_b1, &w[k + k * w_dim1], &c__1, (ftnlen)12);
321 i__1 = k + k * w_dim1;
322 i__2 = k + k * w_dim1;
323 d__1 = w[i__2].r;
324 w[i__1].r = d__1, w[i__1].i = 0.;
325 i__1 = k + k * w_dim1;
326 absakk = (d__1 = w[i__1].r, abs(d__1));
327 if (k < *n) {
328 i__1 = *n - k;
329 imax = k + izamax_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
330 i__1 = imax + k * w_dim1;
331 colmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[imax +
332 k * w_dim1]), abs(d__2));
333 } else {
334 colmax = 0.;
335 }
336 if (max(absakk,colmax) == 0.) {
337 if (*info == 0) {
338 *info = k;
339 }
340 kp = k;
341 i__1 = k + k * a_dim1;
342 i__2 = k + k * a_dim1;
343 d__1 = a[i__2].r;
344 a[i__1].r = d__1, a[i__1].i = 0.;
345 } else {
346 if (absakk >= alpha * colmax) {
347 kp = k;
348 } else {
349 i__1 = imax - k;
350 zcopy_(&i__1, &a[imax + k * a_dim1], lda, &w[k + (k + 1) *
351 w_dim1], &c__1);
352 i__1 = imax - k;
353 zlacgv_(&i__1, &w[k + (k + 1) * w_dim1], &c__1);
354 i__1 = imax + (k + 1) * w_dim1;
355 i__2 = imax + imax * a_dim1;
356 d__1 = a[i__2].r;
357 w[i__1].r = d__1, w[i__1].i = 0.;
358 if (imax < *n) {
359 i__1 = *n - imax;
360 zcopy_(&i__1, &a[imax + 1 + imax * a_dim1], &c__1, &w[
361 imax + 1 + (k + 1) * w_dim1], &c__1);
362 }
363 i__1 = *n - k + 1;
364 i__2 = k - 1;
365 z__1.r = -1., z__1.i = -0.;
366 zgemv_("No transpose", &i__1, &i__2, &z__1, &a[k + a_dim1],
367 lda, &w[imax + w_dim1], ldw, &c_b1, &w[k + (k + 1) *
368 w_dim1], &c__1, (ftnlen)12);
369 i__1 = imax + (k + 1) * w_dim1;
370 i__2 = imax + (k + 1) * w_dim1;
371 d__1 = w[i__2].r;
372 w[i__1].r = d__1, w[i__1].i = 0.;
373 i__1 = imax - k;
374 jmax = k - 1 + izamax_(&i__1, &w[k + (k + 1) * w_dim1], &c__1)
375 ;
376 i__1 = jmax + (k + 1) * w_dim1;
377 rowmax = (d__1 = w[i__1].r, abs(d__1)) + (d__2 = d_imag(&w[
378 jmax + (k + 1) * w_dim1]), abs(d__2));
379 if (imax < *n) {
380 i__1 = *n - imax;
381 jmax = imax + izamax_(&i__1, &w[imax + 1 + (k + 1) *
382 w_dim1], &c__1);
383 /* Computing MAX */
384 i__1 = jmax + (k + 1) * w_dim1;
385 d__3 = rowmax, d__4 = (d__1 = w[i__1].r, abs(d__1)) + (
386 d__2 = d_imag(&w[jmax + (k + 1) * w_dim1]), abs(
387 d__2));
388 rowmax = max(d__3,d__4);
389 }
390 if (absakk >= alpha * colmax * (colmax / rowmax)) {
391 kp = k;
392 } else /* if(complicated condition) */ {
393 i__1 = imax + (k + 1) * w_dim1;
394 if ((d__1 = w[i__1].r, abs(d__1)) >= alpha * rowmax) {
395 kp = imax;
396 i__1 = *n - k + 1;
397 zcopy_(&i__1, &w[k + (k + 1) * w_dim1], &c__1, &w[k +
398 k * w_dim1], &c__1);
399 } else {
400 kp = imax;
401 kstep = 2;
402 }
403 }
404 }
405 kk = k + kstep - 1;
406 if (kp != kk) {
407 i__1 = kp + kp * a_dim1;
408 i__2 = kk + kk * a_dim1;
409 d__1 = a[i__2].r;
410 a[i__1].r = d__1, a[i__1].i = 0.;
411 i__1 = kp - kk - 1;
412 zcopy_(&i__1, &a[kk + 1 + kk * a_dim1], &c__1, &a[kp + (kk +
413 1) * a_dim1], lda);
414 i__1 = kp - kk - 1;
415 zlacgv_(&i__1, &a[kp + (kk + 1) * a_dim1], lda);
416 if (kp < *n) {
417 i__1 = *n - kp;
418 zcopy_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + 1
419 + kp * a_dim1], &c__1);
420 }
421 if (k > 1) {
422 i__1 = k - 1;
423 zswap_(&i__1, &a[kk + a_dim1], lda, &a[kp + a_dim1], lda);
424 }
425 zswap_(&kk, &w[kk + w_dim1], ldw, &w[kp + w_dim1], ldw);
426 }
427 if (kstep == 1) {
428 i__1 = *n - k + 1;
429 zcopy_(&i__1, &w[k + k * w_dim1], &c__1, &a[k + k * a_dim1], &
430 c__1);
431 if (k < *n) {
432 i__1 = k + k * a_dim1;
433 r1 = 1. / a[i__1].r;
434 i__1 = *n - k;
435 zdscal_(&i__1, &r1, &a[k + 1 + k * a_dim1], &c__1);
436 i__1 = *n - k;
437 zlacgv_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
438 }
439 } else {
440 if (k < *n - 1) {
441 i__1 = k + 1 + k * w_dim1;
442 d21.r = w[i__1].r, d21.i = w[i__1].i;
443 z_div(&z__1, &w[k + 1 + (k + 1) * w_dim1], &d21);
444 d11.r = z__1.r, d11.i = z__1.i;
445 d_cnjg(&z__2, &d21);
446 z_div(&z__1, &w[k + k * w_dim1], &z__2);
447 d22.r = z__1.r, d22.i = z__1.i;
448 z__1.r = d11.r * d22.r - d11.i * d22.i, z__1.i = d11.r *
449 d22.i + d11.i * d22.r;
450 t = 1. / (z__1.r - 1.);
451 z__2.r = t, z__2.i = 0.;
452 z_div(&z__1, &z__2, &d21);
453 d21.r = z__1.r, d21.i = z__1.i;
454 i__1 = *n;
455 for (j = k + 2; j <= i__1; ++j) {
456 i__2 = j + k * a_dim1;
457 d_cnjg(&z__2, &d21);
458 i__3 = j + k * w_dim1;
459 z__4.r = d11.r * w[i__3].r - d11.i * w[i__3].i,
460 z__4.i = d11.r * w[i__3].i + d11.i * w[i__3]
461 .r;
462 i__4 = j + (k + 1) * w_dim1;
463 z__3.r = z__4.r - w[i__4].r, z__3.i = z__4.i - w[i__4]
464 .i;
465 z__1.r = z__2.r * z__3.r - z__2.i * z__3.i, z__1.i =
466 z__2.r * z__3.i + z__2.i * z__3.r;
467 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
468 i__2 = j + (k + 1) * a_dim1;
469 i__3 = j + (k + 1) * w_dim1;
470 z__3.r = d22.r * w[i__3].r - d22.i * w[i__3].i,
471 z__3.i = d22.r * w[i__3].i + d22.i * w[i__3]
472 .r;
473 i__4 = j + k * w_dim1;
474 z__2.r = z__3.r - w[i__4].r, z__2.i = z__3.i - w[i__4]
475 .i;
476 z__1.r = d21.r * z__2.r - d21.i * z__2.i, z__1.i =
477 d21.r * z__2.i + d21.i * z__2.r;
478 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
479 /* L80: */
480 }
481 }
482 i__1 = k + k * a_dim1;
483 i__2 = k + k * w_dim1;
484 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
485 i__1 = k + 1 + k * a_dim1;
486 i__2 = k + 1 + k * w_dim1;
487 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
488 i__1 = k + 1 + (k + 1) * a_dim1;
489 i__2 = k + 1 + (k + 1) * w_dim1;
490 a[i__1].r = w[i__2].r, a[i__1].i = w[i__2].i;
491 i__1 = *n - k;
492 zlacgv_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
493 i__1 = *n - k - 1;
494 zlacgv_(&i__1, &w[k + 2 + (k + 1) * w_dim1], &c__1);
495 }
496 }
497 if (kstep == 1) {
498 ipiv[k] = kp;
499 } else {
500 ipiv[k] = -kp;
501 ipiv[k + 1] = -kp;
502 }
503 k += kstep;
504 goto L70;
505 L90:
506 j = k - 1;
507 L120:
508 jj = j;
509 jp = ipiv[j];
510 if (jp < 0) {
511 jp = -jp;
512 --j;
513 }
514 --j;
515 if (jp != jj && j >= 1) {
516 zswap_(&j, &a[jp + a_dim1], lda, &a[jj + a_dim1], lda);
517 }
518 if (j > 1) {
519 goto L120;
520 }
521 *kb = k - 1;
522 }
523 return;
524 }
525