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__20 = 20;
00022
00023 int dchkbl_(integer *nin, integer *nout)
00024 {
00025
00026 static char fmt_9999[] = "(1x,\002.. test output of DGEBAL .. \002)";
00027 static char fmt_9998[] = "(1x,\002value of largest test error "
00028 " = \002,d12.3)";
00029 static char fmt_9997[] = "(1x,\002example number where info is not zero "
00030 " = \002,i4)";
00031 static char fmt_9996[] = "(1x,\002example number where ILO or IHI wrong "
00032 " = \002,i4)";
00033 static char fmt_9995[] = "(1x,\002example number having largest error "
00034 " = \002,i4)";
00035 static char fmt_9994[] = "(1x,\002number of examples where info is not 0"
00036 " = \002,i4)";
00037 static char fmt_9993[] = "(1x,\002total number of examples tested "
00038 " = \002,i4)";
00039
00040
00041 integer i__1, i__2;
00042 doublereal d__1, d__2, d__3;
00043
00044
00045 integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
00046 e_rsle(void), s_wsfe(cilist *), e_wsfe(void), do_fio(integer *,
00047 char *, ftnlen);
00048
00049
00050 doublereal a[400] ;
00051 integer i__, j, n;
00052 doublereal ain[400] ;
00053 integer ihi, ilo, knt, info, lmax[3];
00054 doublereal meps, temp, rmax, vmax, scale[20];
00055 integer ihiin, ninfo, iloin;
00056 doublereal anorm, sfmin, dummy[1];
00057 extern int dgebal_(char *, integer *, doublereal *,
00058 integer *, integer *, integer *, doublereal *, integer *);
00059 extern doublereal dlamch_(char *), dlange_(char *, integer *,
00060 integer *, doublereal *, integer *, doublereal *);
00061 doublereal scalin[20];
00062
00063
00064 static cilist io___8 = { 0, 0, 0, 0, 0 };
00065 static cilist io___11 = { 0, 0, 0, 0, 0 };
00066 static cilist io___14 = { 0, 0, 0, 0, 0 };
00067 static cilist io___17 = { 0, 0, 0, 0, 0 };
00068 static cilist io___19 = { 0, 0, 0, 0, 0 };
00069 static cilist io___28 = { 0, 0, 0, fmt_9999, 0 };
00070 static cilist io___29 = { 0, 0, 0, fmt_9998, 0 };
00071 static cilist io___30 = { 0, 0, 0, fmt_9997, 0 };
00072 static cilist io___31 = { 0, 0, 0, fmt_9996, 0 };
00073 static cilist io___32 = { 0, 0, 0, fmt_9995, 0 };
00074 static cilist io___33 = { 0, 0, 0, fmt_9994, 0 };
00075 static cilist io___34 = { 0, 0, 0, fmt_9993, 0 };
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 lmax[0] = 0;
00118 lmax[1] = 0;
00119 lmax[2] = 0;
00120 ninfo = 0;
00121 knt = 0;
00122 rmax = 0.;
00123 vmax = 0.;
00124 sfmin = dlamch_("S");
00125 meps = dlamch_("E");
00126
00127 L10:
00128
00129 io___8.ciunit = *nin;
00130 s_rsle(&io___8);
00131 do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer));
00132 e_rsle();
00133 if (n == 0) {
00134 goto L70;
00135 }
00136 i__1 = n;
00137 for (i__ = 1; i__ <= i__1; ++i__) {
00138 io___11.ciunit = *nin;
00139 s_rsle(&io___11);
00140 i__2 = n;
00141 for (j = 1; j <= i__2; ++j) {
00142 do_lio(&c__5, &c__1, (char *)&a[i__ + j * 20 - 21], (ftnlen)
00143 sizeof(doublereal));
00144 }
00145 e_rsle();
00146
00147 }
00148
00149 io___14.ciunit = *nin;
00150 s_rsle(&io___14);
00151 do_lio(&c__3, &c__1, (char *)&iloin, (ftnlen)sizeof(integer));
00152 do_lio(&c__3, &c__1, (char *)&ihiin, (ftnlen)sizeof(integer));
00153 e_rsle();
00154 i__1 = n;
00155 for (i__ = 1; i__ <= i__1; ++i__) {
00156 io___17.ciunit = *nin;
00157 s_rsle(&io___17);
00158 i__2 = n;
00159 for (j = 1; j <= i__2; ++j) {
00160 do_lio(&c__5, &c__1, (char *)&ain[i__ + j * 20 - 21], (ftnlen)
00161 sizeof(doublereal));
00162 }
00163 e_rsle();
00164
00165 }
00166 io___19.ciunit = *nin;
00167 s_rsle(&io___19);
00168 i__1 = n;
00169 for (i__ = 1; i__ <= i__1; ++i__) {
00170 do_lio(&c__5, &c__1, (char *)&scalin[i__ - 1], (ftnlen)sizeof(
00171 doublereal));
00172 }
00173 e_rsle();
00174
00175 anorm = dlange_("M", &n, &n, a, &c__20, dummy);
00176 ++knt;
00177
00178 dgebal_("B", &n, a, &c__20, &ilo, &ihi, scale, &info);
00179
00180 if (info != 0) {
00181 ++ninfo;
00182 lmax[0] = knt;
00183 }
00184
00185 if (ilo != iloin || ihi != ihiin) {
00186 ++ninfo;
00187 lmax[1] = knt;
00188 }
00189
00190 i__1 = n;
00191 for (i__ = 1; i__ <= i__1; ++i__) {
00192 i__2 = n;
00193 for (j = 1; j <= i__2; ++j) {
00194
00195 d__1 = a[i__ + j * 20 - 21], d__2 = ain[i__ + j * 20 - 21];
00196 temp = max(d__1,d__2);
00197 temp = max(temp,sfmin);
00198
00199 d__2 = vmax, d__3 = (d__1 = a[i__ + j * 20 - 21] - ain[i__ + j *
00200 20 - 21], abs(d__1)) / temp;
00201 vmax = max(d__2,d__3);
00202
00203 }
00204
00205 }
00206
00207 i__1 = n;
00208 for (i__ = 1; i__ <= i__1; ++i__) {
00209
00210 d__1 = scale[i__ - 1], d__2 = scalin[i__ - 1];
00211 temp = max(d__1,d__2);
00212 temp = max(temp,sfmin);
00213
00214 d__2 = vmax, d__3 = (d__1 = scale[i__ - 1] - scalin[i__ - 1], abs(
00215 d__1)) / temp;
00216 vmax = max(d__2,d__3);
00217
00218 }
00219
00220
00221 if (vmax > rmax) {
00222 lmax[2] = knt;
00223 rmax = vmax;
00224 }
00225
00226 goto L10;
00227
00228 L70:
00229
00230 io___28.ciunit = *nout;
00231 s_wsfe(&io___28);
00232 e_wsfe();
00233
00234 io___29.ciunit = *nout;
00235 s_wsfe(&io___29);
00236 do_fio(&c__1, (char *)&rmax, (ftnlen)sizeof(doublereal));
00237 e_wsfe();
00238 io___30.ciunit = *nout;
00239 s_wsfe(&io___30);
00240 do_fio(&c__1, (char *)&lmax[0], (ftnlen)sizeof(integer));
00241 e_wsfe();
00242 io___31.ciunit = *nout;
00243 s_wsfe(&io___31);
00244 do_fio(&c__1, (char *)&lmax[1], (ftnlen)sizeof(integer));
00245 e_wsfe();
00246 io___32.ciunit = *nout;
00247 s_wsfe(&io___32);
00248 do_fio(&c__1, (char *)&lmax[2], (ftnlen)sizeof(integer));
00249 e_wsfe();
00250 io___33.ciunit = *nout;
00251 s_wsfe(&io___33);
00252 do_fio(&c__1, (char *)&ninfo, (ftnlen)sizeof(integer));
00253 e_wsfe();
00254 io___34.ciunit = *nout;
00255 s_wsfe(&io___34);
00256 do_fio(&c__1, (char *)&knt, (ftnlen)sizeof(integer));
00257 e_wsfe();
00258
00259 return 0;
00260
00261
00262
00263 }