slasy2.c
Go to the documentation of this file.
00001 /* slasy2.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__4 = 4;
00019 static integer c__1 = 1;
00020 static integer c__16 = 16;
00021 static integer c__0 = 0;
00022 
00023 /* Subroutine */ int slasy2_(logical *ltranl, logical *ltranr, integer *isgn, 
00024         integer *n1, integer *n2, real *tl, integer *ldtl, real *tr, integer *
00025         ldtr, real *b, integer *ldb, real *scale, real *x, integer *ldx, real 
00026         *xnorm, integer *info)
00027 {
00028     /* Initialized data */
00029 
00030     static integer locu12[4] = { 3,4,1,2 };
00031     static integer locl21[4] = { 2,1,4,3 };
00032     static integer locu22[4] = { 4,3,2,1 };
00033     static logical xswpiv[4] = { FALSE_,FALSE_,TRUE_,TRUE_ };
00034     static logical bswpiv[4] = { FALSE_,TRUE_,FALSE_,TRUE_ };
00035 
00036     /* System generated locals */
00037     integer b_dim1, b_offset, tl_dim1, tl_offset, tr_dim1, tr_offset, x_dim1, 
00038             x_offset;
00039     real r__1, r__2, r__3, r__4, r__5, r__6, r__7, r__8;
00040 
00041     /* Local variables */
00042     integer i__, j, k;
00043     real x2[2], l21, u11, u12;
00044     integer ip, jp;
00045     real u22, t16[16]   /* was [4][4] */, gam, bet, eps, sgn, tmp[4], tau1, 
00046             btmp[4], smin;
00047     integer ipiv;
00048     real temp;
00049     integer jpiv[4];
00050     real xmax;
00051     integer ipsv, jpsv;
00052     logical bswap;
00053     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
00054             integer *), sswap_(integer *, real *, integer *, real *, integer *
00055 );
00056     logical xswap;
00057     extern doublereal slamch_(char *);
00058     extern integer isamax_(integer *, real *, integer *);
00059     real smlnum;
00060 
00061 
00062 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00063 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00064 /*     November 2006 */
00065 
00066 /*     .. Scalar Arguments .. */
00067 /*     .. */
00068 /*     .. Array Arguments .. */
00069 /*     .. */
00070 
00071 /*  Purpose */
00072 /*  ======= */
00073 
00074 /*  SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in */
00075 
00076 /*         op(TL)*X + ISGN*X*op(TR) = SCALE*B, */
00077 
00078 /*  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or */
00079 /*  -1.  op(T) = T or T', where T' denotes the transpose of T. */
00080 
00081 /*  Arguments */
00082 /*  ========= */
00083 
00084 /*  LTRANL  (input) LOGICAL */
00085 /*          On entry, LTRANL specifies the op(TL): */
00086 /*             = .FALSE., op(TL) = TL, */
00087 /*             = .TRUE., op(TL) = TL'. */
00088 
00089 /*  LTRANR  (input) LOGICAL */
00090 /*          On entry, LTRANR specifies the op(TR): */
00091 /*            = .FALSE., op(TR) = TR, */
00092 /*            = .TRUE., op(TR) = TR'. */
00093 
00094 /*  ISGN    (input) INTEGER */
00095 /*          On entry, ISGN specifies the sign of the equation */
00096 /*          as described before. ISGN may only be 1 or -1. */
00097 
00098 /*  N1      (input) INTEGER */
00099 /*          On entry, N1 specifies the order of matrix TL. */
00100 /*          N1 may only be 0, 1 or 2. */
00101 
00102 /*  N2      (input) INTEGER */
00103 /*          On entry, N2 specifies the order of matrix TR. */
00104 /*          N2 may only be 0, 1 or 2. */
00105 
00106 /*  TL      (input) REAL array, dimension (LDTL,2) */
00107 /*          On entry, TL contains an N1 by N1 matrix. */
00108 
00109 /*  LDTL    (input) INTEGER */
00110 /*          The leading dimension of the matrix TL. LDTL >= max(1,N1). */
00111 
00112 /*  TR      (input) REAL array, dimension (LDTR,2) */
00113 /*          On entry, TR contains an N2 by N2 matrix. */
00114 
00115 /*  LDTR    (input) INTEGER */
00116 /*          The leading dimension of the matrix TR. LDTR >= max(1,N2). */
00117 
00118 /*  B       (input) REAL array, dimension (LDB,2) */
00119 /*          On entry, the N1 by N2 matrix B contains the right-hand */
00120 /*          side of the equation. */
00121 
00122 /*  LDB     (input) INTEGER */
00123 /*          The leading dimension of the matrix B. LDB >= max(1,N1). */
00124 
00125 /*  SCALE   (output) REAL */
00126 /*          On exit, SCALE contains the scale factor. SCALE is chosen */
00127 /*          less than or equal to 1 to prevent the solution overflowing. */
00128 
00129 /*  X       (output) REAL array, dimension (LDX,2) */
00130 /*          On exit, X contains the N1 by N2 solution. */
00131 
00132 /*  LDX     (input) INTEGER */
00133 /*          The leading dimension of the matrix X. LDX >= max(1,N1). */
00134 
00135 /*  XNORM   (output) REAL */
00136 /*          On exit, XNORM is the infinity-norm of the solution. */
00137 
00138 /*  INFO    (output) INTEGER */
00139 /*          On exit, INFO is set to */
00140 /*             0: successful exit. */
00141 /*             1: TL and TR have too close eigenvalues, so TL or */
00142 /*                TR is perturbed to get a nonsingular equation. */
00143 /*          NOTE: In the interests of speed, this routine does not */
00144 /*                check the inputs for errors. */
00145 
00146 /* ===================================================================== */
00147 
00148 /*     .. Parameters .. */
00149 /*     .. */
00150 /*     .. Local Scalars .. */
00151 /*     .. */
00152 /*     .. Local Arrays .. */
00153 /*     .. */
00154 /*     .. External Functions .. */
00155 /*     .. */
00156 /*     .. External Subroutines .. */
00157 /*     .. */
00158 /*     .. Intrinsic Functions .. */
00159 /*     .. */
00160 /*     .. Data statements .. */
00161     /* Parameter adjustments */
00162     tl_dim1 = *ldtl;
00163     tl_offset = 1 + tl_dim1;
00164     tl -= tl_offset;
00165     tr_dim1 = *ldtr;
00166     tr_offset = 1 + tr_dim1;
00167     tr -= tr_offset;
00168     b_dim1 = *ldb;
00169     b_offset = 1 + b_dim1;
00170     b -= b_offset;
00171     x_dim1 = *ldx;
00172     x_offset = 1 + x_dim1;
00173     x -= x_offset;
00174 
00175     /* Function Body */
00176 /*     .. */
00177 /*     .. Executable Statements .. */
00178 
00179 /*     Do not check the input parameters for errors */
00180 
00181     *info = 0;
00182 
00183 /*     Quick return if possible */
00184 
00185     if (*n1 == 0 || *n2 == 0) {
00186         return 0;
00187     }
00188 
00189 /*     Set constants to control overflow */
00190 
00191     eps = slamch_("P");
00192     smlnum = slamch_("S") / eps;
00193     sgn = (real) (*isgn);
00194 
00195     k = *n1 + *n1 + *n2 - 2;
00196     switch (k) {
00197         case 1:  goto L10;
00198         case 2:  goto L20;
00199         case 3:  goto L30;
00200         case 4:  goto L50;
00201     }
00202 
00203 /*     1 by 1: TL11*X + SGN*X*TR11 = B11 */
00204 
00205 L10:
00206     tau1 = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1];
00207     bet = dabs(tau1);
00208     if (bet <= smlnum) {
00209         tau1 = smlnum;
00210         bet = smlnum;
00211         *info = 1;
00212     }
00213 
00214     *scale = 1.f;
00215     gam = (r__1 = b[b_dim1 + 1], dabs(r__1));
00216     if (smlnum * gam > bet) {
00217         *scale = 1.f / gam;
00218     }
00219 
00220     x[x_dim1 + 1] = b[b_dim1 + 1] * *scale / tau1;
00221     *xnorm = (r__1 = x[x_dim1 + 1], dabs(r__1));
00222     return 0;
00223 
00224 /*     1 by 2: */
00225 /*     TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12]  = [B11 B12] */
00226 /*                                       [TR21 TR22] */
00227 
00228 L20:
00229 
00230 /* Computing MAX */
00231 /* Computing MAX */
00232     r__7 = (r__1 = tl[tl_dim1 + 1], dabs(r__1)), r__8 = (r__2 = tr[tr_dim1 + 
00233             1], dabs(r__2)), r__7 = max(r__7,r__8), r__8 = (r__3 = tr[(
00234             tr_dim1 << 1) + 1], dabs(r__3)), r__7 = max(r__7,r__8), r__8 = (
00235             r__4 = tr[tr_dim1 + 2], dabs(r__4)), r__7 = max(r__7,r__8), r__8 =
00236              (r__5 = tr[(tr_dim1 << 1) + 2], dabs(r__5));
00237     r__6 = eps * dmax(r__7,r__8);
00238     smin = dmax(r__6,smlnum);
00239     tmp[0] = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1];
00240     tmp[3] = tl[tl_dim1 + 1] + sgn * tr[(tr_dim1 << 1) + 2];
00241     if (*ltranr) {
00242         tmp[1] = sgn * tr[tr_dim1 + 2];
00243         tmp[2] = sgn * tr[(tr_dim1 << 1) + 1];
00244     } else {
00245         tmp[1] = sgn * tr[(tr_dim1 << 1) + 1];
00246         tmp[2] = sgn * tr[tr_dim1 + 2];
00247     }
00248     btmp[0] = b[b_dim1 + 1];
00249     btmp[1] = b[(b_dim1 << 1) + 1];
00250     goto L40;
00251 
00252 /*     2 by 1: */
00253 /*          op[TL11 TL12]*[X11] + ISGN* [X11]*TR11  = [B11] */
00254 /*            [TL21 TL22] [X21]         [X21]         [B21] */
00255 
00256 L30:
00257 /* Computing MAX */
00258 /* Computing MAX */
00259     r__7 = (r__1 = tr[tr_dim1 + 1], dabs(r__1)), r__8 = (r__2 = tl[tl_dim1 + 
00260             1], dabs(r__2)), r__7 = max(r__7,r__8), r__8 = (r__3 = tl[(
00261             tl_dim1 << 1) + 1], dabs(r__3)), r__7 = max(r__7,r__8), r__8 = (
00262             r__4 = tl[tl_dim1 + 2], dabs(r__4)), r__7 = max(r__7,r__8), r__8 =
00263              (r__5 = tl[(tl_dim1 << 1) + 2], dabs(r__5));
00264     r__6 = eps * dmax(r__7,r__8);
00265     smin = dmax(r__6,smlnum);
00266     tmp[0] = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1];
00267     tmp[3] = tl[(tl_dim1 << 1) + 2] + sgn * tr[tr_dim1 + 1];
00268     if (*ltranl) {
00269         tmp[1] = tl[(tl_dim1 << 1) + 1];
00270         tmp[2] = tl[tl_dim1 + 2];
00271     } else {
00272         tmp[1] = tl[tl_dim1 + 2];
00273         tmp[2] = tl[(tl_dim1 << 1) + 1];
00274     }
00275     btmp[0] = b[b_dim1 + 1];
00276     btmp[1] = b[b_dim1 + 2];
00277 L40:
00278 
00279 /*     Solve 2 by 2 system using complete pivoting. */
00280 /*     Set pivots less than SMIN to SMIN. */
00281 
00282     ipiv = isamax_(&c__4, tmp, &c__1);
00283     u11 = tmp[ipiv - 1];
00284     if (dabs(u11) <= smin) {
00285         *info = 1;
00286         u11 = smin;
00287     }
00288     u12 = tmp[locu12[ipiv - 1] - 1];
00289     l21 = tmp[locl21[ipiv - 1] - 1] / u11;
00290     u22 = tmp[locu22[ipiv - 1] - 1] - u12 * l21;
00291     xswap = xswpiv[ipiv - 1];
00292     bswap = bswpiv[ipiv - 1];
00293     if (dabs(u22) <= smin) {
00294         *info = 1;
00295         u22 = smin;
00296     }
00297     if (bswap) {
00298         temp = btmp[1];
00299         btmp[1] = btmp[0] - l21 * temp;
00300         btmp[0] = temp;
00301     } else {
00302         btmp[1] -= l21 * btmp[0];
00303     }
00304     *scale = 1.f;
00305     if (smlnum * 2.f * dabs(btmp[1]) > dabs(u22) || smlnum * 2.f * dabs(btmp[
00306             0]) > dabs(u11)) {
00307 /* Computing MAX */
00308         r__1 = dabs(btmp[0]), r__2 = dabs(btmp[1]);
00309         *scale = .5f / dmax(r__1,r__2);
00310         btmp[0] *= *scale;
00311         btmp[1] *= *scale;
00312     }
00313     x2[1] = btmp[1] / u22;
00314     x2[0] = btmp[0] / u11 - u12 / u11 * x2[1];
00315     if (xswap) {
00316         temp = x2[1];
00317         x2[1] = x2[0];
00318         x2[0] = temp;
00319     }
00320     x[x_dim1 + 1] = x2[0];
00321     if (*n1 == 1) {
00322         x[(x_dim1 << 1) + 1] = x2[1];
00323         *xnorm = (r__1 = x[x_dim1 + 1], dabs(r__1)) + (r__2 = x[(x_dim1 << 1) 
00324                 + 1], dabs(r__2));
00325     } else {
00326         x[x_dim1 + 2] = x2[1];
00327 /* Computing MAX */
00328         r__3 = (r__1 = x[x_dim1 + 1], dabs(r__1)), r__4 = (r__2 = x[x_dim1 + 
00329                 2], dabs(r__2));
00330         *xnorm = dmax(r__3,r__4);
00331     }
00332     return 0;
00333 
00334 /*     2 by 2: */
00335 /*     op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12] */
00336 /*       [TL21 TL22] [X21 X22]        [X21 X22]   [TR21 TR22]   [B21 B22] */
00337 
00338 /*     Solve equivalent 4 by 4 system using complete pivoting. */
00339 /*     Set pivots less than SMIN to SMIN. */
00340 
00341 L50:
00342 /* Computing MAX */
00343     r__5 = (r__1 = tr[tr_dim1 + 1], dabs(r__1)), r__6 = (r__2 = tr[(tr_dim1 <<
00344              1) + 1], dabs(r__2)), r__5 = max(r__5,r__6), r__6 = (r__3 = tr[
00345             tr_dim1 + 2], dabs(r__3)), r__5 = max(r__5,r__6), r__6 = (r__4 = 
00346             tr[(tr_dim1 << 1) + 2], dabs(r__4));
00347     smin = dmax(r__5,r__6);
00348 /* Computing MAX */
00349     r__5 = smin, r__6 = (r__1 = tl[tl_dim1 + 1], dabs(r__1)), r__5 = max(r__5,
00350             r__6), r__6 = (r__2 = tl[(tl_dim1 << 1) + 1], dabs(r__2)), r__5 = 
00351             max(r__5,r__6), r__6 = (r__3 = tl[tl_dim1 + 2], dabs(r__3)), r__5 
00352             = max(r__5,r__6), r__6 = (r__4 = tl[(tl_dim1 << 1) + 2], dabs(
00353             r__4));
00354     smin = dmax(r__5,r__6);
00355 /* Computing MAX */
00356     r__1 = eps * smin;
00357     smin = dmax(r__1,smlnum);
00358     btmp[0] = 0.f;
00359     scopy_(&c__16, btmp, &c__0, t16, &c__1);
00360     t16[0] = tl[tl_dim1 + 1] + sgn * tr[tr_dim1 + 1];
00361     t16[5] = tl[(tl_dim1 << 1) + 2] + sgn * tr[tr_dim1 + 1];
00362     t16[10] = tl[tl_dim1 + 1] + sgn * tr[(tr_dim1 << 1) + 2];
00363     t16[15] = tl[(tl_dim1 << 1) + 2] + sgn * tr[(tr_dim1 << 1) + 2];
00364     if (*ltranl) {
00365         t16[4] = tl[tl_dim1 + 2];
00366         t16[1] = tl[(tl_dim1 << 1) + 1];
00367         t16[14] = tl[tl_dim1 + 2];
00368         t16[11] = tl[(tl_dim1 << 1) + 1];
00369     } else {
00370         t16[4] = tl[(tl_dim1 << 1) + 1];
00371         t16[1] = tl[tl_dim1 + 2];
00372         t16[14] = tl[(tl_dim1 << 1) + 1];
00373         t16[11] = tl[tl_dim1 + 2];
00374     }
00375     if (*ltranr) {
00376         t16[8] = sgn * tr[(tr_dim1 << 1) + 1];
00377         t16[13] = sgn * tr[(tr_dim1 << 1) + 1];
00378         t16[2] = sgn * tr[tr_dim1 + 2];
00379         t16[7] = sgn * tr[tr_dim1 + 2];
00380     } else {
00381         t16[8] = sgn * tr[tr_dim1 + 2];
00382         t16[13] = sgn * tr[tr_dim1 + 2];
00383         t16[2] = sgn * tr[(tr_dim1 << 1) + 1];
00384         t16[7] = sgn * tr[(tr_dim1 << 1) + 1];
00385     }
00386     btmp[0] = b[b_dim1 + 1];
00387     btmp[1] = b[b_dim1 + 2];
00388     btmp[2] = b[(b_dim1 << 1) + 1];
00389     btmp[3] = b[(b_dim1 << 1) + 2];
00390 
00391 /*     Perform elimination */
00392 
00393     for (i__ = 1; i__ <= 3; ++i__) {
00394         xmax = 0.f;
00395         for (ip = i__; ip <= 4; ++ip) {
00396             for (jp = i__; jp <= 4; ++jp) {
00397                 if ((r__1 = t16[ip + (jp << 2) - 5], dabs(r__1)) >= xmax) {
00398                     xmax = (r__1 = t16[ip + (jp << 2) - 5], dabs(r__1));
00399                     ipsv = ip;
00400                     jpsv = jp;
00401                 }
00402 /* L60: */
00403             }
00404 /* L70: */
00405         }
00406         if (ipsv != i__) {
00407             sswap_(&c__4, &t16[ipsv - 1], &c__4, &t16[i__ - 1], &c__4);
00408             temp = btmp[i__ - 1];
00409             btmp[i__ - 1] = btmp[ipsv - 1];
00410             btmp[ipsv - 1] = temp;
00411         }
00412         if (jpsv != i__) {
00413             sswap_(&c__4, &t16[(jpsv << 2) - 4], &c__1, &t16[(i__ << 2) - 4], 
00414                     &c__1);
00415         }
00416         jpiv[i__ - 1] = jpsv;
00417         if ((r__1 = t16[i__ + (i__ << 2) - 5], dabs(r__1)) < smin) {
00418             *info = 1;
00419             t16[i__ + (i__ << 2) - 5] = smin;
00420         }
00421         for (j = i__ + 1; j <= 4; ++j) {
00422             t16[j + (i__ << 2) - 5] /= t16[i__ + (i__ << 2) - 5];
00423             btmp[j - 1] -= t16[j + (i__ << 2) - 5] * btmp[i__ - 1];
00424             for (k = i__ + 1; k <= 4; ++k) {
00425                 t16[j + (k << 2) - 5] -= t16[j + (i__ << 2) - 5] * t16[i__ + (
00426                         k << 2) - 5];
00427 /* L80: */
00428             }
00429 /* L90: */
00430         }
00431 /* L100: */
00432     }
00433     if (dabs(t16[15]) < smin) {
00434         t16[15] = smin;
00435     }
00436     *scale = 1.f;
00437     if (smlnum * 8.f * dabs(btmp[0]) > dabs(t16[0]) || smlnum * 8.f * dabs(
00438             btmp[1]) > dabs(t16[5]) || smlnum * 8.f * dabs(btmp[2]) > dabs(
00439             t16[10]) || smlnum * 8.f * dabs(btmp[3]) > dabs(t16[15])) {
00440 /* Computing MAX */
00441         r__1 = dabs(btmp[0]), r__2 = dabs(btmp[1]), r__1 = max(r__1,r__2), 
00442                 r__2 = dabs(btmp[2]), r__1 = max(r__1,r__2), r__2 = dabs(btmp[
00443                 3]);
00444         *scale = .125f / dmax(r__1,r__2);
00445         btmp[0] *= *scale;
00446         btmp[1] *= *scale;
00447         btmp[2] *= *scale;
00448         btmp[3] *= *scale;
00449     }
00450     for (i__ = 1; i__ <= 4; ++i__) {
00451         k = 5 - i__;
00452         temp = 1.f / t16[k + (k << 2) - 5];
00453         tmp[k - 1] = btmp[k - 1] * temp;
00454         for (j = k + 1; j <= 4; ++j) {
00455             tmp[k - 1] -= temp * t16[k + (j << 2) - 5] * tmp[j - 1];
00456 /* L110: */
00457         }
00458 /* L120: */
00459     }
00460     for (i__ = 1; i__ <= 3; ++i__) {
00461         if (jpiv[4 - i__ - 1] != 4 - i__) {
00462             temp = tmp[4 - i__ - 1];
00463             tmp[4 - i__ - 1] = tmp[jpiv[4 - i__ - 1] - 1];
00464             tmp[jpiv[4 - i__ - 1] - 1] = temp;
00465         }
00466 /* L130: */
00467     }
00468     x[x_dim1 + 1] = tmp[0];
00469     x[x_dim1 + 2] = tmp[1];
00470     x[(x_dim1 << 1) + 1] = tmp[2];
00471     x[(x_dim1 << 1) + 2] = tmp[3];
00472 /* Computing MAX */
00473     r__1 = dabs(tmp[0]) + dabs(tmp[2]), r__2 = dabs(tmp[1]) + dabs(tmp[3]);
00474     *xnorm = dmax(r__1,r__2);
00475     return 0;
00476 
00477 /*     End of SLASY2 */
00478 
00479 } /* slasy2_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:11