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__3 = 3;
00019 static integer c__1 = 1;
00020 static integer c__5 = 5;
00021 static integer c__50 = 50;
00022 static doublereal c_b52 = 1.;
00023 static doublereal c_b55 = 0.;
00024
00025 int dchkgk_(integer *nin, integer *nout)
00026 {
00027
00028 static char fmt_9999[] = "(1x,\002.. test output of DGGBAK .. \002)";
00029 static char fmt_9998[] = "(\002 value of largest test error "
00030 " =\002,d12.3)";
00031 static char fmt_9997[] = "(\002 example number where DGGBAL info is not "
00032 "0 =\002,i4)";
00033 static char fmt_9996[] = "(\002 example number where DGGBAK(L) info is n"
00034 "ot 0 =\002,i4)";
00035 static char fmt_9995[] = "(\002 example number where DGGBAK(R) info is n"
00036 "ot 0 =\002,i4)";
00037 static char fmt_9994[] = "(\002 example number having largest error "
00038 " =\002,i4)";
00039 static char fmt_9993[] = "(\002 number of examples where info is not 0 "
00040 " =\002,i4)";
00041 static char fmt_9992[] = "(\002 total number of examples tested "
00042 " =\002,i4)";
00043
00044
00045 integer i__1, i__2;
00046 doublereal d__1, d__2, d__3;
00047
00048
00049 integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
00050 e_rsle(void), s_wsfe(cilist *), e_wsfe(void), do_fio(integer *,
00051 char *, ftnlen);
00052
00053
00054 doublereal a[2500] , b[2500] , e[
00055 2500] , f[2500] ;
00056 integer i__, j, m, n;
00057 doublereal af[2500] , bf[2500] ,
00058 vl[2500] , vr[2500] ;
00059 integer ihi, ilo;
00060 doublereal eps, vlf[2500] ;
00061 integer knt;
00062 doublereal vrf[2500] ;
00063 integer info, lmax[4];
00064 doublereal rmax, vmax, work[2500] ;
00065 extern int dgemm_(char *, char *, integer *, integer *,
00066 integer *, doublereal *, doublereal *, integer *, doublereal *,
00067 integer *, doublereal *, doublereal *, integer *);
00068 integer ninfo;
00069 doublereal anorm, bnorm;
00070 extern int dggbak_(char *, char *, integer *, integer *,
00071 integer *, doublereal *, doublereal *, integer *, doublereal *,
00072 integer *, integer *), dggbal_(char *, integer *,
00073 doublereal *, integer *, doublereal *, integer *, integer *,
00074 integer *, doublereal *, doublereal *, doublereal *, integer *);
00075 extern doublereal dlamch_(char *), dlange_(char *, integer *,
00076 integer *, doublereal *, integer *, doublereal *);
00077 doublereal lscale[50], rscale[50];
00078 extern int dlacpy_(char *, integer *, integer *,
00079 doublereal *, integer *, doublereal *, integer *);
00080
00081
00082 static cilist io___6 = { 0, 0, 0, 0, 0 };
00083 static cilist io___10 = { 0, 0, 0, 0, 0 };
00084 static cilist io___13 = { 0, 0, 0, 0, 0 };
00085 static cilist io___15 = { 0, 0, 0, 0, 0 };
00086 static cilist io___17 = { 0, 0, 0, 0, 0 };
00087 static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
00088 static cilist io___35 = { 0, 0, 0, fmt_9998, 0 };
00089 static cilist io___36 = { 0, 0, 0, fmt_9997, 0 };
00090 static cilist io___37 = { 0, 0, 0, fmt_9996, 0 };
00091 static cilist io___38 = { 0, 0, 0, fmt_9995, 0 };
00092 static cilist io___39 = { 0, 0, 0, fmt_9994, 0 };
00093 static cilist io___40 = { 0, 0, 0, fmt_9993, 0 };
00094 static cilist io___41 = { 0, 0, 0, fmt_9992, 0 };
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
00131
00132
00133
00134
00135
00136
00137
00138 lmax[0] = 0;
00139 lmax[1] = 0;
00140 lmax[2] = 0;
00141 lmax[3] = 0;
00142 ninfo = 0;
00143 knt = 0;
00144 rmax = 0.;
00145
00146 eps = dlamch_("Precision");
00147
00148 L10:
00149 io___6.ciunit = *nin;
00150 s_rsle(&io___6);
00151 do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
00152 do_lio(&c__3, &c__1, (char *)&m, (ftnlen)sizeof(integer));
00153 e_rsle();
00154 if (n == 0) {
00155 goto L100;
00156 }
00157
00158 i__1 = n;
00159 for (i__ = 1; i__ <= i__1; ++i__) {
00160 io___10.ciunit = *nin;
00161 s_rsle(&io___10);
00162 i__2 = n;
00163 for (j = 1; j <= i__2; ++j) {
00164 do_lio(&c__5, &c__1, (char *)&a[i__ + j * 50 - 51], (ftnlen)
00165 sizeof(doublereal));
00166 }
00167 e_rsle();
00168
00169 }
00170
00171 i__1 = n;
00172 for (i__ = 1; i__ <= i__1; ++i__) {
00173 io___13.ciunit = *nin;
00174 s_rsle(&io___13);
00175 i__2 = n;
00176 for (j = 1; j <= i__2; ++j) {
00177 do_lio(&c__5, &c__1, (char *)&b[i__ + j * 50 - 51], (ftnlen)
00178 sizeof(doublereal));
00179 }
00180 e_rsle();
00181
00182 }
00183
00184 i__1 = n;
00185 for (i__ = 1; i__ <= i__1; ++i__) {
00186 io___15.ciunit = *nin;
00187 s_rsle(&io___15);
00188 i__2 = m;
00189 for (j = 1; j <= i__2; ++j) {
00190 do_lio(&c__5, &c__1, (char *)&vl[i__ + j * 50 - 51], (ftnlen)
00191 sizeof(doublereal));
00192 }
00193 e_rsle();
00194
00195 }
00196
00197 i__1 = n;
00198 for (i__ = 1; i__ <= i__1; ++i__) {
00199 io___17.ciunit = *nin;
00200 s_rsle(&io___17);
00201 i__2 = m;
00202 for (j = 1; j <= i__2; ++j) {
00203 do_lio(&c__5, &c__1, (char *)&vr[i__ + j * 50 - 51], (ftnlen)
00204 sizeof(doublereal));
00205 }
00206 e_rsle();
00207
00208 }
00209
00210 ++knt;
00211
00212 anorm = dlange_("M", &n, &n, a, &c__50, work);
00213 bnorm = dlange_("M", &n, &n, b, &c__50, work);
00214
00215 dlacpy_("FULL", &n, &n, a, &c__50, af, &c__50);
00216 dlacpy_("FULL", &n, &n, b, &c__50, bf, &c__50);
00217
00218 dggbal_("B", &n, a, &c__50, b, &c__50, &ilo, &ihi, lscale, rscale, work, &
00219 info);
00220 if (info != 0) {
00221 ++ninfo;
00222 lmax[0] = knt;
00223 }
00224
00225 dlacpy_("FULL", &n, &m, vl, &c__50, vlf, &c__50);
00226 dlacpy_("FULL", &n, &m, vr, &c__50, vrf, &c__50);
00227
00228 dggbak_("B", "L", &n, &ilo, &ihi, lscale, rscale, &m, vl, &c__50, &info);
00229 if (info != 0) {
00230 ++ninfo;
00231 lmax[1] = knt;
00232 }
00233
00234 dggbak_("B", "R", &n, &ilo, &ihi, lscale, rscale, &m, vr, &c__50, &info);
00235 if (info != 0) {
00236 ++ninfo;
00237 lmax[2] = knt;
00238 }
00239
00240
00241
00242
00243
00244
00245 dgemm_("N", "N", &n, &m, &n, &c_b52, af, &c__50, vr, &c__50, &c_b55, work,
00246 &c__50);
00247 dgemm_("T", "N", &m, &m, &n, &c_b52, vl, &c__50, work, &c__50, &c_b55, e,
00248 &c__50);
00249
00250 dgemm_("N", "N", &n, &m, &n, &c_b52, a, &c__50, vrf, &c__50, &c_b55, work,
00251 &c__50);
00252 dgemm_("T", "N", &m, &m, &n, &c_b52, vlf, &c__50, work, &c__50, &c_b55, f,
00253 &c__50);
00254
00255 vmax = 0.;
00256 i__1 = m;
00257 for (j = 1; j <= i__1; ++j) {
00258 i__2 = m;
00259 for (i__ = 1; i__ <= i__2; ++i__) {
00260
00261 d__2 = vmax, d__3 = (d__1 = e[i__ + j * 50 - 51] - f[i__ + j * 50
00262 - 51], abs(d__1));
00263 vmax = max(d__2,d__3);
00264
00265 }
00266
00267 }
00268 vmax /= eps * max(anorm,bnorm);
00269 if (vmax > rmax) {
00270 lmax[3] = knt;
00271 rmax = vmax;
00272 }
00273
00274
00275
00276 dgemm_("N", "N", &n, &m, &n, &c_b52, bf, &c__50, vr, &c__50, &c_b55, work,
00277 &c__50);
00278 dgemm_("T", "N", &m, &m, &n, &c_b52, vl, &c__50, work, &c__50, &c_b55, e,
00279 &c__50);
00280
00281 dgemm_("N", "N", &n, &m, &n, &c_b52, b, &c__50, vrf, &c__50, &c_b55, work,
00282 &c__50);
00283 dgemm_("T", "N", &m, &m, &n, &c_b52, vlf, &c__50, work, &c__50, &c_b55, f,
00284 &c__50);
00285
00286 vmax = 0.;
00287 i__1 = m;
00288 for (j = 1; j <= i__1; ++j) {
00289 i__2 = m;
00290 for (i__ = 1; i__ <= i__2; ++i__) {
00291
00292 d__2 = vmax, d__3 = (d__1 = e[i__ + j * 50 - 51] - f[i__ + j * 50
00293 - 51], abs(d__1));
00294 vmax = max(d__2,d__3);
00295
00296 }
00297
00298 }
00299 vmax /= eps * max(anorm,bnorm);
00300 if (vmax > rmax) {
00301 lmax[3] = knt;
00302 rmax = vmax;
00303 }
00304
00305 goto L10;
00306
00307 L100:
00308
00309 io___34.ciunit = *nout;
00310 s_wsfe(&io___34);
00311 e_wsfe();
00312
00313 io___35.ciunit = *nout;
00314 s_wsfe(&io___35);
00315 do_fio(&c__1, (char *)&rmax, (ftnlen)sizeof(doublereal));
00316 e_wsfe();
00317 io___36.ciunit = *nout;
00318 s_wsfe(&io___36);
00319 do_fio(&c__1, (char *)&lmax[0], (ftnlen)sizeof(integer));
00320 e_wsfe();
00321 io___37.ciunit = *nout;
00322 s_wsfe(&io___37);
00323 do_fio(&c__1, (char *)&lmax[1], (ftnlen)sizeof(integer));
00324 e_wsfe();
00325 io___38.ciunit = *nout;
00326 s_wsfe(&io___38);
00327 do_fio(&c__1, (char *)&lmax[2], (ftnlen)sizeof(integer));
00328 e_wsfe();
00329 io___39.ciunit = *nout;
00330 s_wsfe(&io___39);
00331 do_fio(&c__1, (char *)&lmax[3], (ftnlen)sizeof(integer));
00332 e_wsfe();
00333 io___40.ciunit = *nout;
00334 s_wsfe(&io___40);
00335 do_fio(&c__1, (char *)&ninfo, (ftnlen)sizeof(integer));
00336 e_wsfe();
00337 io___41.ciunit = *nout;
00338 s_wsfe(&io___41);
00339 do_fio(&c__1, (char *)&knt, (ftnlen)sizeof(integer));
00340 e_wsfe();
00341
00342 return 0;
00343
00344
00345
00346 }