dget35.c
Go to the documentation of this file.
00001 /* dget35.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__6 = 6;
00019 static doublereal c_b35 = 1.;
00020 
00021 /* Subroutine */ int dget35_(doublereal *rmax, integer *lmax, integer *ninfo, 
00022         integer *knt)
00023 {
00024     /* Initialized data */
00025 
00026     static integer idim[8] = { 1,2,3,4,3,3,6,4 };
00027     static integer ival[288]    /* was [6][6][8] */ = { 1,0,0,0,0,0,0,0,0,0,0,
00028             0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,2,0,0,0,0,-2,
00029             0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
00030             0,0,5,1,2,0,0,0,-8,-2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
00031             3,4,0,0,0,0,-5,3,0,0,0,0,1,2,1,4,0,0,-3,-9,-1,1,0,0,0,0,0,0,0,0,0,
00032             0,0,0,0,0,1,0,0,0,0,0,2,3,0,0,0,0,5,6,7,0,0,0,0,0,0,0,0,0,0,0,0,0,
00033             0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,3,-4,0,0,0,2,5,2,0,0,0,0,0,0,0,0,0,
00034             0,0,0,0,0,0,0,0,0,0,0,0,1,2,0,0,0,0,-2,0,0,0,0,0,5,6,3,4,0,0,-1,
00035             -9,-5,2,0,0,8,8,8,8,5,6,9,9,9,9,-7,5,1,0,0,0,0,0,1,5,2,0,0,0,2,
00036             -21,5,0,0,0,1,2,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
00037 
00038     /* System generated locals */
00039     integer i__1, i__2, i__3;
00040     doublereal d__1, d__2, d__3;
00041 
00042     /* Builtin functions */
00043     double sqrt(doublereal), sin(doublereal);
00044 
00045     /* Local variables */
00046     doublereal a[36]    /* was [6][6] */, b[36] /* was [6][6] */, c__[36]       
00047             /* was [6][6] */;
00048     integer i__, j, m, n;
00049     doublereal cc[36]   /* was [6][6] */, vm1[3], vm2[3];
00050     integer ima, imb;
00051     doublereal dum[1], eps, res, res1;
00052     integer info;
00053     doublereal cnrm;
00054     integer isgn;
00055     doublereal rmul, tnrm, xnrm, scale;
00056     extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
00057             integer *, doublereal *, doublereal *, integer *, doublereal *, 
00058             integer *, doublereal *, doublereal *, integer *);
00059     char trana[1], tranb[1];
00060     integer imlda1, imlda2, imldb1;
00061     extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
00062     extern doublereal dlamch_(char *), dlange_(char *, integer *, 
00063             integer *, doublereal *, integer *, doublereal *);
00064     integer imloff, itrana, itranb;
00065     doublereal bignum, smlnum;
00066     extern /* Subroutine */ int dtrsyl_(char *, char *, integer *, integer *, 
00067             integer *, doublereal *, integer *, doublereal *, integer *, 
00068             doublereal *, integer *, doublereal *, integer *);
00069 
00070 
00071 /*  -- LAPACK test routine (version 3.1) -- */
00072 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00073 /*     November 2006 */
00074 
00075 /*     .. Scalar Arguments .. */
00076 /*     .. */
00077 
00078 /*  Purpose */
00079 /*  ======= */
00080 
00081 /*  DGET35 tests DTRSYL, a routine for solving the Sylvester matrix */
00082 /*  equation */
00083 
00084 /*     op(A)*X + ISGN*X*op(B) = scale*C, */
00085 
00086 /*  A and B are assumed to be in Schur canonical form, op() represents an */
00087 /*  optional transpose, and ISGN can be -1 or +1.  Scale is an output */
00088 /*  less than or equal to 1, chosen to avoid overflow in X. */
00089 
00090 /*  The test code verifies that the following residual is order 1: */
00091 
00092 /*     norm(op(A)*X + ISGN*X*op(B) - scale*C) / */
00093 /*         (EPS*max(norm(A),norm(B))*norm(X)) */
00094 
00095 /*  Arguments */
00096 /*  ========== */
00097 
00098 /*  RMAX    (output) DOUBLE PRECISION */
00099 /*          Value of the largest test ratio. */
00100 
00101 /*  LMAX    (output) INTEGER */
00102 /*          Example number where largest test ratio achieved. */
00103 
00104 /*  NINFO   (output) INTEGER */
00105 /*          Number of examples where INFO is nonzero. */
00106 
00107 /*  KNT     (output) INTEGER */
00108 /*          Total number of examples tested. */
00109 
00110 /*  ===================================================================== */
00111 
00112 /*     .. Parameters .. */
00113 /*     .. */
00114 /*     .. Local Scalars .. */
00115 /*     .. */
00116 /*     .. Local Arrays .. */
00117 /*     .. */
00118 /*     .. External Functions .. */
00119 /*     .. */
00120 /*     .. External Subroutines .. */
00121 /*     .. */
00122 /*     .. Intrinsic Functions .. */
00123 /*     .. */
00124 /*     .. Data statements .. */
00125 /*     .. */
00126 /*     .. Executable Statements .. */
00127 
00128 /*     Get machine parameters */
00129 
00130     eps = dlamch_("P");
00131     smlnum = dlamch_("S") * 4. / eps;
00132     bignum = 1. / smlnum;
00133     dlabad_(&smlnum, &bignum);
00134 
00135 /*     Set up test case parameters */
00136 
00137     vm1[0] = sqrt(smlnum);
00138     vm1[1] = 1.;
00139     vm1[2] = sqrt(bignum);
00140     vm2[0] = 1.;
00141     vm2[1] = eps * 2. + 1.;
00142     vm2[2] = 2.;
00143 
00144     *knt = 0;
00145     *ninfo = 0;
00146     *lmax = 0;
00147     *rmax = 0.;
00148 
00149 /*     Begin test loop */
00150 
00151     for (itrana = 1; itrana <= 2; ++itrana) {
00152         for (itranb = 1; itranb <= 2; ++itranb) {
00153             for (isgn = -1; isgn <= 1; isgn += 2) {
00154                 for (ima = 1; ima <= 8; ++ima) {
00155                     for (imlda1 = 1; imlda1 <= 3; ++imlda1) {
00156                         for (imlda2 = 1; imlda2 <= 3; ++imlda2) {
00157                             for (imloff = 1; imloff <= 2; ++imloff) {
00158                                 for (imb = 1; imb <= 8; ++imb) {
00159                                     for (imldb1 = 1; imldb1 <= 3; ++imldb1) {
00160                                         if (itrana == 1) {
00161                                             *(unsigned char *)trana = 'N';
00162                                         }
00163                                         if (itrana == 2) {
00164                                             *(unsigned char *)trana = 'T';
00165                                         }
00166                                         if (itranb == 1) {
00167                                             *(unsigned char *)tranb = 'N';
00168                                         }
00169                                         if (itranb == 2) {
00170                                             *(unsigned char *)tranb = 'T';
00171                                         }
00172                                         m = idim[ima - 1];
00173                                         n = idim[imb - 1];
00174                                         tnrm = 0.;
00175                                         i__1 = m;
00176                                         for (i__ = 1; i__ <= i__1; ++i__) {
00177                                             i__2 = m;
00178                                             for (j = 1; j <= i__2; ++j) {
00179                           a[i__ + j * 6 - 7] = (doublereal) ival[i__ + (j + 
00180                                   ima * 6) * 6 - 43];
00181                           if ((i__3 = i__ - j, abs(i__3)) <= 1) {
00182                               a[i__ + j * 6 - 7] *= vm1[imlda1 - 1];
00183                               a[i__ + j * 6 - 7] *= vm2[imlda2 - 1];
00184                           } else {
00185                               a[i__ + j * 6 - 7] *= vm1[imloff - 1];
00186                           }
00187 /* Computing MAX */
00188                           d__2 = tnrm, d__3 = (d__1 = a[i__ + j * 6 - 7], abs(
00189                                   d__1));
00190                           tnrm = max(d__2,d__3);
00191 /* L10: */
00192                                             }
00193 /* L20: */
00194                                         }
00195                                         i__1 = n;
00196                                         for (i__ = 1; i__ <= i__1; ++i__) {
00197                                             i__2 = n;
00198                                             for (j = 1; j <= i__2; ++j) {
00199                           b[i__ + j * 6 - 7] = (doublereal) ival[i__ + (j + 
00200                                   imb * 6) * 6 - 43];
00201                           if ((i__3 = i__ - j, abs(i__3)) <= 1) {
00202                               b[i__ + j * 6 - 7] *= vm1[imldb1 - 1];
00203                           } else {
00204                               b[i__ + j * 6 - 7] *= vm1[imloff - 1];
00205                           }
00206 /* Computing MAX */
00207                           d__2 = tnrm, d__3 = (d__1 = b[i__ + j * 6 - 7], abs(
00208                                   d__1));
00209                           tnrm = max(d__2,d__3);
00210 /* L30: */
00211                                             }
00212 /* L40: */
00213                                         }
00214                                         cnrm = 0.;
00215                                         i__1 = m;
00216                                         for (i__ = 1; i__ <= i__1; ++i__) {
00217                                             i__2 = n;
00218                                             for (j = 1; j <= i__2; ++j) {
00219                           c__[i__ + j * 6 - 7] = sin((doublereal) (i__ * j));
00220 /* Computing MAX */
00221                           d__1 = cnrm, d__2 = c__[i__ + j * 6 - 7];
00222                           cnrm = max(d__1,d__2);
00223                           cc[i__ + j * 6 - 7] = c__[i__ + j * 6 - 7];
00224 /* L50: */
00225                                             }
00226 /* L60: */
00227                                         }
00228                                         ++(*knt);
00229                                         dtrsyl_(trana, tranb, &isgn, &m, &n, 
00230                                                 a, &c__6, b, &c__6, c__, &
00231                                                 c__6, &scale, &info);
00232                                         if (info != 0) {
00233                                             ++(*ninfo);
00234                                         }
00235                                         xnrm = dlange_("M", &m, &n, c__, &
00236                                                 c__6, dum);
00237                                         rmul = 1.;
00238                                         if (xnrm > 1. && tnrm > 1.) {
00239                                             if (xnrm > bignum / tnrm) {
00240                           rmul = 1. / max(xnrm,tnrm);
00241                                             }
00242                                         }
00243                                         d__1 = -scale * rmul;
00244                                         dgemm_(trana, "N", &m, &n, &m, &rmul, 
00245                                                 a, &c__6, c__, &c__6, &d__1, 
00246                                                 cc, &c__6);
00247                                         d__1 = (doublereal) isgn * rmul;
00248                                         dgemm_("N", tranb, &m, &n, &n, &d__1, 
00249                                                 c__, &c__6, b, &c__6, &c_b35, 
00250                                                 cc, &c__6);
00251                                         res1 = dlange_("M", &m, &n, cc, &c__6, 
00252                                                  dum);
00253 /* Computing MAX */
00254                                         d__1 = smlnum, d__2 = smlnum * xnrm, 
00255                                                 d__1 = max(d__1,d__2), d__2 = 
00256                                                 rmul * tnrm * eps * xnrm;
00257                                         res = res1 / max(d__1,d__2);
00258                                         if (res > *rmax) {
00259                                             *lmax = *knt;
00260                                             *rmax = res;
00261                                         }
00262 /* L70: */
00263                                     }
00264 /* L80: */
00265                                 }
00266 /* L90: */
00267                             }
00268 /* L100: */
00269                         }
00270 /* L110: */
00271                     }
00272 /* L120: */
00273                 }
00274 /* L130: */
00275             }
00276 /* L140: */
00277         }
00278 /* L150: */
00279     }
00280 
00281     return 0;
00282 
00283 /*     End of DGET35 */
00284 
00285 } /* dget35_ */


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