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 blasint c__1 = 1;
18 static blasint c_false = FALSE_;
19 static blasint c__2 = 2;
20 static float c_b26 = 1.f;
21 static float c_b30 = 0.f;
22 static blasint c_true = TRUE_;
23
RELAPACK_strsyl_rec2(char * trana,char * tranb,blasint * isgn,int * m,blasint * n,float * a,blasint * lda,float * b,blasint * ldb,float * c__,blasint * ldc,float * scale,blasint * info,ftnlen trana_len,ftnlen tranb_len)24 void RELAPACK_strsyl_rec2(char *trana, char *tranb, blasint *isgn, int
25 *m, blasint *n, float *a, blasint *lda, float *b, blasint *ldb, float *
26 c__, blasint *ldc, float *scale, blasint *info, ftnlen trana_len,
27 ftnlen tranb_len)
28 {
29 /* System generated locals */
30 blasint a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
31 i__3, i__4;
32 float r__1, r__2;
33
34 /* Local variables */
35 static blasint j, k, l;
36 static float x[4] /* was [2][2] */;
37 static blasint k1, k2, l1, l2;
38 static float a11, db, da11, vec[4] /* was [2][2] */, dum[1], eps, sgn;
39 static blasint ierr;
40 static float smin;
41 extern float sdot_(int *, float *, blasint *, float *, blasint *);
42 static float suml, sumr;
43 extern blasint lsame_(char *, char *, ftnlen, ftnlen);
44 extern /* Subroutine */ blasint sscal_(int *, float *, float *, blasint *);
45 static blasint knext, lnext;
46 static float xnorm;
47 extern /* Subroutine */ blasint slaln2_(int *, blasint *, blasint *, float
48 *, float *, float *, blasint *, float *, float *, float *, blasint *,
49 float *, float *, float *, blasint *, float *, float *, blasint *),
50 slasy2_(int *, blasint *, blasint *, blasint *, blasint *,
51 float *, blasint *, float *, blasint *, float *, blasint *, float *,
52 float *, blasint *, float *, blasint *), slabad_(float *, float *);
53 static float scaloc;
54 extern float slamch_(char *, ftnlen), slange_(char *, blasint *,
55 blasint *, float *, blasint *, float *, ftnlen);
56 extern /* Subroutine */ blasint xerbla_(char *, blasint *, ftnlen);
57 static float bignum;
58 static blasint notrna, notrnb;
59 static float smlnum;
60
61 /* Parameter adjustments */
62 a_dim1 = *lda;
63 a_offset = 1 + a_dim1;
64 a -= a_offset;
65 b_dim1 = *ldb;
66 b_offset = 1 + b_dim1;
67 b -= b_offset;
68 c_dim1 = *ldc;
69 c_offset = 1 + c_dim1;
70 c__ -= c_offset;
71
72 /* Function Body */
73 notrna = lsame_(trana, "N", (ftnlen)1, (ftnlen)1);
74 notrnb = lsame_(tranb, "N", (ftnlen)1, (ftnlen)1);
75 *info = 0;
76 if (! notrna && ! lsame_(trana, "T", (ftnlen)1, (ftnlen)1) && ! lsame_(
77 trana, "C", (ftnlen)1, (ftnlen)1)) {
78 *info = -1;
79 } else if (! notrnb && ! lsame_(tranb, "T", (ftnlen)1, (ftnlen)1) && !
80 lsame_(tranb, "C", (ftnlen)1, (ftnlen)1)) {
81 *info = -2;
82 } else if (*isgn != 1 && *isgn != -1) {
83 *info = -3;
84 } else if (*m < 0) {
85 *info = -4;
86 } else if (*n < 0) {
87 *info = -5;
88 } else if (*lda < max(1,*m)) {
89 *info = -7;
90 } else if (*ldb < max(1,*n)) {
91 *info = -9;
92 } else if (*ldc < max(1,*m)) {
93 *info = -11;
94 }
95 if (*info != 0) {
96 i__1 = -(*info);
97 xerbla_("STRSYL", &i__1, (ftnlen)6);
98 return;
99 }
100 *scale = 1.f;
101 if (*m == 0 || *n == 0) {
102 return;
103 }
104 eps = slamch_("P", (ftnlen)1);
105 smlnum = slamch_("S", (ftnlen)1);
106 bignum = 1.f / smlnum;
107 slabad_(&smlnum, &bignum);
108 smlnum = smlnum * (float) (*m * *n) / eps;
109 bignum = 1.f / smlnum;
110 /* Computing MAX */
111 r__1 = smlnum, r__2 = eps * slange_("M", m, m, &a[a_offset], lda, dum, (
112 ftnlen)1), r__1 = max(r__1,r__2), r__2 = eps * slange_("M", n, n,
113 &b[b_offset], ldb, dum, (ftnlen)1);
114 smin = dmax(r__1,r__2);
115 sgn = (float) (*isgn);
116 if (notrna && notrnb) {
117 lnext = 1;
118 i__1 = *n;
119 for (l = 1; l <= i__1; ++l) {
120 if (l < lnext) {
121 goto L70;
122 }
123 if (l == *n) {
124 l1 = l;
125 l2 = l;
126 } else {
127 if (b[l + 1 + l * b_dim1] != 0.f) {
128 l1 = l;
129 l2 = l + 1;
130 lnext = l + 2;
131 } else {
132 l1 = l;
133 l2 = l;
134 lnext = l + 1;
135 }
136 }
137 knext = *m;
138 for (k = *m; k >= 1; --k) {
139 if (k > knext) {
140 goto L60;
141 }
142 if (k == 1) {
143 k1 = k;
144 k2 = k;
145 } else {
146 if (a[k + (k - 1) * a_dim1] != 0.f) {
147 k1 = k - 1;
148 k2 = k;
149 knext = k - 2;
150 } else {
151 k1 = k;
152 k2 = k;
153 knext = k - 1;
154 }
155 }
156 if (l1 == l2 && k1 == k2) {
157 i__2 = *m - k1;
158 /* Computing MIN */
159 i__3 = k1 + 1;
160 /* Computing MIN */
161 i__4 = k1 + 1;
162 suml = sdot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, &
163 c__[min(i__4,*m) + l1 * c_dim1], &c__1);
164 i__2 = l1 - 1;
165 sumr = sdot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 *
166 b_dim1 + 1], &c__1);
167 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
168 scaloc = 1.f;
169 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
170 da11 = dabs(a11);
171 if (da11 <= smin) {
172 a11 = smin;
173 da11 = smin;
174 *info = 1;
175 }
176 db = dabs(vec[0]);
177 if (da11 < 1.f && db > 1.f) {
178 if (db > bignum * da11) {
179 scaloc = 1.f / db;
180 }
181 }
182 x[0] = vec[0] * scaloc / a11;
183 if (scaloc != 1.f) {
184 i__2 = *n;
185 for (j = 1; j <= i__2; ++j) {
186 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
187 /* L10: */
188 }
189 *scale *= scaloc;
190 }
191 c__[k1 + l1 * c_dim1] = x[0];
192 } else if (l1 == l2 && k1 != k2) {
193 i__2 = *m - k2;
194 /* Computing MIN */
195 i__3 = k2 + 1;
196 /* Computing MIN */
197 i__4 = k2 + 1;
198 suml = sdot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, &
199 c__[min(i__4,*m) + l1 * c_dim1], &c__1);
200 i__2 = l1 - 1;
201 sumr = sdot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 *
202 b_dim1 + 1], &c__1);
203 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
204 i__2 = *m - k2;
205 /* Computing MIN */
206 i__3 = k2 + 1;
207 /* Computing MIN */
208 i__4 = k2 + 1;
209 suml = sdot_(&i__2, &a[k2 + min(i__3,*m) * a_dim1], lda, &
210 c__[min(i__4,*m) + l1 * c_dim1], &c__1);
211 i__2 = l1 - 1;
212 sumr = sdot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l1 *
213 b_dim1 + 1], &c__1);
214 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
215 r__1 = -sgn * b[l1 + l1 * b_dim1];
216 slaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1
217 * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &r__1,
218 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
219 if (ierr != 0) {
220 *info = 1;
221 }
222 if (scaloc != 1.f) {
223 i__2 = *n;
224 for (j = 1; j <= i__2; ++j) {
225 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
226 /* L20: */
227 }
228 *scale *= scaloc;
229 }
230 c__[k1 + l1 * c_dim1] = x[0];
231 c__[k2 + l1 * c_dim1] = x[1];
232 } else if (l1 != l2 && k1 == k2) {
233 i__2 = *m - k1;
234 /* Computing MIN */
235 i__3 = k1 + 1;
236 /* Computing MIN */
237 i__4 = k1 + 1;
238 suml = sdot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, &
239 c__[min(i__4,*m) + l1 * c_dim1], &c__1);
240 i__2 = l1 - 1;
241 sumr = sdot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 *
242 b_dim1 + 1], &c__1);
243 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn *
244 sumr));
245 i__2 = *m - k1;
246 /* Computing MIN */
247 i__3 = k1 + 1;
248 /* Computing MIN */
249 i__4 = k1 + 1;
250 suml = sdot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, &
251 c__[min(i__4,*m) + l2 * c_dim1], &c__1);
252 i__2 = l1 - 1;
253 sumr = sdot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l2 *
254 b_dim1 + 1], &c__1);
255 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn *
256 sumr));
257 r__1 = -sgn * a[k1 + k1 * a_dim1];
258 slaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 *
259 b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &r__1,
260 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
261 if (ierr != 0) {
262 *info = 1;
263 }
264 if (scaloc != 1.f) {
265 i__2 = *n;
266 for (j = 1; j <= i__2; ++j) {
267 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
268 /* L40: */
269 }
270 *scale *= scaloc;
271 }
272 c__[k1 + l1 * c_dim1] = x[0];
273 c__[k1 + l2 * c_dim1] = x[1];
274 } else if (l1 != l2 && k1 != k2) {
275 i__2 = *m - k2;
276 /* Computing MIN */
277 i__3 = k2 + 1;
278 /* Computing MIN */
279 i__4 = k2 + 1;
280 suml = sdot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, &
281 c__[min(i__4,*m) + l1 * c_dim1], &c__1);
282 i__2 = l1 - 1;
283 sumr = sdot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 *
284 b_dim1 + 1], &c__1);
285 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
286 i__2 = *m - k2;
287 /* Computing MIN */
288 i__3 = k2 + 1;
289 /* Computing MIN */
290 i__4 = k2 + 1;
291 suml = sdot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, &
292 c__[min(i__4,*m) + l2 * c_dim1], &c__1);
293 i__2 = l1 - 1;
294 sumr = sdot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l2 *
295 b_dim1 + 1], &c__1);
296 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
297 i__2 = *m - k2;
298 /* Computing MIN */
299 i__3 = k2 + 1;
300 /* Computing MIN */
301 i__4 = k2 + 1;
302 suml = sdot_(&i__2, &a[k2 + min(i__3,*m) * a_dim1], lda, &
303 c__[min(i__4,*m) + l1 * c_dim1], &c__1);
304 i__2 = l1 - 1;
305 sumr = sdot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l1 *
306 b_dim1 + 1], &c__1);
307 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
308 i__2 = *m - k2;
309 /* Computing MIN */
310 i__3 = k2 + 1;
311 /* Computing MIN */
312 i__4 = k2 + 1;
313 suml = sdot_(&i__2, &a[k2 + min(i__3,*m) * a_dim1], lda, &
314 c__[min(i__4,*m) + l2 * c_dim1], &c__1);
315 i__2 = l1 - 1;
316 sumr = sdot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l2 *
317 b_dim1 + 1], &c__1);
318 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
319 slasy2_(&c_false, &c_false, isgn, &c__2, &c__2, &a[k1 +
320 k1 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec,
321 &c__2, &scaloc, x, &c__2, &xnorm, &ierr);
322 if (ierr != 0) {
323 *info = 1;
324 }
325 if (scaloc != 1.f) {
326 i__2 = *n;
327 for (j = 1; j <= i__2; ++j) {
328 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
329 /* L50: */
330 }
331 *scale *= scaloc;
332 }
333 c__[k1 + l1 * c_dim1] = x[0];
334 c__[k1 + l2 * c_dim1] = x[2];
335 c__[k2 + l1 * c_dim1] = x[1];
336 c__[k2 + l2 * c_dim1] = x[3];
337 }
338 L60:
339 ;
340 }
341 L70:
342 ;
343 }
344 } else if (! notrna && notrnb) {
345 lnext = 1;
346 i__1 = *n;
347 for (l = 1; l <= i__1; ++l) {
348 if (l < lnext) {
349 goto L130;
350 }
351 if (l == *n) {
352 l1 = l;
353 l2 = l;
354 } else {
355 if (b[l + 1 + l * b_dim1] != 0.f) {
356 l1 = l;
357 l2 = l + 1;
358 lnext = l + 2;
359 } else {
360 l1 = l;
361 l2 = l;
362 lnext = l + 1;
363 }
364 }
365 knext = 1;
366 i__2 = *m;
367 for (k = 1; k <= i__2; ++k) {
368 if (k < knext) {
369 goto L120;
370 }
371 if (k == *m) {
372 k1 = k;
373 k2 = k;
374 } else {
375 if (a[k + 1 + k * a_dim1] != 0.f) {
376 k1 = k;
377 k2 = k + 1;
378 knext = k + 2;
379 } else {
380 k1 = k;
381 k2 = k;
382 knext = k + 1;
383 }
384 }
385 if (l1 == l2 && k1 == k2) {
386 i__3 = k1 - 1;
387 suml = sdot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
388 c_dim1 + 1], &c__1);
389 i__3 = l1 - 1;
390 sumr = sdot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 *
391 b_dim1 + 1], &c__1);
392 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
393 scaloc = 1.f;
394 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
395 da11 = dabs(a11);
396 if (da11 <= smin) {
397 a11 = smin;
398 da11 = smin;
399 *info = 1;
400 }
401 db = dabs(vec[0]);
402 if (da11 < 1.f && db > 1.f) {
403 if (db > bignum * da11) {
404 scaloc = 1.f / db;
405 }
406 }
407 x[0] = vec[0] * scaloc / a11;
408 if (scaloc != 1.f) {
409 i__3 = *n;
410 for (j = 1; j <= i__3; ++j) {
411 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
412 /* L80: */
413 }
414 *scale *= scaloc;
415 }
416 c__[k1 + l1 * c_dim1] = x[0];
417 } else if (l1 == l2 && k1 != k2) {
418 i__3 = k1 - 1;
419 suml = sdot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
420 c_dim1 + 1], &c__1);
421 i__3 = l1 - 1;
422 sumr = sdot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 *
423 b_dim1 + 1], &c__1);
424 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
425 i__3 = k1 - 1;
426 suml = sdot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 *
427 c_dim1 + 1], &c__1);
428 i__3 = l1 - 1;
429 sumr = sdot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l1 *
430 b_dim1 + 1], &c__1);
431 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
432 r__1 = -sgn * b[l1 + l1 * b_dim1];
433 slaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 *
434 a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &r__1,
435 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
436 if (ierr != 0) {
437 *info = 1;
438 }
439 if (scaloc != 1.f) {
440 i__3 = *n;
441 for (j = 1; j <= i__3; ++j) {
442 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
443 /* L90: */
444 }
445 *scale *= scaloc;
446 }
447 c__[k1 + l1 * c_dim1] = x[0];
448 c__[k2 + l1 * c_dim1] = x[1];
449 } else if (l1 != l2 && k1 == k2) {
450 i__3 = k1 - 1;
451 suml = sdot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
452 c_dim1 + 1], &c__1);
453 i__3 = l1 - 1;
454 sumr = sdot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 *
455 b_dim1 + 1], &c__1);
456 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn *
457 sumr));
458 i__3 = k1 - 1;
459 suml = sdot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 *
460 c_dim1 + 1], &c__1);
461 i__3 = l1 - 1;
462 sumr = sdot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l2 *
463 b_dim1 + 1], &c__1);
464 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn *
465 sumr));
466 r__1 = -sgn * a[k1 + k1 * a_dim1];
467 slaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 *
468 b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &r__1,
469 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
470 if (ierr != 0) {
471 *info = 1;
472 }
473 if (scaloc != 1.f) {
474 i__3 = *n;
475 for (j = 1; j <= i__3; ++j) {
476 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
477 /* L100: */
478 }
479 *scale *= scaloc;
480 }
481 c__[k1 + l1 * c_dim1] = x[0];
482 c__[k1 + l2 * c_dim1] = x[1];
483 } else if (l1 != l2 && k1 != k2) {
484 i__3 = k1 - 1;
485 suml = sdot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
486 c_dim1 + 1], &c__1);
487 i__3 = l1 - 1;
488 sumr = sdot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 *
489 b_dim1 + 1], &c__1);
490 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
491 i__3 = k1 - 1;
492 suml = sdot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 *
493 c_dim1 + 1], &c__1);
494 i__3 = l1 - 1;
495 sumr = sdot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l2 *
496 b_dim1 + 1], &c__1);
497 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
498 i__3 = k1 - 1;
499 suml = sdot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 *
500 c_dim1 + 1], &c__1);
501 i__3 = l1 - 1;
502 sumr = sdot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l1 *
503 b_dim1 + 1], &c__1);
504 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
505 i__3 = k1 - 1;
506 suml = sdot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l2 *
507 c_dim1 + 1], &c__1);
508 i__3 = l1 - 1;
509 sumr = sdot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l2 *
510 b_dim1 + 1], &c__1);
511 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
512 slasy2_(&c_true, &c_false, isgn, &c__2, &c__2, &a[k1 + k1
513 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &
514 c__2, &scaloc, x, &c__2, &xnorm, &ierr);
515 if (ierr != 0) {
516 *info = 1;
517 }
518 if (scaloc != 1.f) {
519 i__3 = *n;
520 for (j = 1; j <= i__3; ++j) {
521 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
522 /* L110: */
523 }
524 *scale *= scaloc;
525 }
526 c__[k1 + l1 * c_dim1] = x[0];
527 c__[k1 + l2 * c_dim1] = x[2];
528 c__[k2 + l1 * c_dim1] = x[1];
529 c__[k2 + l2 * c_dim1] = x[3];
530 }
531 L120:
532 ;
533 }
534 L130:
535 ;
536 }
537 } else if (! notrna && ! notrnb) {
538 lnext = *n;
539 for (l = *n; l >= 1; --l) {
540 if (l > lnext) {
541 goto L190;
542 }
543 if (l == 1) {
544 l1 = l;
545 l2 = l;
546 } else {
547 if (b[l + (l - 1) * b_dim1] != 0.f) {
548 l1 = l - 1;
549 l2 = l;
550 lnext = l - 2;
551 } else {
552 l1 = l;
553 l2 = l;
554 lnext = l - 1;
555 }
556 }
557 knext = 1;
558 i__1 = *m;
559 for (k = 1; k <= i__1; ++k) {
560 if (k < knext) {
561 goto L180;
562 }
563 if (k == *m) {
564 k1 = k;
565 k2 = k;
566 } else {
567 if (a[k + 1 + k * a_dim1] != 0.f) {
568 k1 = k;
569 k2 = k + 1;
570 knext = k + 2;
571 } else {
572 k1 = k;
573 k2 = k;
574 knext = k + 1;
575 }
576 }
577 if (l1 == l2 && k1 == k2) {
578 i__2 = k1 - 1;
579 suml = sdot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
580 c_dim1 + 1], &c__1);
581 i__2 = *n - l1;
582 /* Computing MIN */
583 i__3 = l1 + 1;
584 /* Computing MIN */
585 i__4 = l1 + 1;
586 sumr = sdot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc,
587 &b[l1 + min(i__4,*n) * b_dim1], ldb);
588 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
589 scaloc = 1.f;
590 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
591 da11 = dabs(a11);
592 if (da11 <= smin) {
593 a11 = smin;
594 da11 = smin;
595 *info = 1;
596 }
597 db = dabs(vec[0]);
598 if (da11 < 1.f && db > 1.f) {
599 if (db > bignum * da11) {
600 scaloc = 1.f / db;
601 }
602 }
603 x[0] = vec[0] * scaloc / a11;
604 if (scaloc != 1.f) {
605 i__2 = *n;
606 for (j = 1; j <= i__2; ++j) {
607 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
608 /* L140: */
609 }
610 *scale *= scaloc;
611 }
612 c__[k1 + l1 * c_dim1] = x[0];
613 } else if (l1 == l2 && k1 != k2) {
614 i__2 = k1 - 1;
615 suml = sdot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
616 c_dim1 + 1], &c__1);
617 i__2 = *n - l2;
618 /* Computing MIN */
619 i__3 = l2 + 1;
620 /* Computing MIN */
621 i__4 = l2 + 1;
622 sumr = sdot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc,
623 &b[l1 + min(i__4,*n) * b_dim1], ldb);
624 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
625 i__2 = k1 - 1;
626 suml = sdot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 *
627 c_dim1 + 1], &c__1);
628 i__2 = *n - l2;
629 /* Computing MIN */
630 i__3 = l2 + 1;
631 /* Computing MIN */
632 i__4 = l2 + 1;
633 sumr = sdot_(&i__2, &c__[k2 + min(i__3,*n) * c_dim1], ldc,
634 &b[l1 + min(i__4,*n) * b_dim1], ldb);
635 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
636 r__1 = -sgn * b[l1 + l1 * b_dim1];
637 slaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 *
638 a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &r__1,
639 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
640 if (ierr != 0) {
641 *info = 1;
642 }
643 if (scaloc != 1.f) {
644 i__2 = *n;
645 for (j = 1; j <= i__2; ++j) {
646 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
647 /* L150: */
648 }
649 *scale *= scaloc;
650 }
651 c__[k1 + l1 * c_dim1] = x[0];
652 c__[k2 + l1 * c_dim1] = x[1];
653 } else if (l1 != l2 && k1 == k2) {
654 i__2 = k1 - 1;
655 suml = sdot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
656 c_dim1 + 1], &c__1);
657 i__2 = *n - l2;
658 /* Computing MIN */
659 i__3 = l2 + 1;
660 /* Computing MIN */
661 i__4 = l2 + 1;
662 sumr = sdot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc,
663 &b[l1 + min(i__4,*n) * b_dim1], ldb);
664 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn *
665 sumr));
666 i__2 = k1 - 1;
667 suml = sdot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 *
668 c_dim1 + 1], &c__1);
669 i__2 = *n - l2;
670 /* Computing MIN */
671 i__3 = l2 + 1;
672 /* Computing MIN */
673 i__4 = l2 + 1;
674 sumr = sdot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc,
675 &b[l2 + min(i__4,*n) * b_dim1], ldb);
676 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn *
677 sumr));
678 r__1 = -sgn * a[k1 + k1 * a_dim1];
679 slaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1
680 * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &r__1,
681 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
682 if (ierr != 0) {
683 *info = 1;
684 }
685 if (scaloc != 1.f) {
686 i__2 = *n;
687 for (j = 1; j <= i__2; ++j) {
688 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
689 /* L160: */
690 }
691 *scale *= scaloc;
692 }
693 c__[k1 + l1 * c_dim1] = x[0];
694 c__[k1 + l2 * c_dim1] = x[1];
695 } else if (l1 != l2 && k1 != k2) {
696 i__2 = k1 - 1;
697 suml = sdot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 *
698 c_dim1 + 1], &c__1);
699 i__2 = *n - l2;
700 /* Computing MIN */
701 i__3 = l2 + 1;
702 /* Computing MIN */
703 i__4 = l2 + 1;
704 sumr = sdot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc,
705 &b[l1 + min(i__4,*n) * b_dim1], ldb);
706 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
707 i__2 = k1 - 1;
708 suml = sdot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 *
709 c_dim1 + 1], &c__1);
710 i__2 = *n - l2;
711 /* Computing MIN */
712 i__3 = l2 + 1;
713 /* Computing MIN */
714 i__4 = l2 + 1;
715 sumr = sdot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc,
716 &b[l2 + min(i__4,*n) * b_dim1], ldb);
717 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
718 i__2 = k1 - 1;
719 suml = sdot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 *
720 c_dim1 + 1], &c__1);
721 i__2 = *n - l2;
722 /* Computing MIN */
723 i__3 = l2 + 1;
724 /* Computing MIN */
725 i__4 = l2 + 1;
726 sumr = sdot_(&i__2, &c__[k2 + min(i__3,*n) * c_dim1], ldc,
727 &b[l1 + min(i__4,*n) * b_dim1], ldb);
728 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
729 i__2 = k1 - 1;
730 suml = sdot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l2 *
731 c_dim1 + 1], &c__1);
732 i__2 = *n - l2;
733 /* Computing MIN */
734 i__3 = l2 + 1;
735 /* Computing MIN */
736 i__4 = l2 + 1;
737 sumr = sdot_(&i__2, &c__[k2 + min(i__3,*n) * c_dim1], ldc,
738 &b[l2 + min(i__4,*n) * b_dim1], ldb);
739 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
740 slasy2_(&c_true, &c_true, isgn, &c__2, &c__2, &a[k1 + k1 *
741 a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &
742 c__2, &scaloc, x, &c__2, &xnorm, &ierr);
743 if (ierr != 0) {
744 *info = 1;
745 }
746 if (scaloc != 1.f) {
747 i__2 = *n;
748 for (j = 1; j <= i__2; ++j) {
749 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
750 /* L170: */
751 }
752 *scale *= scaloc;
753 }
754 c__[k1 + l1 * c_dim1] = x[0];
755 c__[k1 + l2 * c_dim1] = x[2];
756 c__[k2 + l1 * c_dim1] = x[1];
757 c__[k2 + l2 * c_dim1] = x[3];
758 }
759 L180:
760 ;
761 }
762 L190:
763 ;
764 }
765 } else if (notrna && ! notrnb) {
766 lnext = *n;
767 for (l = *n; l >= 1; --l) {
768 if (l > lnext) {
769 goto L250;
770 }
771 if (l == 1) {
772 l1 = l;
773 l2 = l;
774 } else {
775 if (b[l + (l - 1) * b_dim1] != 0.f) {
776 l1 = l - 1;
777 l2 = l;
778 lnext = l - 2;
779 } else {
780 l1 = l;
781 l2 = l;
782 lnext = l - 1;
783 }
784 }
785 knext = *m;
786 for (k = *m; k >= 1; --k) {
787 if (k > knext) {
788 goto L240;
789 }
790 if (k == 1) {
791 k1 = k;
792 k2 = k;
793 } else {
794 if (a[k + (k - 1) * a_dim1] != 0.f) {
795 k1 = k - 1;
796 k2 = k;
797 knext = k - 2;
798 } else {
799 k1 = k;
800 k2 = k;
801 knext = k - 1;
802 }
803 }
804 if (l1 == l2 && k1 == k2) {
805 i__1 = *m - k1;
806 /* Computing MIN */
807 i__2 = k1 + 1;
808 /* Computing MIN */
809 i__3 = k1 + 1;
810 suml = sdot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, &
811 c__[min(i__3,*m) + l1 * c_dim1], &c__1);
812 i__1 = *n - l1;
813 /* Computing MIN */
814 i__2 = l1 + 1;
815 /* Computing MIN */
816 i__3 = l1 + 1;
817 sumr = sdot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc,
818 &b[l1 + min(i__3,*n) * b_dim1], ldb);
819 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
820 scaloc = 1.f;
821 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
822 da11 = dabs(a11);
823 if (da11 <= smin) {
824 a11 = smin;
825 da11 = smin;
826 *info = 1;
827 }
828 db = dabs(vec[0]);
829 if (da11 < 1.f && db > 1.f) {
830 if (db > bignum * da11) {
831 scaloc = 1.f / db;
832 }
833 }
834 x[0] = vec[0] * scaloc / a11;
835 if (scaloc != 1.f) {
836 i__1 = *n;
837 for (j = 1; j <= i__1; ++j) {
838 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
839 /* L200: */
840 }
841 *scale *= scaloc;
842 }
843 c__[k1 + l1 * c_dim1] = x[0];
844 } else if (l1 == l2 && k1 != k2) {
845 i__1 = *m - k2;
846 /* Computing MIN */
847 i__2 = k2 + 1;
848 /* Computing MIN */
849 i__3 = k2 + 1;
850 suml = sdot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, &
851 c__[min(i__3,*m) + l1 * c_dim1], &c__1);
852 i__1 = *n - l2;
853 /* Computing MIN */
854 i__2 = l2 + 1;
855 /* Computing MIN */
856 i__3 = l2 + 1;
857 sumr = sdot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc,
858 &b[l1 + min(i__3,*n) * b_dim1], ldb);
859 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
860 i__1 = *m - k2;
861 /* Computing MIN */
862 i__2 = k2 + 1;
863 /* Computing MIN */
864 i__3 = k2 + 1;
865 suml = sdot_(&i__1, &a[k2 + min(i__2,*m) * a_dim1], lda, &
866 c__[min(i__3,*m) + l1 * c_dim1], &c__1);
867 i__1 = *n - l2;
868 /* Computing MIN */
869 i__2 = l2 + 1;
870 /* Computing MIN */
871 i__3 = l2 + 1;
872 sumr = sdot_(&i__1, &c__[k2 + min(i__2,*n) * c_dim1], ldc,
873 &b[l1 + min(i__3,*n) * b_dim1], ldb);
874 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
875 r__1 = -sgn * b[l1 + l1 * b_dim1];
876 slaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1
877 * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &r__1,
878 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
879 if (ierr != 0) {
880 *info = 1;
881 }
882 if (scaloc != 1.f) {
883 i__1 = *n;
884 for (j = 1; j <= i__1; ++j) {
885 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
886 /* L210: */
887 }
888 *scale *= scaloc;
889 }
890 c__[k1 + l1 * c_dim1] = x[0];
891 c__[k2 + l1 * c_dim1] = x[1];
892 } else if (l1 != l2 && k1 == k2) {
893 i__1 = *m - k1;
894 /* Computing MIN */
895 i__2 = k1 + 1;
896 /* Computing MIN */
897 i__3 = k1 + 1;
898 suml = sdot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, &
899 c__[min(i__3,*m) + l1 * c_dim1], &c__1);
900 i__1 = *n - l2;
901 /* Computing MIN */
902 i__2 = l2 + 1;
903 /* Computing MIN */
904 i__3 = l2 + 1;
905 sumr = sdot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc,
906 &b[l1 + min(i__3,*n) * b_dim1], ldb);
907 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn *
908 sumr));
909 i__1 = *m - k1;
910 /* Computing MIN */
911 i__2 = k1 + 1;
912 /* Computing MIN */
913 i__3 = k1 + 1;
914 suml = sdot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, &
915 c__[min(i__3,*m) + l2 * c_dim1], &c__1);
916 i__1 = *n - l2;
917 /* Computing MIN */
918 i__2 = l2 + 1;
919 /* Computing MIN */
920 i__3 = l2 + 1;
921 sumr = sdot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc,
922 &b[l2 + min(i__3,*n) * b_dim1], ldb);
923 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn *
924 sumr));
925 r__1 = -sgn * a[k1 + k1 * a_dim1];
926 slaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1
927 * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &r__1,
928 &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
929 if (ierr != 0) {
930 *info = 1;
931 }
932 if (scaloc != 1.f) {
933 i__1 = *n;
934 for (j = 1; j <= i__1; ++j) {
935 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
936 /* L220: */
937 }
938 *scale *= scaloc;
939 }
940 c__[k1 + l1 * c_dim1] = x[0];
941 c__[k1 + l2 * c_dim1] = x[1];
942 } else if (l1 != l2 && k1 != k2) {
943 i__1 = *m - k2;
944 /* Computing MIN */
945 i__2 = k2 + 1;
946 /* Computing MIN */
947 i__3 = k2 + 1;
948 suml = sdot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, &
949 c__[min(i__3,*m) + l1 * c_dim1], &c__1);
950 i__1 = *n - l2;
951 /* Computing MIN */
952 i__2 = l2 + 1;
953 /* Computing MIN */
954 i__3 = l2 + 1;
955 sumr = sdot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc,
956 &b[l1 + min(i__3,*n) * b_dim1], ldb);
957 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
958 i__1 = *m - k2;
959 /* Computing MIN */
960 i__2 = k2 + 1;
961 /* Computing MIN */
962 i__3 = k2 + 1;
963 suml = sdot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, &
964 c__[min(i__3,*m) + l2 * c_dim1], &c__1);
965 i__1 = *n - l2;
966 /* Computing MIN */
967 i__2 = l2 + 1;
968 /* Computing MIN */
969 i__3 = l2 + 1;
970 sumr = sdot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc,
971 &b[l2 + min(i__3,*n) * b_dim1], ldb);
972 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
973 i__1 = *m - k2;
974 /* Computing MIN */
975 i__2 = k2 + 1;
976 /* Computing MIN */
977 i__3 = k2 + 1;
978 suml = sdot_(&i__1, &a[k2 + min(i__2,*m) * a_dim1], lda, &
979 c__[min(i__3,*m) + l1 * c_dim1], &c__1);
980 i__1 = *n - l2;
981 /* Computing MIN */
982 i__2 = l2 + 1;
983 /* Computing MIN */
984 i__3 = l2 + 1;
985 sumr = sdot_(&i__1, &c__[k2 + min(i__2,*n) * c_dim1], ldc,
986 &b[l1 + min(i__3,*n) * b_dim1], ldb);
987 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
988 i__1 = *m - k2;
989 /* Computing MIN */
990 i__2 = k2 + 1;
991 /* Computing MIN */
992 i__3 = k2 + 1;
993 suml = sdot_(&i__1, &a[k2 + min(i__2,*m) * a_dim1], lda, &
994 c__[min(i__3,*m) + l2 * c_dim1], &c__1);
995 i__1 = *n - l2;
996 /* Computing MIN */
997 i__2 = l2 + 1;
998 /* Computing MIN */
999 i__3 = l2 + 1;
1000 sumr = sdot_(&i__1, &c__[k2 + min(i__2,*n) * c_dim1], ldc,
1001 &b[l2 + min(i__3,*n) * b_dim1], ldb);
1002 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
1003 slasy2_(&c_false, &c_true, isgn, &c__2, &c__2, &a[k1 + k1
1004 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &
1005 c__2, &scaloc, x, &c__2, &xnorm, &ierr);
1006 if (ierr != 0) {
1007 *info = 1;
1008 }
1009 if (scaloc != 1.f) {
1010 i__1 = *n;
1011 for (j = 1; j <= i__1; ++j) {
1012 sscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
1013 /* L230: */
1014 }
1015 *scale *= scaloc;
1016 }
1017 c__[k1 + l1 * c_dim1] = x[0];
1018 c__[k1 + l2 * c_dim1] = x[2];
1019 c__[k2 + l1 * c_dim1] = x[1];
1020 c__[k2 + l2 * c_dim1] = x[3];
1021 }
1022 L240:
1023 ;
1024 }
1025 L250:
1026 ;
1027 }
1028 }
1029 }
1030