00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 static integer c__6 = 6;
00019 static doublereal c_b35 = 1.;
00020
00021 int dget35_(doublereal *rmax, integer *lmax, integer *ninfo,
00022 integer *knt)
00023 {
00024
00025
00026 static integer idim[8] = { 1,2,3,4,3,3,6,4 };
00027 static integer ival[288] = { 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
00039 integer i__1, i__2, i__3;
00040 doublereal d__1, d__2, d__3;
00041
00042
00043 double sqrt(doublereal), sin(doublereal);
00044
00045
00046 doublereal a[36] , b[36] , c__[36]
00047 ;
00048 integer i__, j, m, n;
00049 doublereal cc[36] , 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 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 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 int dtrsyl_(char *, char *, integer *, integer *,
00067 integer *, doublereal *, integer *, doublereal *, integer *,
00068 doublereal *, integer *, doublereal *, integer *);
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130 eps = dlamch_("P");
00131 smlnum = dlamch_("S") * 4. / eps;
00132 bignum = 1. / smlnum;
00133 dlabad_(&smlnum, &bignum);
00134
00135
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
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
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
00192 }
00193
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
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
00211 }
00212
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
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
00225 }
00226
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
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
00263 }
00264
00265 }
00266
00267 }
00268
00269 }
00270
00271 }
00272
00273 }
00274
00275 }
00276
00277 }
00278
00279 }
00280
00281 return 0;
00282
00283
00284
00285 }