dtrsyl.c
Go to the documentation of this file.
00001 /* dtrsyl.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 static logical c_false = FALSE_;
00020 static integer c__2 = 2;
00021 static doublereal c_b26 = 1.;
00022 static doublereal c_b30 = 0.;
00023 static logical c_true = TRUE_;
00024 
00025 /* Subroutine */ int dtrsyl_(char *trana, char *tranb, integer *isgn, integer 
00026         *m, integer *n, doublereal *a, integer *lda, doublereal *b, integer *
00027         ldb, doublereal *c__, integer *ldc, doublereal *scale, integer *info)
00028 {
00029     /* System generated locals */
00030     integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2, 
00031             i__3, i__4;
00032     doublereal d__1, d__2;
00033 
00034     /* Local variables */
00035     integer j, k, l;
00036     doublereal x[4]     /* was [2][2] */;
00037     integer k1, k2, l1, l2;
00038     doublereal a11, db, da11, vec[4]    /* was [2][2] */, dum[1], eps, sgn;
00039     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
00040             integer *);
00041     integer ierr;
00042     doublereal smin, suml, sumr;
00043     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
00044             integer *);
00045     extern logical lsame_(char *, char *);
00046     integer knext, lnext;
00047     doublereal xnorm;
00048     extern /* Subroutine */ int dlaln2_(logical *, integer *, integer *, 
00049             doublereal *, doublereal *, doublereal *, integer *, doublereal *, 
00050              doublereal *, doublereal *, integer *, doublereal *, doublereal *
00051 , doublereal *, integer *, doublereal *, doublereal *, integer *),
00052              dlasy2_(logical *, logical *, integer *, integer *, integer *, 
00053             doublereal *, integer *, doublereal *, integer *, doublereal *, 
00054             integer *, doublereal *, doublereal *, integer *, doublereal *, 
00055             integer *), dlabad_(doublereal *, doublereal *);
00056     extern doublereal dlamch_(char *), dlange_(char *, integer *, 
00057             integer *, doublereal *, integer *, doublereal *);
00058     doublereal scaloc;
00059     extern /* Subroutine */ int xerbla_(char *, integer *);
00060     doublereal bignum;
00061     logical notrna, notrnb;
00062     doublereal smlnum;
00063 
00064 
00065 /*  -- LAPACK routine (version 3.2) -- */
00066 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00067 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00068 /*     November 2006 */
00069 
00070 /*     .. Scalar Arguments .. */
00071 /*     .. */
00072 /*     .. Array Arguments .. */
00073 /*     .. */
00074 
00075 /*  Purpose */
00076 /*  ======= */
00077 
00078 /*  DTRSYL solves the real Sylvester matrix equation: */
00079 
00080 /*     op(A)*X + X*op(B) = scale*C or */
00081 /*     op(A)*X - X*op(B) = scale*C, */
00082 
00083 /*  where op(A) = A or A**T, and  A and B are both upper quasi- */
00084 /*  triangular. A is M-by-M and B is N-by-N; the right hand side C and */
00085 /*  the solution X are M-by-N; and scale is an output scale factor, set */
00086 /*  <= 1 to avoid overflow in X. */
00087 
00088 /*  A and B must be in Schur canonical form (as returned by DHSEQR), that */
00089 /*  is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; */
00090 /*  each 2-by-2 diagonal block has its diagonal elements equal and its */
00091 /*  off-diagonal elements of opposite sign. */
00092 
00093 /*  Arguments */
00094 /*  ========= */
00095 
00096 /*  TRANA   (input) CHARACTER*1 */
00097 /*          Specifies the option op(A): */
00098 /*          = 'N': op(A) = A    (No transpose) */
00099 /*          = 'T': op(A) = A**T (Transpose) */
00100 /*          = 'C': op(A) = A**H (Conjugate transpose = Transpose) */
00101 
00102 /*  TRANB   (input) CHARACTER*1 */
00103 /*          Specifies the option op(B): */
00104 /*          = 'N': op(B) = B    (No transpose) */
00105 /*          = 'T': op(B) = B**T (Transpose) */
00106 /*          = 'C': op(B) = B**H (Conjugate transpose = Transpose) */
00107 
00108 /*  ISGN    (input) INTEGER */
00109 /*          Specifies the sign in the equation: */
00110 /*          = +1: solve op(A)*X + X*op(B) = scale*C */
00111 /*          = -1: solve op(A)*X - X*op(B) = scale*C */
00112 
00113 /*  M       (input) INTEGER */
00114 /*          The order of the matrix A, and the number of rows in the */
00115 /*          matrices X and C. M >= 0. */
00116 
00117 /*  N       (input) INTEGER */
00118 /*          The order of the matrix B, and the number of columns in the */
00119 /*          matrices X and C. N >= 0. */
00120 
00121 /*  A       (input) DOUBLE PRECISION array, dimension (LDA,M) */
00122 /*          The upper quasi-triangular matrix A, in Schur canonical form. */
00123 
00124 /*  LDA     (input) INTEGER */
00125 /*          The leading dimension of the array A. LDA >= max(1,M). */
00126 
00127 /*  B       (input) DOUBLE PRECISION array, dimension (LDB,N) */
00128 /*          The upper quasi-triangular matrix B, in Schur canonical form. */
00129 
00130 /*  LDB     (input) INTEGER */
00131 /*          The leading dimension of the array B. LDB >= max(1,N). */
00132 
00133 /*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
00134 /*          On entry, the M-by-N right hand side matrix C. */
00135 /*          On exit, C is overwritten by the solution matrix X. */
00136 
00137 /*  LDC     (input) INTEGER */
00138 /*          The leading dimension of the array C. LDC >= max(1,M) */
00139 
00140 /*  SCALE   (output) DOUBLE PRECISION */
00141 /*          The scale factor, scale, set <= 1 to avoid overflow in X. */
00142 
00143 /*  INFO    (output) INTEGER */
00144 /*          = 0: successful exit */
00145 /*          < 0: if INFO = -i, the i-th argument had an illegal value */
00146 /*          = 1: A and B have common or very close eigenvalues; perturbed */
00147 /*               values were used to solve the equation (but the matrices */
00148 /*               A and B are unchanged). */
00149 
00150 /*  ===================================================================== */
00151 
00152 /*     .. Parameters .. */
00153 /*     .. */
00154 /*     .. Local Scalars .. */
00155 /*     .. */
00156 /*     .. Local Arrays .. */
00157 /*     .. */
00158 /*     .. External Functions .. */
00159 /*     .. */
00160 /*     .. External Subroutines .. */
00161 /*     .. */
00162 /*     .. Intrinsic Functions .. */
00163 /*     .. */
00164 /*     .. Executable Statements .. */
00165 
00166 /*     Decode and Test input parameters */
00167 
00168     /* Parameter adjustments */
00169     a_dim1 = *lda;
00170     a_offset = 1 + a_dim1;
00171     a -= a_offset;
00172     b_dim1 = *ldb;
00173     b_offset = 1 + b_dim1;
00174     b -= b_offset;
00175     c_dim1 = *ldc;
00176     c_offset = 1 + c_dim1;
00177     c__ -= c_offset;
00178 
00179     /* Function Body */
00180     notrna = lsame_(trana, "N");
00181     notrnb = lsame_(tranb, "N");
00182 
00183     *info = 0;
00184     if (! notrna && ! lsame_(trana, "T") && ! lsame_(
00185             trana, "C")) {
00186         *info = -1;
00187     } else if (! notrnb && ! lsame_(tranb, "T") && ! 
00188             lsame_(tranb, "C")) {
00189         *info = -2;
00190     } else if (*isgn != 1 && *isgn != -1) {
00191         *info = -3;
00192     } else if (*m < 0) {
00193         *info = -4;
00194     } else if (*n < 0) {
00195         *info = -5;
00196     } else if (*lda < max(1,*m)) {
00197         *info = -7;
00198     } else if (*ldb < max(1,*n)) {
00199         *info = -9;
00200     } else if (*ldc < max(1,*m)) {
00201         *info = -11;
00202     }
00203     if (*info != 0) {
00204         i__1 = -(*info);
00205         xerbla_("DTRSYL", &i__1);
00206         return 0;
00207     }
00208 
00209 /*     Quick return if possible */
00210 
00211     *scale = 1.;
00212     if (*m == 0 || *n == 0) {
00213         return 0;
00214     }
00215 
00216 /*     Set constants to control overflow */
00217 
00218     eps = dlamch_("P");
00219     smlnum = dlamch_("S");
00220     bignum = 1. / smlnum;
00221     dlabad_(&smlnum, &bignum);
00222     smlnum = smlnum * (doublereal) (*m * *n) / eps;
00223     bignum = 1. / smlnum;
00224 
00225 /* Computing MAX */
00226     d__1 = smlnum, d__2 = eps * dlange_("M", m, m, &a[a_offset], lda, dum), d__1 = max(d__1,d__2), d__2 = eps * dlange_("M", n, n, 
00227             &b[b_offset], ldb, dum);
00228     smin = max(d__1,d__2);
00229 
00230     sgn = (doublereal) (*isgn);
00231 
00232     if (notrna && notrnb) {
00233 
00234 /*        Solve    A*X + ISGN*X*B = scale*C. */
00235 
00236 /*        The (K,L)th block of X is determined starting from */
00237 /*        bottom-left corner column by column by */
00238 
00239 /*         A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) */
00240 
00241 /*        Where */
00242 /*                  M                         L-1 */
00243 /*        R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. */
00244 /*                I=K+1                       J=1 */
00245 
00246 /*        Start column loop (index = L) */
00247 /*        L1 (L2) : column index of the first (first) row of X(K,L). */
00248 
00249         lnext = 1;
00250         i__1 = *n;
00251         for (l = 1; l <= i__1; ++l) {
00252             if (l < lnext) {
00253                 goto L60;
00254             }
00255             if (l == *n) {
00256                 l1 = l;
00257                 l2 = l;
00258             } else {
00259                 if (b[l + 1 + l * b_dim1] != 0.) {
00260                     l1 = l;
00261                     l2 = l + 1;
00262                     lnext = l + 2;
00263                 } else {
00264                     l1 = l;
00265                     l2 = l;
00266                     lnext = l + 1;
00267                 }
00268             }
00269 
00270 /*           Start row loop (index = K) */
00271 /*           K1 (K2): row index of the first (last) row of X(K,L). */
00272 
00273             knext = *m;
00274             for (k = *m; k >= 1; --k) {
00275                 if (k > knext) {
00276                     goto L50;
00277                 }
00278                 if (k == 1) {
00279                     k1 = k;
00280                     k2 = k;
00281                 } else {
00282                     if (a[k + (k - 1) * a_dim1] != 0.) {
00283                         k1 = k - 1;
00284                         k2 = k;
00285                         knext = k - 2;
00286                     } else {
00287                         k1 = k;
00288                         k2 = k;
00289                         knext = k - 1;
00290                     }
00291                 }
00292 
00293                 if (l1 == l2 && k1 == k2) {
00294                     i__2 = *m - k1;
00295 /* Computing MIN */
00296                     i__3 = k1 + 1;
00297 /* Computing MIN */
00298                     i__4 = k1 + 1;
00299                     suml = ddot_(&i__2, &a[k1 + min(i__3, *m)* a_dim1], lda, &
00300                             c__[min(i__4, *m)+ l1 * c_dim1], &c__1);
00301                     i__2 = l1 - 1;
00302                     sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * 
00303                             b_dim1 + 1], &c__1);
00304                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00305                     scaloc = 1.;
00306 
00307                     a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
00308                     da11 = abs(a11);
00309                     if (da11 <= smin) {
00310                         a11 = smin;
00311                         da11 = smin;
00312                         *info = 1;
00313                     }
00314                     db = abs(vec[0]);
00315                     if (da11 < 1. && db > 1.) {
00316                         if (db > bignum * da11) {
00317                             scaloc = 1. / db;
00318                         }
00319                     }
00320                     x[0] = vec[0] * scaloc / a11;
00321 
00322                     if (scaloc != 1.) {
00323                         i__2 = *n;
00324                         for (j = 1; j <= i__2; ++j) {
00325                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00326 /* L10: */
00327                         }
00328                         *scale *= scaloc;
00329                     }
00330                     c__[k1 + l1 * c_dim1] = x[0];
00331 
00332                 } else if (l1 == l2 && k1 != k2) {
00333 
00334                     i__2 = *m - k2;
00335 /* Computing MIN */
00336                     i__3 = k2 + 1;
00337 /* Computing MIN */
00338                     i__4 = k2 + 1;
00339                     suml = ddot_(&i__2, &a[k1 + min(i__3, *m)* a_dim1], lda, &
00340                             c__[min(i__4, *m)+ l1 * c_dim1], &c__1);
00341                     i__2 = l1 - 1;
00342                     sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * 
00343                             b_dim1 + 1], &c__1);
00344                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00345 
00346                     i__2 = *m - k2;
00347 /* Computing MIN */
00348                     i__3 = k2 + 1;
00349 /* Computing MIN */
00350                     i__4 = k2 + 1;
00351                     suml = ddot_(&i__2, &a[k2 + min(i__3, *m)* a_dim1], lda, &
00352                             c__[min(i__4, *m)+ l1 * c_dim1], &c__1);
00353                     i__2 = l1 - 1;
00354                     sumr = ddot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l1 * 
00355                             b_dim1 + 1], &c__1);
00356                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
00357 
00358                     d__1 = -sgn * b[l1 + l1 * b_dim1];
00359                     dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 
00360                             * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, 
00361                              &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
00362                     if (ierr != 0) {
00363                         *info = 1;
00364                     }
00365 
00366                     if (scaloc != 1.) {
00367                         i__2 = *n;
00368                         for (j = 1; j <= i__2; ++j) {
00369                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00370 /* L20: */
00371                         }
00372                         *scale *= scaloc;
00373                     }
00374                     c__[k1 + l1 * c_dim1] = x[0];
00375                     c__[k2 + l1 * c_dim1] = x[1];
00376 
00377                 } else if (l1 != l2 && k1 == k2) {
00378 
00379                     i__2 = *m - k1;
00380 /* Computing MIN */
00381                     i__3 = k1 + 1;
00382 /* Computing MIN */
00383                     i__4 = k1 + 1;
00384                     suml = ddot_(&i__2, &a[k1 + min(i__3, *m)* a_dim1], lda, &
00385                             c__[min(i__4, *m)+ l1 * c_dim1], &c__1);
00386                     i__2 = l1 - 1;
00387                     sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * 
00388                             b_dim1 + 1], &c__1);
00389                     vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * 
00390                             sumr));
00391 
00392                     i__2 = *m - k1;
00393 /* Computing MIN */
00394                     i__3 = k1 + 1;
00395 /* Computing MIN */
00396                     i__4 = k1 + 1;
00397                     suml = ddot_(&i__2, &a[k1 + min(i__3, *m)* a_dim1], lda, &
00398                             c__[min(i__4, *m)+ l2 * c_dim1], &c__1);
00399                     i__2 = l1 - 1;
00400                     sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l2 * 
00401                             b_dim1 + 1], &c__1);
00402                     vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * 
00403                             sumr));
00404 
00405                     d__1 = -sgn * a[k1 + k1 * a_dim1];
00406                     dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 *
00407                              b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, 
00408                             &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
00409                     if (ierr != 0) {
00410                         *info = 1;
00411                     }
00412 
00413                     if (scaloc != 1.) {
00414                         i__2 = *n;
00415                         for (j = 1; j <= i__2; ++j) {
00416                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00417 /* L30: */
00418                         }
00419                         *scale *= scaloc;
00420                     }
00421                     c__[k1 + l1 * c_dim1] = x[0];
00422                     c__[k1 + l2 * c_dim1] = x[1];
00423 
00424                 } else if (l1 != l2 && k1 != k2) {
00425 
00426                     i__2 = *m - k2;
00427 /* Computing MIN */
00428                     i__3 = k2 + 1;
00429 /* Computing MIN */
00430                     i__4 = k2 + 1;
00431                     suml = ddot_(&i__2, &a[k1 + min(i__3, *m)* a_dim1], lda, &
00432                             c__[min(i__4, *m)+ l1 * c_dim1], &c__1);
00433                     i__2 = l1 - 1;
00434                     sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l1 * 
00435                             b_dim1 + 1], &c__1);
00436                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00437 
00438                     i__2 = *m - k2;
00439 /* Computing MIN */
00440                     i__3 = k2 + 1;
00441 /* Computing MIN */
00442                     i__4 = k2 + 1;
00443                     suml = ddot_(&i__2, &a[k1 + min(i__3, *m)* a_dim1], lda, &
00444                             c__[min(i__4, *m)+ l2 * c_dim1], &c__1);
00445                     i__2 = l1 - 1;
00446                     sumr = ddot_(&i__2, &c__[k1 + c_dim1], ldc, &b[l2 * 
00447                             b_dim1 + 1], &c__1);
00448                     vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
00449 
00450                     i__2 = *m - k2;
00451 /* Computing MIN */
00452                     i__3 = k2 + 1;
00453 /* Computing MIN */
00454                     i__4 = k2 + 1;
00455                     suml = ddot_(&i__2, &a[k2 + min(i__3, *m)* a_dim1], lda, &
00456                             c__[min(i__4, *m)+ l1 * c_dim1], &c__1);
00457                     i__2 = l1 - 1;
00458                     sumr = ddot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l1 * 
00459                             b_dim1 + 1], &c__1);
00460                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
00461 
00462                     i__2 = *m - k2;
00463 /* Computing MIN */
00464                     i__3 = k2 + 1;
00465 /* Computing MIN */
00466                     i__4 = k2 + 1;
00467                     suml = ddot_(&i__2, &a[k2 + min(i__3, *m)* a_dim1], lda, &
00468                             c__[min(i__4, *m)+ l2 * c_dim1], &c__1);
00469                     i__2 = l1 - 1;
00470                     sumr = ddot_(&i__2, &c__[k2 + c_dim1], ldc, &b[l2 * 
00471                             b_dim1 + 1], &c__1);
00472                     vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
00473 
00474                     dlasy2_(&c_false, &c_false, isgn, &c__2, &c__2, &a[k1 + 
00475                             k1 * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, 
00476                              &c__2, &scaloc, x, &c__2, &xnorm, &ierr);
00477                     if (ierr != 0) {
00478                         *info = 1;
00479                     }
00480 
00481                     if (scaloc != 1.) {
00482                         i__2 = *n;
00483                         for (j = 1; j <= i__2; ++j) {
00484                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00485 /* L40: */
00486                         }
00487                         *scale *= scaloc;
00488                     }
00489                     c__[k1 + l1 * c_dim1] = x[0];
00490                     c__[k1 + l2 * c_dim1] = x[2];
00491                     c__[k2 + l1 * c_dim1] = x[1];
00492                     c__[k2 + l2 * c_dim1] = x[3];
00493                 }
00494 
00495 L50:
00496                 ;
00497             }
00498 
00499 L60:
00500             ;
00501         }
00502 
00503     } else if (! notrna && notrnb) {
00504 
00505 /*        Solve    A' *X + ISGN*X*B = scale*C. */
00506 
00507 /*        The (K,L)th block of X is determined starting from */
00508 /*        upper-left corner column by column by */
00509 
00510 /*          A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) */
00511 
00512 /*        Where */
00513 /*                   K-1                        L-1 */
00514 /*          R(K,L) = SUM [A(I,K)'*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] */
00515 /*                   I=1                        J=1 */
00516 
00517 /*        Start column loop (index = L) */
00518 /*        L1 (L2): column index of the first (last) row of X(K,L) */
00519 
00520         lnext = 1;
00521         i__1 = *n;
00522         for (l = 1; l <= i__1; ++l) {
00523             if (l < lnext) {
00524                 goto L120;
00525             }
00526             if (l == *n) {
00527                 l1 = l;
00528                 l2 = l;
00529             } else {
00530                 if (b[l + 1 + l * b_dim1] != 0.) {
00531                     l1 = l;
00532                     l2 = l + 1;
00533                     lnext = l + 2;
00534                 } else {
00535                     l1 = l;
00536                     l2 = l;
00537                     lnext = l + 1;
00538                 }
00539             }
00540 
00541 /*           Start row loop (index = K) */
00542 /*           K1 (K2): row index of the first (last) row of X(K,L) */
00543 
00544             knext = 1;
00545             i__2 = *m;
00546             for (k = 1; k <= i__2; ++k) {
00547                 if (k < knext) {
00548                     goto L110;
00549                 }
00550                 if (k == *m) {
00551                     k1 = k;
00552                     k2 = k;
00553                 } else {
00554                     if (a[k + 1 + k * a_dim1] != 0.) {
00555                         k1 = k;
00556                         k2 = k + 1;
00557                         knext = k + 2;
00558                     } else {
00559                         k1 = k;
00560                         k2 = k;
00561                         knext = k + 1;
00562                     }
00563                 }
00564 
00565                 if (l1 == l2 && k1 == k2) {
00566                     i__3 = k1 - 1;
00567                     suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00568                             c_dim1 + 1], &c__1);
00569                     i__3 = l1 - 1;
00570                     sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * 
00571                             b_dim1 + 1], &c__1);
00572                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00573                     scaloc = 1.;
00574 
00575                     a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
00576                     da11 = abs(a11);
00577                     if (da11 <= smin) {
00578                         a11 = smin;
00579                         da11 = smin;
00580                         *info = 1;
00581                     }
00582                     db = abs(vec[0]);
00583                     if (da11 < 1. && db > 1.) {
00584                         if (db > bignum * da11) {
00585                             scaloc = 1. / db;
00586                         }
00587                     }
00588                     x[0] = vec[0] * scaloc / a11;
00589 
00590                     if (scaloc != 1.) {
00591                         i__3 = *n;
00592                         for (j = 1; j <= i__3; ++j) {
00593                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00594 /* L70: */
00595                         }
00596                         *scale *= scaloc;
00597                     }
00598                     c__[k1 + l1 * c_dim1] = x[0];
00599 
00600                 } else if (l1 == l2 && k1 != k2) {
00601 
00602                     i__3 = k1 - 1;
00603                     suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00604                             c_dim1 + 1], &c__1);
00605                     i__3 = l1 - 1;
00606                     sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * 
00607                             b_dim1 + 1], &c__1);
00608                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00609 
00610                     i__3 = k1 - 1;
00611                     suml = ddot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * 
00612                             c_dim1 + 1], &c__1);
00613                     i__3 = l1 - 1;
00614                     sumr = ddot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l1 * 
00615                             b_dim1 + 1], &c__1);
00616                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
00617 
00618                     d__1 = -sgn * b[l1 + l1 * b_dim1];
00619                     dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 *
00620                              a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, 
00621                             &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
00622                     if (ierr != 0) {
00623                         *info = 1;
00624                     }
00625 
00626                     if (scaloc != 1.) {
00627                         i__3 = *n;
00628                         for (j = 1; j <= i__3; ++j) {
00629                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00630 /* L80: */
00631                         }
00632                         *scale *= scaloc;
00633                     }
00634                     c__[k1 + l1 * c_dim1] = x[0];
00635                     c__[k2 + l1 * c_dim1] = x[1];
00636 
00637                 } else if (l1 != l2 && k1 == k2) {
00638 
00639                     i__3 = k1 - 1;
00640                     suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00641                             c_dim1 + 1], &c__1);
00642                     i__3 = l1 - 1;
00643                     sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * 
00644                             b_dim1 + 1], &c__1);
00645                     vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * 
00646                             sumr));
00647 
00648                     i__3 = k1 - 1;
00649                     suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * 
00650                             c_dim1 + 1], &c__1);
00651                     i__3 = l1 - 1;
00652                     sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l2 * 
00653                             b_dim1 + 1], &c__1);
00654                     vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * 
00655                             sumr));
00656 
00657                     d__1 = -sgn * a[k1 + k1 * a_dim1];
00658                     dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 *
00659                              b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, 
00660                             &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
00661                     if (ierr != 0) {
00662                         *info = 1;
00663                     }
00664 
00665                     if (scaloc != 1.) {
00666                         i__3 = *n;
00667                         for (j = 1; j <= i__3; ++j) {
00668                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00669 /* L90: */
00670                         }
00671                         *scale *= scaloc;
00672                     }
00673                     c__[k1 + l1 * c_dim1] = x[0];
00674                     c__[k1 + l2 * c_dim1] = x[1];
00675 
00676                 } else if (l1 != l2 && k1 != k2) {
00677 
00678                     i__3 = k1 - 1;
00679                     suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00680                             c_dim1 + 1], &c__1);
00681                     i__3 = l1 - 1;
00682                     sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l1 * 
00683                             b_dim1 + 1], &c__1);
00684                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00685 
00686                     i__3 = k1 - 1;
00687                     suml = ddot_(&i__3, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * 
00688                             c_dim1 + 1], &c__1);
00689                     i__3 = l1 - 1;
00690                     sumr = ddot_(&i__3, &c__[k1 + c_dim1], ldc, &b[l2 * 
00691                             b_dim1 + 1], &c__1);
00692                     vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
00693 
00694                     i__3 = k1 - 1;
00695                     suml = ddot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * 
00696                             c_dim1 + 1], &c__1);
00697                     i__3 = l1 - 1;
00698                     sumr = ddot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l1 * 
00699                             b_dim1 + 1], &c__1);
00700                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
00701 
00702                     i__3 = k1 - 1;
00703                     suml = ddot_(&i__3, &a[k2 * a_dim1 + 1], &c__1, &c__[l2 * 
00704                             c_dim1 + 1], &c__1);
00705                     i__3 = l1 - 1;
00706                     sumr = ddot_(&i__3, &c__[k2 + c_dim1], ldc, &b[l2 * 
00707                             b_dim1 + 1], &c__1);
00708                     vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
00709 
00710                     dlasy2_(&c_true, &c_false, isgn, &c__2, &c__2, &a[k1 + k1 
00711                             * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &
00712                             c__2, &scaloc, x, &c__2, &xnorm, &ierr);
00713                     if (ierr != 0) {
00714                         *info = 1;
00715                     }
00716 
00717                     if (scaloc != 1.) {
00718                         i__3 = *n;
00719                         for (j = 1; j <= i__3; ++j) {
00720                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00721 /* L100: */
00722                         }
00723                         *scale *= scaloc;
00724                     }
00725                     c__[k1 + l1 * c_dim1] = x[0];
00726                     c__[k1 + l2 * c_dim1] = x[2];
00727                     c__[k2 + l1 * c_dim1] = x[1];
00728                     c__[k2 + l2 * c_dim1] = x[3];
00729                 }
00730 
00731 L110:
00732                 ;
00733             }
00734 L120:
00735             ;
00736         }
00737 
00738     } else if (! notrna && ! notrnb) {
00739 
00740 /*        Solve    A'*X + ISGN*X*B' = scale*C. */
00741 
00742 /*        The (K,L)th block of X is determined starting from */
00743 /*        top-right corner column by column by */
00744 
00745 /*           A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L) */
00746 
00747 /*        Where */
00748 /*                     K-1                          N */
00749 /*            R(K,L) = SUM [A(I,K)'*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)']. */
00750 /*                     I=1                        J=L+1 */
00751 
00752 /*        Start column loop (index = L) */
00753 /*        L1 (L2): column index of the first (last) row of X(K,L) */
00754 
00755         lnext = *n;
00756         for (l = *n; l >= 1; --l) {
00757             if (l > lnext) {
00758                 goto L180;
00759             }
00760             if (l == 1) {
00761                 l1 = l;
00762                 l2 = l;
00763             } else {
00764                 if (b[l + (l - 1) * b_dim1] != 0.) {
00765                     l1 = l - 1;
00766                     l2 = l;
00767                     lnext = l - 2;
00768                 } else {
00769                     l1 = l;
00770                     l2 = l;
00771                     lnext = l - 1;
00772                 }
00773             }
00774 
00775 /*           Start row loop (index = K) */
00776 /*           K1 (K2): row index of the first (last) row of X(K,L) */
00777 
00778             knext = 1;
00779             i__1 = *m;
00780             for (k = 1; k <= i__1; ++k) {
00781                 if (k < knext) {
00782                     goto L170;
00783                 }
00784                 if (k == *m) {
00785                     k1 = k;
00786                     k2 = k;
00787                 } else {
00788                     if (a[k + 1 + k * a_dim1] != 0.) {
00789                         k1 = k;
00790                         k2 = k + 1;
00791                         knext = k + 2;
00792                     } else {
00793                         k1 = k;
00794                         k2 = k;
00795                         knext = k + 1;
00796                     }
00797                 }
00798 
00799                 if (l1 == l2 && k1 == k2) {
00800                     i__2 = k1 - 1;
00801                     suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00802                             c_dim1 + 1], &c__1);
00803                     i__2 = *n - l1;
00804 /* Computing MIN */
00805                     i__3 = l1 + 1;
00806 /* Computing MIN */
00807                     i__4 = l1 + 1;
00808                     sumr = ddot_(&i__2, &c__[k1 + min(i__3, *n)* c_dim1], ldc, 
00809                              &b[l1 + min(i__4, *n)* b_dim1], ldb);
00810                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00811                     scaloc = 1.;
00812 
00813                     a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
00814                     da11 = abs(a11);
00815                     if (da11 <= smin) {
00816                         a11 = smin;
00817                         da11 = smin;
00818                         *info = 1;
00819                     }
00820                     db = abs(vec[0]);
00821                     if (da11 < 1. && db > 1.) {
00822                         if (db > bignum * da11) {
00823                             scaloc = 1. / db;
00824                         }
00825                     }
00826                     x[0] = vec[0] * scaloc / a11;
00827 
00828                     if (scaloc != 1.) {
00829                         i__2 = *n;
00830                         for (j = 1; j <= i__2; ++j) {
00831                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00832 /* L130: */
00833                         }
00834                         *scale *= scaloc;
00835                     }
00836                     c__[k1 + l1 * c_dim1] = x[0];
00837 
00838                 } else if (l1 == l2 && k1 != k2) {
00839 
00840                     i__2 = k1 - 1;
00841                     suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00842                             c_dim1 + 1], &c__1);
00843                     i__2 = *n - l2;
00844 /* Computing MIN */
00845                     i__3 = l2 + 1;
00846 /* Computing MIN */
00847                     i__4 = l2 + 1;
00848                     sumr = ddot_(&i__2, &c__[k1 + min(i__3, *n)* c_dim1], ldc, 
00849                              &b[l1 + min(i__4, *n)* b_dim1], ldb);
00850                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00851 
00852                     i__2 = k1 - 1;
00853                     suml = ddot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * 
00854                             c_dim1 + 1], &c__1);
00855                     i__2 = *n - l2;
00856 /* Computing MIN */
00857                     i__3 = l2 + 1;
00858 /* Computing MIN */
00859                     i__4 = l2 + 1;
00860                     sumr = ddot_(&i__2, &c__[k2 + min(i__3, *n)* c_dim1], ldc, 
00861                              &b[l1 + min(i__4, *n)* b_dim1], ldb);
00862                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
00863 
00864                     d__1 = -sgn * b[l1 + l1 * b_dim1];
00865                     dlaln2_(&c_true, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 *
00866                              a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, 
00867                             &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
00868                     if (ierr != 0) {
00869                         *info = 1;
00870                     }
00871 
00872                     if (scaloc != 1.) {
00873                         i__2 = *n;
00874                         for (j = 1; j <= i__2; ++j) {
00875                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00876 /* L140: */
00877                         }
00878                         *scale *= scaloc;
00879                     }
00880                     c__[k1 + l1 * c_dim1] = x[0];
00881                     c__[k2 + l1 * c_dim1] = x[1];
00882 
00883                 } else if (l1 != l2 && k1 == k2) {
00884 
00885                     i__2 = k1 - 1;
00886                     suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00887                             c_dim1 + 1], &c__1);
00888                     i__2 = *n - l2;
00889 /* Computing MIN */
00890                     i__3 = l2 + 1;
00891 /* Computing MIN */
00892                     i__4 = l2 + 1;
00893                     sumr = ddot_(&i__2, &c__[k1 + min(i__3, *n)* c_dim1], ldc, 
00894                              &b[l1 + min(i__4, *n)* b_dim1], ldb);
00895                     vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * 
00896                             sumr));
00897 
00898                     i__2 = k1 - 1;
00899                     suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * 
00900                             c_dim1 + 1], &c__1);
00901                     i__2 = *n - l2;
00902 /* Computing MIN */
00903                     i__3 = l2 + 1;
00904 /* Computing MIN */
00905                     i__4 = l2 + 1;
00906                     sumr = ddot_(&i__2, &c__[k1 + min(i__3, *n)* c_dim1], ldc, 
00907                              &b[l2 + min(i__4, *n)* b_dim1], ldb);
00908                     vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * 
00909                             sumr));
00910 
00911                     d__1 = -sgn * a[k1 + k1 * a_dim1];
00912                     dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 
00913                             * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, 
00914                              &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
00915                     if (ierr != 0) {
00916                         *info = 1;
00917                     }
00918 
00919                     if (scaloc != 1.) {
00920                         i__2 = *n;
00921                         for (j = 1; j <= i__2; ++j) {
00922                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00923 /* L150: */
00924                         }
00925                         *scale *= scaloc;
00926                     }
00927                     c__[k1 + l1 * c_dim1] = x[0];
00928                     c__[k1 + l2 * c_dim1] = x[1];
00929 
00930                 } else if (l1 != l2 && k1 != k2) {
00931 
00932                     i__2 = k1 - 1;
00933                     suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l1 * 
00934                             c_dim1 + 1], &c__1);
00935                     i__2 = *n - l2;
00936 /* Computing MIN */
00937                     i__3 = l2 + 1;
00938 /* Computing MIN */
00939                     i__4 = l2 + 1;
00940                     sumr = ddot_(&i__2, &c__[k1 + min(i__3, *n)* c_dim1], ldc, 
00941                              &b[l1 + min(i__4, *n)* b_dim1], ldb);
00942                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
00943 
00944                     i__2 = k1 - 1;
00945                     suml = ddot_(&i__2, &a[k1 * a_dim1 + 1], &c__1, &c__[l2 * 
00946                             c_dim1 + 1], &c__1);
00947                     i__2 = *n - l2;
00948 /* Computing MIN */
00949                     i__3 = l2 + 1;
00950 /* Computing MIN */
00951                     i__4 = l2 + 1;
00952                     sumr = ddot_(&i__2, &c__[k1 + min(i__3, *n)* c_dim1], ldc, 
00953                              &b[l2 + min(i__4, *n)* b_dim1], ldb);
00954                     vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
00955 
00956                     i__2 = k1 - 1;
00957                     suml = ddot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l1 * 
00958                             c_dim1 + 1], &c__1);
00959                     i__2 = *n - l2;
00960 /* Computing MIN */
00961                     i__3 = l2 + 1;
00962 /* Computing MIN */
00963                     i__4 = l2 + 1;
00964                     sumr = ddot_(&i__2, &c__[k2 + min(i__3, *n)* c_dim1], ldc, 
00965                              &b[l1 + min(i__4, *n)* b_dim1], ldb);
00966                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
00967 
00968                     i__2 = k1 - 1;
00969                     suml = ddot_(&i__2, &a[k2 * a_dim1 + 1], &c__1, &c__[l2 * 
00970                             c_dim1 + 1], &c__1);
00971                     i__2 = *n - l2;
00972 /* Computing MIN */
00973                     i__3 = l2 + 1;
00974 /* Computing MIN */
00975                     i__4 = l2 + 1;
00976                     sumr = ddot_(&i__2, &c__[k2 + min(i__3, *n)* c_dim1], ldc, 
00977                              &b[l2 + min(i__4, *n)* b_dim1], ldb);
00978                     vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
00979 
00980                     dlasy2_(&c_true, &c_true, isgn, &c__2, &c__2, &a[k1 + k1 *
00981                              a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &
00982                             c__2, &scaloc, x, &c__2, &xnorm, &ierr);
00983                     if (ierr != 0) {
00984                         *info = 1;
00985                     }
00986 
00987                     if (scaloc != 1.) {
00988                         i__2 = *n;
00989                         for (j = 1; j <= i__2; ++j) {
00990                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
00991 /* L160: */
00992                         }
00993                         *scale *= scaloc;
00994                     }
00995                     c__[k1 + l1 * c_dim1] = x[0];
00996                     c__[k1 + l2 * c_dim1] = x[2];
00997                     c__[k2 + l1 * c_dim1] = x[1];
00998                     c__[k2 + l2 * c_dim1] = x[3];
00999                 }
01000 
01001 L170:
01002                 ;
01003             }
01004 L180:
01005             ;
01006         }
01007 
01008     } else if (notrna && ! notrnb) {
01009 
01010 /*        Solve    A*X + ISGN*X*B' = scale*C. */
01011 
01012 /*        The (K,L)th block of X is determined starting from */
01013 /*        bottom-right corner column by column by */
01014 
01015 /*            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L) */
01016 
01017 /*        Where */
01018 /*                      M                          N */
01019 /*            R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)']. */
01020 /*                    I=K+1                      J=L+1 */
01021 
01022 /*        Start column loop (index = L) */
01023 /*        L1 (L2): column index of the first (last) row of X(K,L) */
01024 
01025         lnext = *n;
01026         for (l = *n; l >= 1; --l) {
01027             if (l > lnext) {
01028                 goto L240;
01029             }
01030             if (l == 1) {
01031                 l1 = l;
01032                 l2 = l;
01033             } else {
01034                 if (b[l + (l - 1) * b_dim1] != 0.) {
01035                     l1 = l - 1;
01036                     l2 = l;
01037                     lnext = l - 2;
01038                 } else {
01039                     l1 = l;
01040                     l2 = l;
01041                     lnext = l - 1;
01042                 }
01043             }
01044 
01045 /*           Start row loop (index = K) */
01046 /*           K1 (K2): row index of the first (last) row of X(K,L) */
01047 
01048             knext = *m;
01049             for (k = *m; k >= 1; --k) {
01050                 if (k > knext) {
01051                     goto L230;
01052                 }
01053                 if (k == 1) {
01054                     k1 = k;
01055                     k2 = k;
01056                 } else {
01057                     if (a[k + (k - 1) * a_dim1] != 0.) {
01058                         k1 = k - 1;
01059                         k2 = k;
01060                         knext = k - 2;
01061                     } else {
01062                         k1 = k;
01063                         k2 = k;
01064                         knext = k - 1;
01065                     }
01066                 }
01067 
01068                 if (l1 == l2 && k1 == k2) {
01069                     i__1 = *m - k1;
01070 /* Computing MIN */
01071                     i__2 = k1 + 1;
01072 /* Computing MIN */
01073                     i__3 = k1 + 1;
01074                     suml = ddot_(&i__1, &a[k1 + min(i__2, *m)* a_dim1], lda, &
01075                             c__[min(i__3, *m)+ l1 * c_dim1], &c__1);
01076                     i__1 = *n - l1;
01077 /* Computing MIN */
01078                     i__2 = l1 + 1;
01079 /* Computing MIN */
01080                     i__3 = l1 + 1;
01081                     sumr = ddot_(&i__1, &c__[k1 + min(i__2, *n)* c_dim1], ldc, 
01082                              &b[l1 + min(i__3, *n)* b_dim1], ldb);
01083                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
01084                     scaloc = 1.;
01085 
01086                     a11 = a[k1 + k1 * a_dim1] + sgn * b[l1 + l1 * b_dim1];
01087                     da11 = abs(a11);
01088                     if (da11 <= smin) {
01089                         a11 = smin;
01090                         da11 = smin;
01091                         *info = 1;
01092                     }
01093                     db = abs(vec[0]);
01094                     if (da11 < 1. && db > 1.) {
01095                         if (db > bignum * da11) {
01096                             scaloc = 1. / db;
01097                         }
01098                     }
01099                     x[0] = vec[0] * scaloc / a11;
01100 
01101                     if (scaloc != 1.) {
01102                         i__1 = *n;
01103                         for (j = 1; j <= i__1; ++j) {
01104                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
01105 /* L190: */
01106                         }
01107                         *scale *= scaloc;
01108                     }
01109                     c__[k1 + l1 * c_dim1] = x[0];
01110 
01111                 } else if (l1 == l2 && k1 != k2) {
01112 
01113                     i__1 = *m - k2;
01114 /* Computing MIN */
01115                     i__2 = k2 + 1;
01116 /* Computing MIN */
01117                     i__3 = k2 + 1;
01118                     suml = ddot_(&i__1, &a[k1 + min(i__2, *m)* a_dim1], lda, &
01119                             c__[min(i__3, *m)+ l1 * c_dim1], &c__1);
01120                     i__1 = *n - l2;
01121 /* Computing MIN */
01122                     i__2 = l2 + 1;
01123 /* Computing MIN */
01124                     i__3 = l2 + 1;
01125                     sumr = ddot_(&i__1, &c__[k1 + min(i__2, *n)* c_dim1], ldc, 
01126                              &b[l1 + min(i__3, *n)* b_dim1], ldb);
01127                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
01128 
01129                     i__1 = *m - k2;
01130 /* Computing MIN */
01131                     i__2 = k2 + 1;
01132 /* Computing MIN */
01133                     i__3 = k2 + 1;
01134                     suml = ddot_(&i__1, &a[k2 + min(i__2, *m)* a_dim1], lda, &
01135                             c__[min(i__3, *m)+ l1 * c_dim1], &c__1);
01136                     i__1 = *n - l2;
01137 /* Computing MIN */
01138                     i__2 = l2 + 1;
01139 /* Computing MIN */
01140                     i__3 = l2 + 1;
01141                     sumr = ddot_(&i__1, &c__[k2 + min(i__2, *n)* c_dim1], ldc, 
01142                              &b[l1 + min(i__3, *n)* b_dim1], ldb);
01143                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
01144 
01145                     d__1 = -sgn * b[l1 + l1 * b_dim1];
01146                     dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &a[k1 + k1 
01147                             * a_dim1], lda, &c_b26, &c_b26, vec, &c__2, &d__1, 
01148                              &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
01149                     if (ierr != 0) {
01150                         *info = 1;
01151                     }
01152 
01153                     if (scaloc != 1.) {
01154                         i__1 = *n;
01155                         for (j = 1; j <= i__1; ++j) {
01156                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
01157 /* L200: */
01158                         }
01159                         *scale *= scaloc;
01160                     }
01161                     c__[k1 + l1 * c_dim1] = x[0];
01162                     c__[k2 + l1 * c_dim1] = x[1];
01163 
01164                 } else if (l1 != l2 && k1 == k2) {
01165 
01166                     i__1 = *m - k1;
01167 /* Computing MIN */
01168                     i__2 = k1 + 1;
01169 /* Computing MIN */
01170                     i__3 = k1 + 1;
01171                     suml = ddot_(&i__1, &a[k1 + min(i__2, *m)* a_dim1], lda, &
01172                             c__[min(i__3, *m)+ l1 * c_dim1], &c__1);
01173                     i__1 = *n - l2;
01174 /* Computing MIN */
01175                     i__2 = l2 + 1;
01176 /* Computing MIN */
01177                     i__3 = l2 + 1;
01178                     sumr = ddot_(&i__1, &c__[k1 + min(i__2, *n)* c_dim1], ldc, 
01179                              &b[l1 + min(i__3, *n)* b_dim1], ldb);
01180                     vec[0] = sgn * (c__[k1 + l1 * c_dim1] - (suml + sgn * 
01181                             sumr));
01182 
01183                     i__1 = *m - k1;
01184 /* Computing MIN */
01185                     i__2 = k1 + 1;
01186 /* Computing MIN */
01187                     i__3 = k1 + 1;
01188                     suml = ddot_(&i__1, &a[k1 + min(i__2, *m)* a_dim1], lda, &
01189                             c__[min(i__3, *m)+ l2 * c_dim1], &c__1);
01190                     i__1 = *n - l2;
01191 /* Computing MIN */
01192                     i__2 = l2 + 1;
01193 /* Computing MIN */
01194                     i__3 = l2 + 1;
01195                     sumr = ddot_(&i__1, &c__[k1 + min(i__2, *n)* c_dim1], ldc, 
01196                              &b[l2 + min(i__3, *n)* b_dim1], ldb);
01197                     vec[1] = sgn * (c__[k1 + l2 * c_dim1] - (suml + sgn * 
01198                             sumr));
01199 
01200                     d__1 = -sgn * a[k1 + k1 * a_dim1];
01201                     dlaln2_(&c_false, &c__2, &c__1, &smin, &c_b26, &b[l1 + l1 
01202                             * b_dim1], ldb, &c_b26, &c_b26, vec, &c__2, &d__1, 
01203                              &c_b30, x, &c__2, &scaloc, &xnorm, &ierr);
01204                     if (ierr != 0) {
01205                         *info = 1;
01206                     }
01207 
01208                     if (scaloc != 1.) {
01209                         i__1 = *n;
01210                         for (j = 1; j <= i__1; ++j) {
01211                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
01212 /* L210: */
01213                         }
01214                         *scale *= scaloc;
01215                     }
01216                     c__[k1 + l1 * c_dim1] = x[0];
01217                     c__[k1 + l2 * c_dim1] = x[1];
01218 
01219                 } else if (l1 != l2 && k1 != k2) {
01220 
01221                     i__1 = *m - k2;
01222 /* Computing MIN */
01223                     i__2 = k2 + 1;
01224 /* Computing MIN */
01225                     i__3 = k2 + 1;
01226                     suml = ddot_(&i__1, &a[k1 + min(i__2, *m)* a_dim1], lda, &
01227                             c__[min(i__3, *m)+ l1 * c_dim1], &c__1);
01228                     i__1 = *n - l2;
01229 /* Computing MIN */
01230                     i__2 = l2 + 1;
01231 /* Computing MIN */
01232                     i__3 = l2 + 1;
01233                     sumr = ddot_(&i__1, &c__[k1 + min(i__2, *n)* c_dim1], ldc, 
01234                              &b[l1 + min(i__3, *n)* b_dim1], ldb);
01235                     vec[0] = c__[k1 + l1 * c_dim1] - (suml + sgn * sumr);
01236 
01237                     i__1 = *m - k2;
01238 /* Computing MIN */
01239                     i__2 = k2 + 1;
01240 /* Computing MIN */
01241                     i__3 = k2 + 1;
01242                     suml = ddot_(&i__1, &a[k1 + min(i__2, *m)* a_dim1], lda, &
01243                             c__[min(i__3, *m)+ l2 * c_dim1], &c__1);
01244                     i__1 = *n - l2;
01245 /* Computing MIN */
01246                     i__2 = l2 + 1;
01247 /* Computing MIN */
01248                     i__3 = l2 + 1;
01249                     sumr = ddot_(&i__1, &c__[k1 + min(i__2, *n)* c_dim1], ldc, 
01250                              &b[l2 + min(i__3, *n)* b_dim1], ldb);
01251                     vec[2] = c__[k1 + l2 * c_dim1] - (suml + sgn * sumr);
01252 
01253                     i__1 = *m - k2;
01254 /* Computing MIN */
01255                     i__2 = k2 + 1;
01256 /* Computing MIN */
01257                     i__3 = k2 + 1;
01258                     suml = ddot_(&i__1, &a[k2 + min(i__2, *m)* a_dim1], lda, &
01259                             c__[min(i__3, *m)+ l1 * c_dim1], &c__1);
01260                     i__1 = *n - l2;
01261 /* Computing MIN */
01262                     i__2 = l2 + 1;
01263 /* Computing MIN */
01264                     i__3 = l2 + 1;
01265                     sumr = ddot_(&i__1, &c__[k2 + min(i__2, *n)* c_dim1], ldc, 
01266                              &b[l1 + min(i__3, *n)* b_dim1], ldb);
01267                     vec[1] = c__[k2 + l1 * c_dim1] - (suml + sgn * sumr);
01268 
01269                     i__1 = *m - k2;
01270 /* Computing MIN */
01271                     i__2 = k2 + 1;
01272 /* Computing MIN */
01273                     i__3 = k2 + 1;
01274                     suml = ddot_(&i__1, &a[k2 + min(i__2, *m)* a_dim1], lda, &
01275                             c__[min(i__3, *m)+ l2 * c_dim1], &c__1);
01276                     i__1 = *n - l2;
01277 /* Computing MIN */
01278                     i__2 = l2 + 1;
01279 /* Computing MIN */
01280                     i__3 = l2 + 1;
01281                     sumr = ddot_(&i__1, &c__[k2 + min(i__2, *n)* c_dim1], ldc, 
01282                              &b[l2 + min(i__3, *n)* b_dim1], ldb);
01283                     vec[3] = c__[k2 + l2 * c_dim1] - (suml + sgn * sumr);
01284 
01285                     dlasy2_(&c_false, &c_true, isgn, &c__2, &c__2, &a[k1 + k1 
01286                             * a_dim1], lda, &b[l1 + l1 * b_dim1], ldb, vec, &
01287                             c__2, &scaloc, x, &c__2, &xnorm, &ierr);
01288                     if (ierr != 0) {
01289                         *info = 1;
01290                     }
01291 
01292                     if (scaloc != 1.) {
01293                         i__1 = *n;
01294                         for (j = 1; j <= i__1; ++j) {
01295                             dscal_(m, &scaloc, &c__[j * c_dim1 + 1], &c__1);
01296 /* L220: */
01297                         }
01298                         *scale *= scaloc;
01299                     }
01300                     c__[k1 + l1 * c_dim1] = x[0];
01301                     c__[k1 + l2 * c_dim1] = x[2];
01302                     c__[k2 + l1 * c_dim1] = x[1];
01303                     c__[k2 + l2 * c_dim1] = x[3];
01304                 }
01305 
01306 L230:
01307                 ;
01308             }
01309 L240:
01310             ;
01311         }
01312 
01313     }
01314 
01315     return 0;
01316 
01317 /*     End of DTRSYL */
01318 
01319 } /* dtrsyl_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:50