1 /* ../netlib/dtrsyl.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static integer c__1 = 1;
5 static logical c_false = FALSE_;
6 static integer c__2 = 2;
7 static doublereal c_b26 = 1.;
8 static doublereal c_b30 = 0.;
9 static logical c_true = TRUE_;
10 /* > \brief \b DTRSYL */
11 /* =========== DOCUMENTATION =========== */
12 /* Online html documentation available at */
13 /* http://www.netlib.org/lapack/explore-html/ */
14 /* > \htmlonly */
15 /* > Download DTRSYL + dependencies */
16 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl. f"> */
17 /* > [TGZ]</a> */
18 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl. f"> */
19 /* > [ZIP]</a> */
20 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl. f"> */
21 /* > [TXT]</a> */
22 /* > \endhtmlonly */
23 /* Definition: */
24 /* =========== */
25 /* SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, */
26 /* LDC, SCALE, INFO ) */
27 /* .. Scalar Arguments .. */
28 /* CHARACTER TRANA, TRANB */
29 /* INTEGER INFO, ISGN, LDA, LDB, LDC, M, N */
30 /* DOUBLE PRECISION SCALE */
31 /* .. */
32 /* .. Array Arguments .. */
33 /* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) */
34 /* .. */
35 /* > \par Purpose: */
36 /* ============= */
37 /* > */
38 /* > \verbatim */
39 /* > */
40 /* > DTRSYL solves the real Sylvester matrix equation: */
41 /* > */
42 /* > op(A)*X + X*op(B) = scale*C or */
43 /* > op(A)*X - X*op(B) = scale*C, */
44 /* > */
45 /* > where op(A) = A or A**T, and A and B are both upper quasi- */
46 /* > triangular. A is M-by-M and B is N-by-N;
47 the right hand side C and */
48 /* > the solution X are M-by-N;
49 and scale is an output scale factor, set */
50 /* > <= 1 to avoid overflow in X. */
51 /* > */
52 /* > A and B must be in Schur canonical form (as returned by DHSEQR), that */
53 /* > is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
54 */
55 /* > each 2-by-2 diagonal block has its diagonal elements equal and its */
56 /* > off-diagonal elements of opposite sign. */
57 /* > \endverbatim */
58 /* Arguments: */
59 /* ========== */
60 /* > \param[in] TRANA */
61 /* > \verbatim */
62 /* > TRANA is CHARACTER*1 */
63 /* > Specifies the option op(A): */
64 /* > = 'N': op(A) = A (No transpose) */
65 /* > = 'T': op(A) = A**T (Transpose) */
66 /* > = 'C': op(A) = A**H (Conjugate transpose = Transpose) */
67 /* > \endverbatim */
68 /* > */
69 /* > \param[in] TRANB */
70 /* > \verbatim */
71 /* > TRANB is CHARACTER*1 */
72 /* > Specifies the option op(B): */
73 /* > = 'N': op(B) = B (No transpose) */
74 /* > = 'T': op(B) = B**T (Transpose) */
75 /* > = 'C': op(B) = B**H (Conjugate transpose = Transpose) */
76 /* > \endverbatim */
77 /* > */
78 /* > \param[in] ISGN */
79 /* > \verbatim */
80 /* > ISGN is INTEGER */
81 /* > Specifies the sign in the equation: */
82 /* > = +1: solve op(A)*X + X*op(B) = scale*C */
83 /* > = -1: solve op(A)*X - X*op(B) = scale*C */
84 /* > \endverbatim */
85 /* > */
86 /* > \param[in] M */
87 /* > \verbatim */
88 /* > M is INTEGER */
89 /* > The order of the matrix A, and the number of rows in the */
90 /* > matrices X and C. M >= 0. */
91 /* > \endverbatim */
92 /* > */
93 /* > \param[in] N */
94 /* > \verbatim */
95 /* > N is INTEGER */
96 /* > The order of the matrix B, and the number of columns in the */
97 /* > matrices X and C. N >= 0. */
98 /* > \endverbatim */
99 /* > */
100 /* > \param[in] A */
101 /* > \verbatim */
102 /* > A is DOUBLE PRECISION array, dimension (LDA,M) */
103 /* > The upper quasi-triangular matrix A, in Schur canonical form. */
104 /* > \endverbatim */
105 /* > */
106 /* > \param[in] LDA */
107 /* > \verbatim */
108 /* > LDA is INTEGER */
109 /* > The leading dimension of the array A. LDA >= max(1,M). */
110 /* > \endverbatim */
111 /* > */
112 /* > \param[in] B */
113 /* > \verbatim */
114 /* > B is DOUBLE PRECISION array, dimension (LDB,N) */
115 /* > The upper quasi-triangular matrix B, in Schur canonical form. */
116 /* > \endverbatim */
117 /* > */
118 /* > \param[in] LDB */
119 /* > \verbatim */
120 /* > LDB is INTEGER */
121 /* > The leading dimension of the array B. LDB >= max(1,N). */
122 /* > \endverbatim */
123 /* > */
124 /* > \param[in,out] C */
125 /* > \verbatim */
126 /* > C is DOUBLE PRECISION array, dimension (LDC,N) */
127 /* > On entry, the M-by-N right hand side matrix C. */
128 /* > On exit, C is overwritten by the solution matrix X. */
129 /* > \endverbatim */
130 /* > */
131 /* > \param[in] LDC */
132 /* > \verbatim */
133 /* > LDC is INTEGER */
134 /* > The leading dimension of the array C. LDC >= max(1,M) */
135 /* > \endverbatim */
136 /* > */
137 /* > \param[out] SCALE */
138 /* > \verbatim */
139 /* > SCALE is DOUBLE PRECISION */
140 /* > The scale factor, scale, set <= 1 to avoid overflow in X. */
141 /* > \endverbatim */
142 /* > */
143 /* > \param[out] INFO */
144 /* > \verbatim */
145 /* > INFO is INTEGER */
146 /* > = 0: successful exit */
147 /* > < 0: if INFO = -i, the i-th argument had an illegal value */
148 /* > = 1: A and B have common or very close eigenvalues;
149 perturbed */
150 /* > values were used to solve the equation (but the matrices */
151 /* > A and B are unchanged). */
152 /* > \endverbatim */
153 /* Authors: */
154 /* ======== */
155 /* > \author Univ. of Tennessee */
156 /* > \author Univ. of California Berkeley */
157 /* > \author Univ. of Colorado Denver */
158 /* > \author NAG Ltd. */
159 /* > \date November 2011 */
160 /* > \ingroup doubleSYcomputational */
161 /* ===================================================================== */
162 /* Subroutine */
dtrsyl_(char * trana,char * tranb,integer * isgn,integer * m,integer * n,doublereal * a,integer * lda,doublereal * b,integer * ldb,doublereal * c__,integer * ldc,doublereal * scale,integer * info)163 int dtrsyl_(char *trana, char *tranb, integer *isgn, integer *m, integer *n, doublereal *a, integer *lda, doublereal *b, integer * ldb, doublereal *c__, integer *ldc, doublereal *scale, integer *info)
164 {
165 /* System generated locals */
166 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, i__3, i__4;
167 doublereal d__1, d__2;
168 /* Local variables */
169 integer j, k, l;
170 doublereal x[4] /* was [2][2] */
171 ;
172 integer k1, k2, l1, l2;
173 doublereal a11, db, da11, vec[4] /* was [2][2] */
174 , dum[1], eps, sgn;
175 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, integer *);
176 integer ierr;
177 doublereal smin, suml, sumr;
178 extern /* Subroutine */
179 int dscal_(integer *, doublereal *, doublereal *, integer *);
180 extern logical lsame_(char *, char *);
181 integer knext, lnext;
182 doublereal xnorm;
183 extern /* Subroutine */
184 int dlaln2_(logical *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal * , doublereal *, integer *, doublereal *, doublereal *, integer *), dlasy2_(logical *, logical *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *), dlabad_(doublereal *, doublereal *);
185 extern doublereal dlamch_(char *), dlange_(char *, integer *, integer *, doublereal *, integer *, doublereal *);
186 doublereal scaloc;
187 extern /* Subroutine */
188 int xerbla_(char *, integer *);
189 doublereal bignum;
190 logical notrna, notrnb;
191 doublereal smlnum;
192 /* -- LAPACK computational routine (version 3.4.0) -- */
193 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
194 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
195 /* November 2011 */
196 /* .. Scalar Arguments .. */
197 /* .. */
198 /* .. Array Arguments .. */
199 /* .. */
200 /* ===================================================================== */
201 /* .. Parameters .. */
202 /* .. */
203 /* .. Local Scalars .. */
204 /* .. */
205 /* .. Local Arrays .. */
206 /* .. */
207 /* .. External Functions .. */
208 /* .. */
209 /* .. External Subroutines .. */
210 /* .. */
211 /* .. Intrinsic Functions .. */
212 /* .. */
213 /* .. Executable Statements .. */
214 /* Decode and Test input parameters */
215 /* Parameter adjustments */
216 a_dim1 = *lda;
217 a_offset = 1 + a_dim1;
218 a -= a_offset;
219 b_dim1 = *ldb;
220 b_offset = 1 + b_dim1;
221 b -= b_offset;
222 c_dim1 = *ldc;
223 c_offset = 1 + c_dim1;
224 c__ -= c_offset;
225 /* Function Body */
226 notrna = lsame_(trana, "N");
227 notrnb = lsame_(tranb, "N");
228 *info = 0;
229 if (! notrna && ! lsame_(trana, "T") && ! lsame_( trana, "C"))
230 {
231 *info = -1;
232 }
233 else if (! notrnb && ! lsame_(tranb, "T") && ! lsame_(tranb, "C"))
234 {
235 *info = -2;
236 }
237 else if (*isgn != 1 && *isgn != -1)
238 {
239 *info = -3;
240 }
241 else if (*m < 0)
242 {
243 *info = -4;
244 }
245 else if (*n < 0)
246 {
247 *info = -5;
248 }
249 else if (*lda < max(1,*m))
250 {
251 *info = -7;
252 }
253 else if (*ldb < max(1,*n))
254 {
255 *info = -9;
256 }
257 else if (*ldc < max(1,*m))
258 {
259 *info = -11;
260 }
261 if (*info != 0)
262 {
263 i__1 = -(*info);
264 xerbla_("DTRSYL", &i__1);
265 return 0;
266 }
267 /* Quick return if possible */
268 *scale = 1.;
269 if (*m == 0 || *n == 0)
270 {
271 return 0;
272 }
273 /* Set constants to control overflow */
274 eps = dlamch_("P");
275 smlnum = dlamch_("S");
276 bignum = 1. / smlnum;
277 dlabad_(&smlnum, &bignum);
278 smlnum = smlnum * (doublereal) (*m * *n) / eps;
279 bignum = 1. / smlnum;
280 /* Computing MAX */
281 d__1 = smlnum, d__2 = eps * dlange_("M", m, m, &a[a_offset], lda, dum);
282 d__1 = max(d__1,d__2);
283 d__2 = eps * dlange_("M", n, n, &b[b_offset], ldb, dum); // ; expr subst
284 smin = max(d__1,d__2);
285 sgn = (doublereal) (*isgn);
286 if (notrna && notrnb)
287 {
288 /* Solve A*X + ISGN*X*B = scale*C. */
289 /* The (K,L)th block of X is determined starting from */
290 /* bottom-left corner column by column by */
291 /* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) */
292 /* Where */
293 /* M L-1 */
294 /* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. */
295 /* I=K+1 J=1 */
296 /* Start column loop (index = L) */
297 /* L1 (L2) : column index of the first (first) row of X(K,L). */
298 lnext = 1;
299 i__1 = *n;
300 for (l = 1;
301 l <= i__1;
302 ++l)
303 {
304 if (l < lnext)
305 {
306 goto L60;
307 }
308 if (l == *n)
309 {
310 l1 = l;
311 l2 = l;
312 }
313 else
314 {
315 if (b[l + 1 + l * b_dim1] != 0.)
316 {
317 l1 = l;
318 l2 = l + 1;
319 lnext = l + 2;
320 }
321 else
322 {
323 l1 = l;
324 l2 = l;
325 lnext = l + 1;
326 }
327 }
328 /* Start row loop (index = K) */
329 /* K1 (K2): row index of the first (last) row of X(K,L). */
330 knext = *m;
331 for (k = *m;
332 k >= 1;
333 --k)
334 {
335 if (k > knext)
336 {
337 goto L50;
338 }
339 if (k == 1)
340 {
341 k1 = k;
342 k2 = k;
343 }
344 else
345 {
346 if (a[k + (k - 1) * a_dim1] != 0.)
347 {
348 k1 = k - 1;
349 k2 = k;
350 knext = k - 2;
351 }
352 else
353 {
354 k1 = k;
355 k2 = k;
356 knext = k - 1;
357 }
358 }
359 if (l1 == l2 && k1 == k2)
360 {
361 i__2 = *m - k1;
362 /* Computing MIN */
363 i__3 = k1 + 1;
364 /* Computing MIN */
365 i__4 = k1 + 1;
366 suml = ddot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l1 * c_dim1], &c__1);
367 i__2 = l1 - 1;
368 sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
369 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
370 scaloc = 1.;
371 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
372 da11 = f2c_abs(a11);
373 if (da11 <= smin)
374 {
375 a11 = smin;
376 da11 = smin;
377 *info = 1;
378 }
379 db = f2c_abs(vec[0]);
380 if (da11 < 1. && db > 1.)
381 {
382 if (db > bignum * da11)
383 {
384 scaloc = 1. / db;
385 }
386 }
387 x[0] = vec[0] * scaloc / a11;
388 if (scaloc != 1.)
389 {
390 i__2 = *n;
391 for (j = 1;
392 j <= i__2;
393 ++j)
394 {
395 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
396 /* L10: */
397 }
398 *scale *= scaloc;
399 }
400 c__[k1 + l1 * c_dim1] = x[0];
401 }
402 else if (l1 == l2 && k1 != k2)
403 {
404 i__2 = *m - k2;
405 /* Computing MIN */
406 i__3 = k2 + 1;
407 /* Computing MIN */
408 i__4 = k2 + 1;
409 suml = ddot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l1 * c_dim1], &c__1);
410 i__2 = l1 - 1;
411 sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
412 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
413 i__2 = *m - k2;
414 /* Computing MIN */
415 i__3 = k2 + 1;
416 /* Computing MIN */
417 i__4 = k2 + 1;
418 suml = ddot_(&i__2, &a[k2 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l1 * c_dim1], &c__1);
419 i__2 = l1 - 1;
420 sumr = ddot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
421 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
422 d__1 = -sgn * b[l1 + l1 * b_dim1];
423 dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
424 if (ierr != 0)
425 {
426 *info = 1;
427 }
428 if (scaloc != 1.)
429 {
430 i__2 = *n;
431 for (j = 1;
432 j <= i__2;
433 ++j)
434 {
435 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
436 /* L20: */
437 }
438 *scale *= scaloc;
439 }
440 c__[k1 + l1 * c_dim1] = x[0];
441 c__[k2 + l1 * c_dim1] = x[1];
442 }
443 else if (l1 != l2 && k1 == k2)
444 {
445 i__2 = *m - k1;
446 /* Computing MIN */
447 i__3 = k1 + 1;
448 /* Computing MIN */
449 i__4 = k1 + 1;
450 suml = ddot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l1 * c_dim1], &c__1);
451 i__2 = l1 - 1;
452 sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
453 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * sumr));
454 i__2 = *m - k1;
455 /* Computing MIN */
456 i__3 = k1 + 1;
457 /* Computing MIN */
458 i__4 = k1 + 1;
459 suml = ddot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l2 * c_dim1], &c__1);
460 i__2 = l1 - 1;
461 sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l2 * b_dim1 + 1], &c__1);
462 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * sumr));
463 d__1 = -sgn * a[k1 + k1 * a_dim1];
464 dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
465 if (ierr != 0)
466 {
467 *info = 1;
468 }
469 if (scaloc != 1.)
470 {
471 i__2 = *n;
472 for (j = 1;
473 j <= i__2;
474 ++j)
475 {
476 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
477 /* L30: */
478 }
479 *scale *= scaloc;
480 }
481 c__[k1 + l1 * c_dim1] = x[0];
482 c__[k1 + l2 * c_dim1] = x[1];
483 }
484 else if (l1 != l2 && k1 != k2)
485 {
486 i__2 = *m - k2;
487 /* Computing MIN */
488 i__3 = k2 + 1;
489 /* Computing MIN */
490 i__4 = k2 + 1;
491 suml = ddot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l1 * c_dim1], &c__1);
492 i__2 = l1 - 1;
493 sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
494 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
495 i__2 = *m - k2;
496 /* Computing MIN */
497 i__3 = k2 + 1;
498 /* Computing MIN */
499 i__4 = k2 + 1;
500 suml = ddot_(&i__2, &a[k1 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l2 * c_dim1], &c__1);
501 i__2 = l1 - 1;
502 sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l2 * b_dim1 + 1], &c__1);
503 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
504 i__2 = *m - k2;
505 /* Computing MIN */
506 i__3 = k2 + 1;
507 /* Computing MIN */
508 i__4 = k2 + 1;
509 suml = ddot_(&i__2, &a[k2 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l1 * c_dim1], &c__1);
510 i__2 = l1 - 1;
511 sumr = ddot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
512 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
513 i__2 = *m - k2;
514 /* Computing MIN */
515 i__3 = k2 + 1;
516 /* Computing MIN */
517 i__4 = k2 + 1;
518 suml = ddot_(&i__2, &a[k2 + min(i__3,*m) * a_dim1], lda, & c__[min(i__4,*m) + l2 * c_dim1], &c__1);
519 i__2 = l1 - 1;
520 sumr = ddot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l2 * b_dim1 + 1], &c__1);
521 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
522 dlasy2_(&c_false, &c_false, isgn, &c__2, &c__2, &a[k1 + k1 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &c__2, &scaloc, x, &c__2, &xnorm, &ierr);
523 if (ierr != 0)
524 {
525 *info = 1;
526 }
527 if (scaloc != 1.)
528 {
529 i__2 = *n;
530 for (j = 1;
531 j <= i__2;
532 ++j)
533 {
534 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
535 /* L40: */
536 }
537 *scale *= scaloc;
538 }
539 c__[k1 + l1 * c_dim1] = x[0];
540 c__[k1 + l2 * c_dim1] = x[2];
541 c__[k2 + l1 * c_dim1] = x[1];
542 c__[k2 + l2 * c_dim1] = x[3];
543 }
544 L50:
545 ;
546 }
547 L60:
548 ;
549 }
550 }
551 else if (! notrna && notrnb)
552 {
553 /* Solve A**T *X + ISGN*X*B = scale*C. */
554 /* The (K,L)th block of X is determined starting from */
555 /* upper-left corner column by column by */
556 /* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) */
557 /* Where */
558 /* K-1 T L-1 */
559 /* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] */
560 /* I=1 J=1 */
561 /* Start column loop (index = L) */
562 /* L1 (L2): column index of the first (last) row of X(K,L) */
563 lnext = 1;
564 i__1 = *n;
565 for (l = 1;
566 l <= i__1;
567 ++l)
568 {
569 if (l < lnext)
570 {
571 goto L120;
572 }
573 if (l == *n)
574 {
575 l1 = l;
576 l2 = l;
577 }
578 else
579 {
580 if (b[l + 1 + l * b_dim1] != 0.)
581 {
582 l1 = l;
583 l2 = l + 1;
584 lnext = l + 2;
585 }
586 else
587 {
588 l1 = l;
589 l2 = l;
590 lnext = l + 1;
591 }
592 }
593 /* Start row loop (index = K) */
594 /* K1 (K2): row index of the first (last) row of X(K,L) */
595 knext = 1;
596 i__2 = *m;
597 for (k = 1;
598 k <= i__2;
599 ++k)
600 {
601 if (k < knext)
602 {
603 goto L110;
604 }
605 if (k == *m)
606 {
607 k1 = k;
608 k2 = k;
609 }
610 else
611 {
612 if (a[k + 1 + k * a_dim1] != 0.)
613 {
614 k1 = k;
615 k2 = k + 1;
616 knext = k + 2;
617 }
618 else
619 {
620 k1 = k;
621 k2 = k;
622 knext = k + 1;
623 }
624 }
625 if (l1 == l2 && k1 == k2)
626 {
627 i__3 = k1 - 1;
628 suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
629 i__3 = l1 - 1;
630 sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
631 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
632 scaloc = 1.;
633 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
634 da11 = f2c_abs(a11);
635 if (da11 <= smin)
636 {
637 a11 = smin;
638 da11 = smin;
639 *info = 1;
640 }
641 db = f2c_abs(vec[0]);
642 if (da11 < 1. && db > 1.)
643 {
644 if (db > bignum * da11)
645 {
646 scaloc = 1. / db;
647 }
648 }
649 x[0] = vec[0] * scaloc / a11;
650 if (scaloc != 1.)
651 {
652 i__3 = *n;
653 for (j = 1;
654 j <= i__3;
655 ++j)
656 {
657 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
658 /* L70: */
659 }
660 *scale *= scaloc;
661 }
662 c__[k1 + l1 * c_dim1] = x[0];
663 }
664 else if (l1 == l2 && k1 != k2)
665 {
666 i__3 = k1 - 1;
667 suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
668 i__3 = l1 - 1;
669 sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
670 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
671 i__3 = k1 - 1;
672 suml = ddot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
673 i__3 = l1 - 1;
674 sumr = ddot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
675 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
676 d__1 = -sgn * b[l1 + l1 * b_dim1];
677 dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
678 if (ierr != 0)
679 {
680 *info = 1;
681 }
682 if (scaloc != 1.)
683 {
684 i__3 = *n;
685 for (j = 1;
686 j <= i__3;
687 ++j)
688 {
689 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
690 /* L80: */
691 }
692 *scale *= scaloc;
693 }
694 c__[k1 + l1 * c_dim1] = x[0];
695 c__[k2 + l1 * c_dim1] = x[1];
696 }
697 else if (l1 != l2 && k1 == k2)
698 {
699 i__3 = k1 - 1;
700 suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
701 i__3 = l1 - 1;
702 sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
703 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * sumr));
704 i__3 = k1 - 1;
705 suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * c_dim1 + 1], &c__1);
706 i__3 = l1 - 1;
707 sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l2 * b_dim1 + 1], &c__1);
708 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * sumr));
709 d__1 = -sgn * a[k1 + k1 * a_dim1];
710 dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
711 if (ierr != 0)
712 {
713 *info = 1;
714 }
715 if (scaloc != 1.)
716 {
717 i__3 = *n;
718 for (j = 1;
719 j <= i__3;
720 ++j)
721 {
722 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
723 /* L90: */
724 }
725 *scale *= scaloc;
726 }
727 c__[k1 + l1 * c_dim1] = x[0];
728 c__[k1 + l2 * c_dim1] = x[1];
729 }
730 else if (l1 != l2 && k1 != k2)
731 {
732 i__3 = k1 - 1;
733 suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
734 i__3 = l1 - 1;
735 sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
736 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
737 i__3 = k1 - 1;
738 suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * c_dim1 + 1], &c__1);
739 i__3 = l1 - 1;
740 sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l2 * b_dim1 + 1], &c__1);
741 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
742 i__3 = k1 - 1;
743 suml = ddot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
744 i__3 = l1 - 1;
745 sumr = ddot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l1 * b_dim1 + 1], &c__1);
746 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
747 i__3 = k1 - 1;
748 suml = ddot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l2 * c_dim1 + 1], &c__1);
749 i__3 = l1 - 1;
750 sumr = ddot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l2 * b_dim1 + 1], &c__1);
751 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
752 dlasy2_(&c_true, &c_false, isgn, &c__2, &c__2, &a[k1 + k1 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, & c__2, &scaloc, x, &c__2, &xnorm, &ierr);
753 if (ierr != 0)
754 {
755 *info = 1;
756 }
757 if (scaloc != 1.)
758 {
759 i__3 = *n;
760 for (j = 1;
761 j <= i__3;
762 ++j)
763 {
764 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
765 /* L100: */
766 }
767 *scale *= scaloc;
768 }
769 c__[k1 + l1 * c_dim1] = x[0];
770 c__[k1 + l2 * c_dim1] = x[2];
771 c__[k2 + l1 * c_dim1] = x[1];
772 c__[k2 + l2 * c_dim1] = x[3];
773 }
774 L110:
775 ;
776 }
777 L120:
778 ;
779 }
780 }
781 else if (! notrna && ! notrnb)
782 {
783 /* Solve A**T*X + ISGN*X*B**T = scale*C. */
784 /* The (K,L)th block of X is determined starting from */
785 /* top-right corner column by column by */
786 /* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) */
787 /* Where */
788 /* K-1 N */
789 /* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. */
790 /* I=1 J=L+1 */
791 /* Start column loop (index = L) */
792 /* L1 (L2): column index of the first (last) row of X(K,L) */
793 lnext = *n;
794 for (l = *n;
795 l >= 1;
796 --l)
797 {
798 if (l > lnext)
799 {
800 goto L180;
801 }
802 if (l == 1)
803 {
804 l1 = l;
805 l2 = l;
806 }
807 else
808 {
809 if (b[l + (l - 1) * b_dim1] != 0.)
810 {
811 l1 = l - 1;
812 l2 = l;
813 lnext = l - 2;
814 }
815 else
816 {
817 l1 = l;
818 l2 = l;
819 lnext = l - 1;
820 }
821 }
822 /* Start row loop (index = K) */
823 /* K1 (K2): row index of the first (last) row of X(K,L) */
824 knext = 1;
825 i__1 = *m;
826 for (k = 1;
827 k <= i__1;
828 ++k)
829 {
830 if (k < knext)
831 {
832 goto L170;
833 }
834 if (k == *m)
835 {
836 k1 = k;
837 k2 = k;
838 }
839 else
840 {
841 if (a[k + 1 + k * a_dim1] != 0.)
842 {
843 k1 = k;
844 k2 = k + 1;
845 knext = k + 2;
846 }
847 else
848 {
849 k1 = k;
850 k2 = k;
851 knext = k + 1;
852 }
853 }
854 if (l1 == l2 && k1 == k2)
855 {
856 i__2 = k1 - 1;
857 suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
858 i__2 = *n - l1;
859 /* Computing MIN */
860 i__3 = l1 + 1;
861 /* Computing MIN */
862 i__4 = l1 + 1;
863 sumr = ddot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc, &b[l1 + min(i__4,*n) * b_dim1], ldb);
864 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
865 scaloc = 1.;
866 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
867 da11 = f2c_abs(a11);
868 if (da11 <= smin)
869 {
870 a11 = smin;
871 da11 = smin;
872 *info = 1;
873 }
874 db = f2c_abs(vec[0]);
875 if (da11 < 1. && db > 1.)
876 {
877 if (db > bignum * da11)
878 {
879 scaloc = 1. / db;
880 }
881 }
882 x[0] = vec[0] * scaloc / a11;
883 if (scaloc != 1.)
884 {
885 i__2 = *n;
886 for (j = 1;
887 j <= i__2;
888 ++j)
889 {
890 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
891 /* L130: */
892 }
893 *scale *= scaloc;
894 }
895 c__[k1 + l1 * c_dim1] = x[0];
896 }
897 else if (l1 == l2 && k1 != k2)
898 {
899 i__2 = k1 - 1;
900 suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
901 i__2 = *n - l2;
902 /* Computing MIN */
903 i__3 = l2 + 1;
904 /* Computing MIN */
905 i__4 = l2 + 1;
906 sumr = ddot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc, &b[l1 + min(i__4,*n) * b_dim1], ldb);
907 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
908 i__2 = k1 - 1;
909 suml = ddot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
910 i__2 = *n - l2;
911 /* Computing MIN */
912 i__3 = l2 + 1;
913 /* Computing MIN */
914 i__4 = l2 + 1;
915 sumr = ddot_(&i__2, &c__[k2 + min(i__3,*n) * c_dim1], ldc, &b[l1 + min(i__4,*n) * b_dim1], ldb);
916 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
917 d__1 = -sgn * b[l1 + l1 * b_dim1];
918 dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
919 if (ierr != 0)
920 {
921 *info = 1;
922 }
923 if (scaloc != 1.)
924 {
925 i__2 = *n;
926 for (j = 1;
927 j <= i__2;
928 ++j)
929 {
930 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
931 /* L140: */
932 }
933 *scale *= scaloc;
934 }
935 c__[k1 + l1 * c_dim1] = x[0];
936 c__[k2 + l1 * c_dim1] = x[1];
937 }
938 else if (l1 != l2 && k1 == k2)
939 {
940 i__2 = k1 - 1;
941 suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
942 i__2 = *n - l2;
943 /* Computing MIN */
944 i__3 = l2 + 1;
945 /* Computing MIN */
946 i__4 = l2 + 1;
947 sumr = ddot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc, &b[l1 + min(i__4,*n) * b_dim1], ldb);
948 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * sumr));
949 i__2 = k1 - 1;
950 suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * c_dim1 + 1], &c__1);
951 i__2 = *n - l2;
952 /* Computing MIN */
953 i__3 = l2 + 1;
954 /* Computing MIN */
955 i__4 = l2 + 1;
956 sumr = ddot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc, &b[l2 + min(i__4,*n) * b_dim1], ldb);
957 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * sumr));
958 d__1 = -sgn * a[k1 + k1 * a_dim1];
959 dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
960 if (ierr != 0)
961 {
962 *info = 1;
963 }
964 if (scaloc != 1.)
965 {
966 i__2 = *n;
967 for (j = 1;
968 j <= i__2;
969 ++j)
970 {
971 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
972 /* L150: */
973 }
974 *scale *= scaloc;
975 }
976 c__[k1 + l1 * c_dim1] = x[0];
977 c__[k1 + l2 * c_dim1] = x[1];
978 }
979 else if (l1 != l2 && k1 != k2)
980 {
981 i__2 = k1 - 1;
982 suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
983 i__2 = *n - l2;
984 /* Computing MIN */
985 i__3 = l2 + 1;
986 /* Computing MIN */
987 i__4 = l2 + 1;
988 sumr = ddot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc, &b[l1 + min(i__4,*n) * b_dim1], ldb);
989 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
990 i__2 = k1 - 1;
991 suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * c_dim1 + 1], &c__1);
992 i__2 = *n - l2;
993 /* Computing MIN */
994 i__3 = l2 + 1;
995 /* Computing MIN */
996 i__4 = l2 + 1;
997 sumr = ddot_(&i__2, &c__[k1 + min(i__3,*n) * c_dim1], ldc, &b[l2 + min(i__4,*n) * b_dim1], ldb);
998 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
999 i__2 = k1 - 1;
1000 suml = ddot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * c_dim1 + 1], &c__1);
1001 i__2 = *n - l2;
1002 /* Computing MIN */
1003 i__3 = l2 + 1;
1004 /* Computing MIN */
1005 i__4 = l2 + 1;
1006 sumr = ddot_(&i__2, &c__[k2 + min(i__3,*n) * c_dim1], ldc, &b[l1 + min(i__4,*n) * b_dim1], ldb);
1007 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
1008 i__2 = k1 - 1;
1009 suml = ddot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l2 * c_dim1 + 1], &c__1);
1010 i__2 = *n - l2;
1011 /* Computing MIN */
1012 i__3 = l2 + 1;
1013 /* Computing MIN */
1014 i__4 = l2 + 1;
1015 sumr = ddot_(&i__2, &c__[k2 + min(i__3,*n) * c_dim1], ldc, &b[l2 + min(i__4,*n) * b_dim1], ldb);
1016 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
1017 dlasy2_(&c_true, &c_true, isgn, &c__2, &c__2, &a[k1 + k1 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, & c__2, &scaloc, x, &c__2, &xnorm, &ierr);
1018 if (ierr != 0)
1019 {
1020 *info = 1;
1021 }
1022 if (scaloc != 1.)
1023 {
1024 i__2 = *n;
1025 for (j = 1;
1026 j <= i__2;
1027 ++j)
1028 {
1029 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
1030 /* L160: */
1031 }
1032 *scale *= scaloc;
1033 }
1034 c__[k1 + l1 * c_dim1] = x[0];
1035 c__[k1 + l2 * c_dim1] = x[2];
1036 c__[k2 + l1 * c_dim1] = x[1];
1037 c__[k2 + l2 * c_dim1] = x[3];
1038 }
1039 L170:
1040 ;
1041 }
1042 L180:
1043 ;
1044 }
1045 }
1046 else if (notrna && ! notrnb)
1047 {
1048 /* Solve A*X + ISGN*X*B**T = scale*C. */
1049 /* The (K,L)th block of X is determined starting from */
1050 /* bottom-right corner column by column by */
1051 /* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L) */
1052 /* Where */
1053 /* M N */
1054 /* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T]. */
1055 /* I=K+1 J=L+1 */
1056 /* Start column loop (index = L) */
1057 /* L1 (L2): column index of the first (last) row of X(K,L) */
1058 lnext = *n;
1059 for (l = *n;
1060 l >= 1;
1061 --l)
1062 {
1063 if (l > lnext)
1064 {
1065 goto L240;
1066 }
1067 if (l == 1)
1068 {
1069 l1 = l;
1070 l2 = l;
1071 }
1072 else
1073 {
1074 if (b[l + (l - 1) * b_dim1] != 0.)
1075 {
1076 l1 = l - 1;
1077 l2 = l;
1078 lnext = l - 2;
1079 }
1080 else
1081 {
1082 l1 = l;
1083 l2 = l;
1084 lnext = l - 1;
1085 }
1086 }
1087 /* Start row loop (index = K) */
1088 /* K1 (K2): row index of the first (last) row of X(K,L) */
1089 knext = *m;
1090 for (k = *m;
1091 k >= 1;
1092 --k)
1093 {
1094 if (k > knext)
1095 {
1096 goto L230;
1097 }
1098 if (k == 1)
1099 {
1100 k1 = k;
1101 k2 = k;
1102 }
1103 else
1104 {
1105 if (a[k + (k - 1) * a_dim1] != 0.)
1106 {
1107 k1 = k - 1;
1108 k2 = k;
1109 knext = k - 2;
1110 }
1111 else
1112 {
1113 k1 = k;
1114 k2 = k;
1115 knext = k - 1;
1116 }
1117 }
1118 if (l1 == l2 && k1 == k2)
1119 {
1120 i__1 = *m - k1;
1121 /* Computing MIN */
1122 i__2 = k1 + 1;
1123 /* Computing MIN */
1124 i__3 = k1 + 1;
1125 suml = ddot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l1 * c_dim1], &c__1);
1126 i__1 = *n - l1;
1127 /* Computing MIN */
1128 i__2 = l1 + 1;
1129 /* Computing MIN */
1130 i__3 = l1 + 1;
1131 sumr = ddot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc, &b[l1 + min(i__3,*n) * b_dim1], ldb);
1132 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
1133 scaloc = 1.;
1134 a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
1135 da11 = f2c_abs(a11);
1136 if (da11 <= smin)
1137 {
1138 a11 = smin;
1139 da11 = smin;
1140 *info = 1;
1141 }
1142 db = f2c_abs(vec[0]);
1143 if (da11 < 1. && db > 1.)
1144 {
1145 if (db > bignum * da11)
1146 {
1147 scaloc = 1. / db;
1148 }
1149 }
1150 x[0] = vec[0] * scaloc / a11;
1151 if (scaloc != 1.)
1152 {
1153 i__1 = *n;
1154 for (j = 1;
1155 j <= i__1;
1156 ++j)
1157 {
1158 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
1159 /* L190: */
1160 }
1161 *scale *= scaloc;
1162 }
1163 c__[k1 + l1 * c_dim1] = x[0];
1164 }
1165 else if (l1 == l2 && k1 != k2)
1166 {
1167 i__1 = *m - k2;
1168 /* Computing MIN */
1169 i__2 = k2 + 1;
1170 /* Computing MIN */
1171 i__3 = k2 + 1;
1172 suml = ddot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l1 * c_dim1], &c__1);
1173 i__1 = *n - l2;
1174 /* Computing MIN */
1175 i__2 = l2 + 1;
1176 /* Computing MIN */
1177 i__3 = l2 + 1;
1178 sumr = ddot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc, &b[l1 + min(i__3,*n) * b_dim1], ldb);
1179 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
1180 i__1 = *m - k2;
1181 /* Computing MIN */
1182 i__2 = k2 + 1;
1183 /* Computing MIN */
1184 i__3 = k2 + 1;
1185 suml = ddot_(&i__1, &a[k2 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l1 * c_dim1], &c__1);
1186 i__1 = *n - l2;
1187 /* Computing MIN */
1188 i__2 = l2 + 1;
1189 /* Computing MIN */
1190 i__3 = l2 + 1;
1191 sumr = ddot_(&i__1, &c__[k2 + min(i__2,*n) * c_dim1], ldc, &b[l1 + min(i__3,*n) * b_dim1], ldb);
1192 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
1193 d__1 = -sgn * b[l1 + l1 * b_dim1];
1194 dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
1195 if (ierr != 0)
1196 {
1197 *info = 1;
1198 }
1199 if (scaloc != 1.)
1200 {
1201 i__1 = *n;
1202 for (j = 1;
1203 j <= i__1;
1204 ++j)
1205 {
1206 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
1207 /* L200: */
1208 }
1209 *scale *= scaloc;
1210 }
1211 c__[k1 + l1 * c_dim1] = x[0];
1212 c__[k2 + l1 * c_dim1] = x[1];
1213 }
1214 else if (l1 != l2 && k1 == k2)
1215 {
1216 i__1 = *m - k1;
1217 /* Computing MIN */
1218 i__2 = k1 + 1;
1219 /* Computing MIN */
1220 i__3 = k1 + 1;
1221 suml = ddot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l1 * c_dim1], &c__1);
1222 i__1 = *n - l2;
1223 /* Computing MIN */
1224 i__2 = l2 + 1;
1225 /* Computing MIN */
1226 i__3 = l2 + 1;
1227 sumr = ddot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc, &b[l1 + min(i__3,*n) * b_dim1], ldb);
1228 vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * sumr));
1229 i__1 = *m - k1;
1230 /* Computing MIN */
1231 i__2 = k1 + 1;
1232 /* Computing MIN */
1233 i__3 = k1 + 1;
1234 suml = ddot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l2 * c_dim1], &c__1);
1235 i__1 = *n - l2;
1236 /* Computing MIN */
1237 i__2 = l2 + 1;
1238 /* Computing MIN */
1239 i__3 = l2 + 1;
1240 sumr = ddot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc, &b[l2 + min(i__3,*n) * b_dim1], ldb);
1241 vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * sumr));
1242 d__1 = -sgn * a[k1 + k1 * a_dim1];
1243 dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
1244 if (ierr != 0)
1245 {
1246 *info = 1;
1247 }
1248 if (scaloc != 1.)
1249 {
1250 i__1 = *n;
1251 for (j = 1;
1252 j <= i__1;
1253 ++j)
1254 {
1255 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
1256 /* L210: */
1257 }
1258 *scale *= scaloc;
1259 }
1260 c__[k1 + l1 * c_dim1] = x[0];
1261 c__[k1 + l2 * c_dim1] = x[1];
1262 }
1263 else if (l1 != l2 && k1 != k2)
1264 {
1265 i__1 = *m - k2;
1266 /* Computing MIN */
1267 i__2 = k2 + 1;
1268 /* Computing MIN */
1269 i__3 = k2 + 1;
1270 suml = ddot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l1 * c_dim1], &c__1);
1271 i__1 = *n - l2;
1272 /* Computing MIN */
1273 i__2 = l2 + 1;
1274 /* Computing MIN */
1275 i__3 = l2 + 1;
1276 sumr = ddot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc, &b[l1 + min(i__3,*n) * b_dim1], ldb);
1277 vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
1278 i__1 = *m - k2;
1279 /* Computing MIN */
1280 i__2 = k2 + 1;
1281 /* Computing MIN */
1282 i__3 = k2 + 1;
1283 suml = ddot_(&i__1, &a[k1 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l2 * c_dim1], &c__1);
1284 i__1 = *n - l2;
1285 /* Computing MIN */
1286 i__2 = l2 + 1;
1287 /* Computing MIN */
1288 i__3 = l2 + 1;
1289 sumr = ddot_(&i__1, &c__[k1 + min(i__2,*n) * c_dim1], ldc, &b[l2 + min(i__3,*n) * b_dim1], ldb);
1290 vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
1291 i__1 = *m - k2;
1292 /* Computing MIN */
1293 i__2 = k2 + 1;
1294 /* Computing MIN */
1295 i__3 = k2 + 1;
1296 suml = ddot_(&i__1, &a[k2 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l1 * c_dim1], &c__1);
1297 i__1 = *n - l2;
1298 /* Computing MIN */
1299 i__2 = l2 + 1;
1300 /* Computing MIN */
1301 i__3 = l2 + 1;
1302 sumr = ddot_(&i__1, &c__[k2 + min(i__2,*n) * c_dim1], ldc, &b[l1 + min(i__3,*n) * b_dim1], ldb);
1303 vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
1304 i__1 = *m - k2;
1305 /* Computing MIN */
1306 i__2 = k2 + 1;
1307 /* Computing MIN */
1308 i__3 = k2 + 1;
1309 suml = ddot_(&i__1, &a[k2 + min(i__2,*m) * a_dim1], lda, & c__[min(i__3,*m) + l2 * c_dim1], &c__1);
1310 i__1 = *n - l2;
1311 /* Computing MIN */
1312 i__2 = l2 + 1;
1313 /* Computing MIN */
1314 i__3 = l2 + 1;
1315 sumr = ddot_(&i__1, &c__[k2 + min(i__2,*n) * c_dim1], ldc, &b[l2 + min(i__3,*n) * b_dim1], ldb);
1316 vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
1317 dlasy2_(&c_false, &c_true, isgn, &c__2, &c__2, &a[k1 + k1 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, & c__2, &scaloc, x, &c__2, &xnorm, &ierr);
1318 if (ierr != 0)
1319 {
1320 *info = 1;
1321 }
1322 if (scaloc != 1.)
1323 {
1324 i__1 = *n;
1325 for (j = 1;
1326 j <= i__1;
1327 ++j)
1328 {
1329 dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
1330 /* L220: */
1331 }
1332 *scale *= scaloc;
1333 }
1334 c__[k1 + l1 * c_dim1] = x[0];
1335 c__[k1 + l2 * c_dim1] = x[2];
1336 c__[k2 + l1 * c_dim1] = x[1];
1337 c__[k2 + l2 * c_dim1] = x[3];
1338 }
1339 L230:
1340 ;
1341 }
1342 L240:
1343 ;
1344 }
1345 }
1346 return 0;
1347 /* End of DTRSYL */
1348 }
1349 /* dtrsyl_ */
1350