1 /* ./src_f77/stgsja.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 real c_b13 = 0.f;
11 static real c_b14 = 1.f;
12 static integer c__1 = 1;
13 static real c_b43 = -1.f;
14
stgsja_(char * jobu,char * jobv,char * jobq,integer * m,integer * p,integer * n,integer * k,integer * l,real * a,integer * lda,real * b,integer * ldb,real * tola,real * tolb,real * alpha,real * beta,real * u,integer * ldu,real * v,integer * ldv,real * q,integer * ldq,real * work,integer * ncycle,integer * info,ftnlen jobu_len,ftnlen jobv_len,ftnlen jobq_len)15 /* Subroutine */ int stgsja_(char *jobu, char *jobv, char *jobq, integer *m,
16 integer *p, integer *n, integer *k, integer *l, real *a, integer *lda,
17 real *b, integer *ldb, real *tola, real *tolb, real *alpha, real *
18 beta, real *u, integer *ldu, real *v, integer *ldv, real *q, integer *
19 ldq, real *work, integer *ncycle, integer *info, ftnlen jobu_len,
20 ftnlen jobv_len, ftnlen jobq_len)
21 {
22 /* System generated locals */
23 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, u_dim1,
24 u_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;
25 real r__1;
26
27 /* Local variables */
28 static integer i__, j;
29 static real a1, a2, a3, b1, b2, b3, csq, csu, csv, snq, rwk, snu, snv;
30 extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
31 integer *, real *, real *);
32 static real gamma;
33 extern logical lsame_(char *, char *, ftnlen, ftnlen);
34 extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
35 static logical initq, initu, initv, wantq, upper;
36 static real error, ssmin;
37 static logical wantu, wantv;
38 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
39 integer *), slags2_(logical *, real *, real *, real *, real *,
40 real *, real *, real *, real *, real *, real *, real *, real *);
41 static integer kcycle;
42 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), slapll_(
43 integer *, real *, integer *, real *, integer *, real *), slartg_(
44 real *, real *, real *, real *, real *), slaset_(char *, integer *
45 , integer *, real *, real *, real *, integer *, ftnlen);
46
47
48 /* -- LAPACK routine (version 3.0) -- */
49 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
50 /* Courant Institute, Argonne National Lab, and Rice University */
51 /* June 30, 1999 */
52
53 /* .. Scalar Arguments .. */
54 /* .. */
55 /* .. Array Arguments .. */
56 /* .. */
57
58 /* Purpose */
59 /* ======= */
60
61 /* STGSJA computes the generalized singular value decomposition (GSVD) */
62 /* of two real upper triangular (or trapezoidal) matrices A and B. */
63
64 /* On entry, it is assumed that matrices A and B have the following */
65 /* forms, which may be obtained by the preprocessing subroutine SGGSVP */
66 /* from a general M-by-N matrix A and P-by-N matrix B: */
67
68 /* N-K-L K L */
69 /* A = K ( 0 A12 A13 ) if M-K-L >= 0; */
70 /* L ( 0 0 A23 ) */
71 /* M-K-L ( 0 0 0 ) */
72
73 /* N-K-L K L */
74 /* A = K ( 0 A12 A13 ) if M-K-L < 0; */
75 /* M-K ( 0 0 A23 ) */
76
77 /* N-K-L K L */
78 /* B = L ( 0 0 B13 ) */
79 /* P-L ( 0 0 0 ) */
80
81 /* where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular */
82 /* upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, */
83 /* otherwise A23 is (M-K)-by-L upper trapezoidal. */
84
85 /* On exit, */
86
87 /* U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ), */
88
89 /* where U, V and Q are orthogonal matrices, Z' denotes the transpose */
90 /* of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are */
91 /* ``diagonal'' matrices, which are of the following structures: */
92
93 /* If M-K-L >= 0, */
94
95 /* K L */
96 /* D1 = K ( I 0 ) */
97 /* L ( 0 C ) */
98 /* M-K-L ( 0 0 ) */
99
100 /* K L */
101 /* D2 = L ( 0 S ) */
102 /* P-L ( 0 0 ) */
103
104 /* N-K-L K L */
105 /* ( 0 R ) = K ( 0 R11 R12 ) K */
106 /* L ( 0 0 R22 ) L */
107
108 /* where */
109
110 /* C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), */
111 /* S = diag( BETA(K+1), ... , BETA(K+L) ), */
112 /* C**2 + S**2 = I. */
113
114 /* R is stored in A(1:K+L,N-K-L+1:N) on exit. */
115
116 /* If M-K-L < 0, */
117
118 /* K M-K K+L-M */
119 /* D1 = K ( I 0 0 ) */
120 /* M-K ( 0 C 0 ) */
121
122 /* K M-K K+L-M */
123 /* D2 = M-K ( 0 S 0 ) */
124 /* K+L-M ( 0 0 I ) */
125 /* P-L ( 0 0 0 ) */
126
127 /* N-K-L K M-K K+L-M */
128 /* ( 0 R ) = K ( 0 R11 R12 R13 ) */
129 /* M-K ( 0 0 R22 R23 ) */
130 /* K+L-M ( 0 0 0 R33 ) */
131
132 /* where */
133 /* C = diag( ALPHA(K+1), ... , ALPHA(M) ), */
134 /* S = diag( BETA(K+1), ... , BETA(M) ), */
135 /* C**2 + S**2 = I. */
136
137 /* R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored */
138 /* ( 0 R22 R23 ) */
139 /* in B(M-K+1:L,N+M-K-L+1:N) on exit. */
140
141 /* The computation of the orthogonal transformation matrices U, V or Q */
142 /* is optional. These matrices may either be formed explicitly, or they */
143 /* may be postmultiplied into input matrices U1, V1, or Q1. */
144
145 /* Arguments */
146 /* ========= */
147
148 /* JOBU (input) CHARACTER*1 */
149 /* = 'U': U must contain an orthogonal matrix U1 on entry, and */
150 /* the product U1*U is returned; */
151 /* = 'I': U is initialized to the unit matrix, and the */
152 /* orthogonal matrix U is returned; */
153 /* = 'N': U is not computed. */
154
155 /* JOBV (input) CHARACTER*1 */
156 /* = 'V': V must contain an orthogonal matrix V1 on entry, and */
157 /* the product V1*V is returned; */
158 /* = 'I': V is initialized to the unit matrix, and the */
159 /* orthogonal matrix V is returned; */
160 /* = 'N': V is not computed. */
161
162 /* JOBQ (input) CHARACTER*1 */
163 /* = 'Q': Q must contain an orthogonal matrix Q1 on entry, and */
164 /* the product Q1*Q is returned; */
165 /* = 'I': Q is initialized to the unit matrix, and the */
166 /* orthogonal matrix Q is returned; */
167 /* = 'N': Q is not computed. */
168
169 /* M (input) INTEGER */
170 /* The number of rows of the matrix A. M >= 0. */
171
172 /* P (input) INTEGER */
173 /* The number of rows of the matrix B. P >= 0. */
174
175 /* N (input) INTEGER */
176 /* The number of columns of the matrices A and B. N >= 0. */
177
178 /* K (input) INTEGER */
179 /* L (input) INTEGER */
180 /* K and L specify the subblocks in the input matrices A and B: */
181 /* A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N) */
182 /* of A and B, whose GSVD is going to be computed by STGSJA. */
183 /* See Further details. */
184
185 /* A (input/output) REAL array, dimension (LDA,N) */
186 /* On entry, the M-by-N matrix A. */
187 /* On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular */
188 /* matrix R or part of R. See Purpose for details. */
189
190 /* LDA (input) INTEGER */
191 /* The leading dimension of the array A. LDA >= max(1,M). */
192
193 /* B (input/output) REAL array, dimension (LDB,N) */
194 /* On entry, the P-by-N matrix B. */
195 /* On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains */
196 /* a part of R. See Purpose for details. */
197
198 /* LDB (input) INTEGER */
199 /* The leading dimension of the array B. LDB >= max(1,P). */
200
201 /* TOLA (input) REAL */
202 /* TOLB (input) REAL */
203 /* TOLA and TOLB are the convergence criteria for the Jacobi- */
204 /* Kogbetliantz iteration procedure. Generally, they are the */
205 /* same as used in the preprocessing step, say */
206 /* TOLA = max(M,N)*norm(A)*MACHEPS, */
207 /* TOLB = max(P,N)*norm(B)*MACHEPS. */
208
209 /* ALPHA (output) REAL array, dimension (N) */
210 /* BETA (output) REAL array, dimension (N) */
211 /* On exit, ALPHA and BETA contain the generalized singular */
212 /* value pairs of A and B; */
213 /* ALPHA(1:K) = 1, */
214 /* BETA(1:K) = 0, */
215 /* and if M-K-L >= 0, */
216 /* ALPHA(K+1:K+L) = diag(C), */
217 /* BETA(K+1:K+L) = diag(S), */
218 /* or if M-K-L < 0, */
219 /* ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 */
220 /* BETA(K+1:M) = S, BETA(M+1:K+L) = 1. */
221 /* Furthermore, if K+L < N, */
222 /* ALPHA(K+L+1:N) = 0 and */
223 /* BETA(K+L+1:N) = 0. */
224
225 /* U (input/output) REAL array, dimension (LDU,M) */
226 /* On entry, if JOBU = 'U', U must contain a matrix U1 (usually */
227 /* the orthogonal matrix returned by SGGSVP). */
228 /* On exit, */
229 /* if JOBU = 'I', U contains the orthogonal matrix U; */
230 /* if JOBU = 'U', U contains the product U1*U. */
231 /* If JOBU = 'N', U is not referenced. */
232
233 /* LDU (input) INTEGER */
234 /* The leading dimension of the array U. LDU >= max(1,M) if */
235 /* JOBU = 'U'; LDU >= 1 otherwise. */
236
237 /* V (input/output) REAL array, dimension (LDV,P) */
238 /* On entry, if JOBV = 'V', V must contain a matrix V1 (usually */
239 /* the orthogonal matrix returned by SGGSVP). */
240 /* On exit, */
241 /* if JOBV = 'I', V contains the orthogonal matrix V; */
242 /* if JOBV = 'V', V contains the product V1*V. */
243 /* If JOBV = 'N', V is not referenced. */
244
245 /* LDV (input) INTEGER */
246 /* The leading dimension of the array V. LDV >= max(1,P) if */
247 /* JOBV = 'V'; LDV >= 1 otherwise. */
248
249 /* Q (input/output) REAL array, dimension (LDQ,N) */
250 /* On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually */
251 /* the orthogonal matrix returned by SGGSVP). */
252 /* On exit, */
253 /* if JOBQ = 'I', Q contains the orthogonal matrix Q; */
254 /* if JOBQ = 'Q', Q contains the product Q1*Q. */
255 /* If JOBQ = 'N', Q is not referenced. */
256
257 /* LDQ (input) INTEGER */
258 /* The leading dimension of the array Q. LDQ >= max(1,N) if */
259 /* JOBQ = 'Q'; LDQ >= 1 otherwise. */
260
261 /* WORK (workspace) REAL array, dimension (2*N) */
262
263 /* NCYCLE (output) INTEGER */
264 /* The number of cycles required for convergence. */
265
266 /* INFO (output) INTEGER */
267 /* = 0: successful exit */
268 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
269 /* = 1: the procedure does not converge after MAXIT cycles. */
270
271 /* Internal Parameters */
272 /* =================== */
273
274 /* MAXIT INTEGER */
275 /* MAXIT specifies the total loops that the iterative procedure */
276 /* may take. If after MAXIT cycles, the routine fails to */
277 /* converge, we return INFO = 1. */
278
279 /* Further Details */
280 /* =============== */
281
282 /* STGSJA essentially uses a variant of Kogbetliantz algorithm to reduce */
283 /* min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L */
284 /* matrix B13 to the form: */
285
286 /* U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1, */
287
288 /* where U1, V1 and Q1 are orthogonal matrix, and Z' is the transpose */
289 /* of Z. C1 and S1 are diagonal matrices satisfying */
290
291 /* C1**2 + S1**2 = I, */
292
293 /* and R1 is an L-by-L nonsingular upper triangular matrix. */
294
295 /* ===================================================================== */
296
297 /* .. Parameters .. */
298 /* .. */
299 /* .. Local Scalars .. */
300
301 /* .. */
302 /* .. External Functions .. */
303 /* .. */
304 /* .. External Subroutines .. */
305 /* .. */
306 /* .. Intrinsic Functions .. */
307 /* .. */
308 /* .. Executable Statements .. */
309
310 /* Decode and test the input parameters */
311
312 /* Parameter adjustments */
313 a_dim1 = *lda;
314 a_offset = 1 + a_dim1;
315 a -= a_offset;
316 b_dim1 = *ldb;
317 b_offset = 1 + b_dim1;
318 b -= b_offset;
319 --alpha;
320 --beta;
321 u_dim1 = *ldu;
322 u_offset = 1 + u_dim1;
323 u -= u_offset;
324 v_dim1 = *ldv;
325 v_offset = 1 + v_dim1;
326 v -= v_offset;
327 q_dim1 = *ldq;
328 q_offset = 1 + q_dim1;
329 q -= q_offset;
330 --work;
331
332 /* Function Body */
333 initu = lsame_(jobu, "I", (ftnlen)1, (ftnlen)1);
334 wantu = initu || lsame_(jobu, "U", (ftnlen)1, (ftnlen)1);
335
336 initv = lsame_(jobv, "I", (ftnlen)1, (ftnlen)1);
337 wantv = initv || lsame_(jobv, "V", (ftnlen)1, (ftnlen)1);
338
339 initq = lsame_(jobq, "I", (ftnlen)1, (ftnlen)1);
340 wantq = initq || lsame_(jobq, "Q", (ftnlen)1, (ftnlen)1);
341
342 *info = 0;
343 if (! (initu || wantu || lsame_(jobu, "N", (ftnlen)1, (ftnlen)1))) {
344 *info = -1;
345 } else if (! (initv || wantv || lsame_(jobv, "N", (ftnlen)1, (ftnlen)1)))
346 {
347 *info = -2;
348 } else if (! (initq || wantq || lsame_(jobq, "N", (ftnlen)1, (ftnlen)1)))
349 {
350 *info = -3;
351 } else if (*m < 0) {
352 *info = -4;
353 } else if (*p < 0) {
354 *info = -5;
355 } else if (*n < 0) {
356 *info = -6;
357 } else if (*lda < max(1,*m)) {
358 *info = -10;
359 } else if (*ldb < max(1,*p)) {
360 *info = -12;
361 } else if (*ldu < 1 || wantu && *ldu < *m) {
362 *info = -18;
363 } else if (*ldv < 1 || wantv && *ldv < *p) {
364 *info = -20;
365 } else if (*ldq < 1 || wantq && *ldq < *n) {
366 *info = -22;
367 }
368 if (*info != 0) {
369 i__1 = -(*info);
370 xerbla_("STGSJA", &i__1, (ftnlen)6);
371 return 0;
372 }
373
374 /* Initialize U, V and Q, if necessary */
375
376 if (initu) {
377 slaset_("Full", m, m, &c_b13, &c_b14, &u[u_offset], ldu, (ftnlen)4);
378 }
379 if (initv) {
380 slaset_("Full", p, p, &c_b13, &c_b14, &v[v_offset], ldv, (ftnlen)4);
381 }
382 if (initq) {
383 slaset_("Full", n, n, &c_b13, &c_b14, &q[q_offset], ldq, (ftnlen)4);
384 }
385
386 /* Loop until convergence */
387
388 upper = FALSE_;
389 for (kcycle = 1; kcycle <= 40; ++kcycle) {
390
391 upper = ! upper;
392
393 i__1 = *l - 1;
394 for (i__ = 1; i__ <= i__1; ++i__) {
395 i__2 = *l;
396 for (j = i__ + 1; j <= i__2; ++j) {
397
398 a1 = 0.f;
399 a2 = 0.f;
400 a3 = 0.f;
401 if (*k + i__ <= *m) {
402 a1 = a[*k + i__ + (*n - *l + i__) * a_dim1];
403 }
404 if (*k + j <= *m) {
405 a3 = a[*k + j + (*n - *l + j) * a_dim1];
406 }
407
408 b1 = b[i__ + (*n - *l + i__) * b_dim1];
409 b3 = b[j + (*n - *l + j) * b_dim1];
410
411 if (upper) {
412 if (*k + i__ <= *m) {
413 a2 = a[*k + i__ + (*n - *l + j) * a_dim1];
414 }
415 b2 = b[i__ + (*n - *l + j) * b_dim1];
416 } else {
417 if (*k + j <= *m) {
418 a2 = a[*k + j + (*n - *l + i__) * a_dim1];
419 }
420 b2 = b[j + (*n - *l + i__) * b_dim1];
421 }
422
423 slags2_(&upper, &a1, &a2, &a3, &b1, &b2, &b3, &csu, &snu, &
424 csv, &snv, &csq, &snq);
425
426 /* Update (K+I)-th and (K+J)-th rows of matrix A: U'*A */
427
428 if (*k + j <= *m) {
429 srot_(l, &a[*k + j + (*n - *l + 1) * a_dim1], lda, &a[*k
430 + i__ + (*n - *l + 1) * a_dim1], lda, &csu, &snu);
431 }
432
433 /* Update I-th and J-th rows of matrix B: V'*B */
434
435 srot_(l, &b[j + (*n - *l + 1) * b_dim1], ldb, &b[i__ + (*n - *
436 l + 1) * b_dim1], ldb, &csv, &snv);
437
438 /* Update (N-L+I)-th and (N-L+J)-th columns of matrices */
439 /* A and B: A*Q and B*Q */
440
441 /* Computing MIN */
442 i__4 = *k + *l;
443 i__3 = min(i__4,*m);
444 srot_(&i__3, &a[(*n - *l + j) * a_dim1 + 1], &c__1, &a[(*n - *
445 l + i__) * a_dim1 + 1], &c__1, &csq, &snq);
446
447 srot_(l, &b[(*n - *l + j) * b_dim1 + 1], &c__1, &b[(*n - *l +
448 i__) * b_dim1 + 1], &c__1, &csq, &snq);
449
450 if (upper) {
451 if (*k + i__ <= *m) {
452 a[*k + i__ + (*n - *l + j) * a_dim1] = 0.f;
453 }
454 b[i__ + (*n - *l + j) * b_dim1] = 0.f;
455 } else {
456 if (*k + j <= *m) {
457 a[*k + j + (*n - *l + i__) * a_dim1] = 0.f;
458 }
459 b[j + (*n - *l + i__) * b_dim1] = 0.f;
460 }
461
462 /* Update orthogonal matrices U, V, Q, if desired. */
463
464 if (wantu && *k + j <= *m) {
465 srot_(m, &u[(*k + j) * u_dim1 + 1], &c__1, &u[(*k + i__) *
466 u_dim1 + 1], &c__1, &csu, &snu);
467 }
468
469 if (wantv) {
470 srot_(p, &v[j * v_dim1 + 1], &c__1, &v[i__ * v_dim1 + 1],
471 &c__1, &csv, &snv);
472 }
473
474 if (wantq) {
475 srot_(n, &q[(*n - *l + j) * q_dim1 + 1], &c__1, &q[(*n - *
476 l + i__) * q_dim1 + 1], &c__1, &csq, &snq);
477 }
478
479 /* L10: */
480 }
481 /* L20: */
482 }
483
484 if (! upper) {
485
486 /* The matrices A13 and B13 were lower triangular at the start */
487 /* of the cycle, and are now upper triangular. */
488
489 /* Convergence test: test the parallelism of the corresponding */
490 /* rows of A and B. */
491
492 error = 0.f;
493 /* Computing MIN */
494 i__2 = *l, i__3 = *m - *k;
495 i__1 = min(i__2,i__3);
496 for (i__ = 1; i__ <= i__1; ++i__) {
497 i__2 = *l - i__ + 1;
498 scopy_(&i__2, &a[*k + i__ + (*n - *l + i__) * a_dim1], lda, &
499 work[1], &c__1);
500 i__2 = *l - i__ + 1;
501 scopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &work[*
502 l + 1], &c__1);
503 i__2 = *l - i__ + 1;
504 slapll_(&i__2, &work[1], &c__1, &work[*l + 1], &c__1, &ssmin);
505 error = dmax(error,ssmin);
506 /* L30: */
507 }
508
509 if (dabs(error) <= dmin(*tola,*tolb)) {
510 goto L50;
511 }
512 }
513
514 /* End of cycle loop */
515
516 /* L40: */
517 }
518
519 /* The algorithm has not converged after MAXIT cycles. */
520
521 *info = 1;
522 goto L100;
523
524 L50:
525
526 /* If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. */
527 /* Compute the generalized singular value pairs (ALPHA, BETA), and */
528 /* set the triangular matrix R to array A. */
529
530 i__1 = *k;
531 for (i__ = 1; i__ <= i__1; ++i__) {
532 alpha[i__] = 1.f;
533 beta[i__] = 0.f;
534 /* L60: */
535 }
536
537 /* Computing MIN */
538 i__2 = *l, i__3 = *m - *k;
539 i__1 = min(i__2,i__3);
540 for (i__ = 1; i__ <= i__1; ++i__) {
541
542 a1 = a[*k + i__ + (*n - *l + i__) * a_dim1];
543 b1 = b[i__ + (*n - *l + i__) * b_dim1];
544
545 if (a1 != 0.f) {
546 gamma = b1 / a1;
547
548 /* change sign if necessary */
549
550 if (gamma < 0.f) {
551 i__2 = *l - i__ + 1;
552 sscal_(&i__2, &c_b43, &b[i__ + (*n - *l + i__) * b_dim1], ldb)
553 ;
554 if (wantv) {
555 sscal_(p, &c_b43, &v[i__ * v_dim1 + 1], &c__1);
556 }
557 }
558
559 r__1 = dabs(gamma);
560 slartg_(&r__1, &c_b14, &beta[*k + i__], &alpha[*k + i__], &rwk);
561
562 if (alpha[*k + i__] >= beta[*k + i__]) {
563 i__2 = *l - i__ + 1;
564 r__1 = 1.f / alpha[*k + i__];
565 sscal_(&i__2, &r__1, &a[*k + i__ + (*n - *l + i__) * a_dim1],
566 lda);
567 } else {
568 i__2 = *l - i__ + 1;
569 r__1 = 1.f / beta[*k + i__];
570 sscal_(&i__2, &r__1, &b[i__ + (*n - *l + i__) * b_dim1], ldb);
571 i__2 = *l - i__ + 1;
572 scopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &a[*k
573 + i__ + (*n - *l + i__) * a_dim1], lda);
574 }
575
576 } else {
577
578 alpha[*k + i__] = 0.f;
579 beta[*k + i__] = 1.f;
580 i__2 = *l - i__ + 1;
581 scopy_(&i__2, &b[i__ + (*n - *l + i__) * b_dim1], ldb, &a[*k +
582 i__ + (*n - *l + i__) * a_dim1], lda);
583
584 }
585
586 /* L70: */
587 }
588
589 /* Post-assignment */
590
591 i__1 = *k + *l;
592 for (i__ = *m + 1; i__ <= i__1; ++i__) {
593 alpha[i__] = 0.f;
594 beta[i__] = 1.f;
595 /* L80: */
596 }
597
598 if (*k + *l < *n) {
599 i__1 = *n;
600 for (i__ = *k + *l + 1; i__ <= i__1; ++i__) {
601 alpha[i__] = 0.f;
602 beta[i__] = 0.f;
603 /* L90: */
604 }
605 }
606
607 L100:
608 *ncycle = kcycle;
609 return 0;
610
611 /* End of STGSJA */
612
613 } /* stgsja_ */
614
615