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__1 = 1;
00019
00020 int dgtt05_(char *trans, integer *n, integer *nrhs,
00021 doublereal *dl, doublereal *d__, doublereal *du, doublereal *b,
00022 integer *ldb, doublereal *x, integer *ldx, doublereal *xact, integer *
00023 ldxact, doublereal *ferr, doublereal *berr, doublereal *reslts)
00024 {
00025
00026 integer b_dim1, b_offset, x_dim1, x_offset, xact_dim1, xact_offset, i__1,
00027 i__2;
00028 doublereal d__1, d__2, d__3, d__4;
00029
00030
00031 integer i__, j, k, nz;
00032 doublereal eps, tmp, diff, axbi;
00033 integer imax;
00034 doublereal unfl, ovfl;
00035 extern logical lsame_(char *, char *);
00036 doublereal xnorm;
00037 extern doublereal dlamch_(char *);
00038 extern integer idamax_(integer *, doublereal *, integer *);
00039 doublereal errbnd;
00040 logical notran;
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
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
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 --dl;
00147 --d__;
00148 --du;
00149 b_dim1 = *ldb;
00150 b_offset = 1 + b_dim1;
00151 b -= b_offset;
00152 x_dim1 = *ldx;
00153 x_offset = 1 + x_dim1;
00154 x -= x_offset;
00155 xact_dim1 = *ldxact;
00156 xact_offset = 1 + xact_dim1;
00157 xact -= xact_offset;
00158 --ferr;
00159 --berr;
00160 --reslts;
00161
00162
00163 if (*n <= 0 || *nrhs <= 0) {
00164 reslts[1] = 0.;
00165 reslts[2] = 0.;
00166 return 0;
00167 }
00168
00169 eps = dlamch_("Epsilon");
00170 unfl = dlamch_("Safe minimum");
00171 ovfl = 1. / unfl;
00172 notran = lsame_(trans, "N");
00173 nz = 4;
00174
00175
00176
00177
00178
00179 errbnd = 0.;
00180 i__1 = *nrhs;
00181 for (j = 1; j <= i__1; ++j) {
00182 imax = idamax_(n, &x[j * x_dim1 + 1], &c__1);
00183
00184 d__2 = (d__1 = x[imax + j * x_dim1], abs(d__1));
00185 xnorm = max(d__2,unfl);
00186 diff = 0.;
00187 i__2 = *n;
00188 for (i__ = 1; i__ <= i__2; ++i__) {
00189
00190 d__2 = diff, d__3 = (d__1 = x[i__ + j * x_dim1] - xact[i__ + j *
00191 xact_dim1], abs(d__1));
00192 diff = max(d__2,d__3);
00193
00194 }
00195
00196 if (xnorm > 1.) {
00197 goto L20;
00198 } else if (diff <= ovfl * xnorm) {
00199 goto L20;
00200 } else {
00201 errbnd = 1. / eps;
00202 goto L30;
00203 }
00204
00205 L20:
00206 if (diff / xnorm <= ferr[j]) {
00207
00208 d__1 = errbnd, d__2 = diff / xnorm / ferr[j];
00209 errbnd = max(d__1,d__2);
00210 } else {
00211 errbnd = 1. / eps;
00212 }
00213 L30:
00214 ;
00215 }
00216 reslts[1] = errbnd;
00217
00218
00219
00220
00221 i__1 = *nrhs;
00222 for (k = 1; k <= i__1; ++k) {
00223 if (notran) {
00224 if (*n == 1) {
00225 axbi = (d__1 = b[k * b_dim1 + 1], abs(d__1)) + (d__2 = d__[1]
00226 * x[k * x_dim1 + 1], abs(d__2));
00227 } else {
00228 axbi = (d__1 = b[k * b_dim1 + 1], abs(d__1)) + (d__2 = d__[1]
00229 * x[k * x_dim1 + 1], abs(d__2)) + (d__3 = du[1] * x[k
00230 * x_dim1 + 2], abs(d__3));
00231 i__2 = *n - 1;
00232 for (i__ = 2; i__ <= i__2; ++i__) {
00233 tmp = (d__1 = b[i__ + k * b_dim1], abs(d__1)) + (d__2 =
00234 dl[i__ - 1] * x[i__ - 1 + k * x_dim1], abs(d__2))
00235 + (d__3 = d__[i__] * x[i__ + k * x_dim1], abs(
00236 d__3)) + (d__4 = du[i__] * x[i__ + 1 + k * x_dim1]
00237 , abs(d__4));
00238 axbi = min(axbi,tmp);
00239
00240 }
00241 tmp = (d__1 = b[*n + k * b_dim1], abs(d__1)) + (d__2 = dl[*n
00242 - 1] * x[*n - 1 + k * x_dim1], abs(d__2)) + (d__3 =
00243 d__[*n] * x[*n + k * x_dim1], abs(d__3));
00244 axbi = min(axbi,tmp);
00245 }
00246 } else {
00247 if (*n == 1) {
00248 axbi = (d__1 = b[k * b_dim1 + 1], abs(d__1)) + (d__2 = d__[1]
00249 * x[k * x_dim1 + 1], abs(d__2));
00250 } else {
00251 axbi = (d__1 = b[k * b_dim1 + 1], abs(d__1)) + (d__2 = d__[1]
00252 * x[k * x_dim1 + 1], abs(d__2)) + (d__3 = dl[1] * x[k
00253 * x_dim1 + 2], abs(d__3));
00254 i__2 = *n - 1;
00255 for (i__ = 2; i__ <= i__2; ++i__) {
00256 tmp = (d__1 = b[i__ + k * b_dim1], abs(d__1)) + (d__2 =
00257 du[i__ - 1] * x[i__ - 1 + k * x_dim1], abs(d__2))
00258 + (d__3 = d__[i__] * x[i__ + k * x_dim1], abs(
00259 d__3)) + (d__4 = dl[i__] * x[i__ + 1 + k * x_dim1]
00260 , abs(d__4));
00261 axbi = min(axbi,tmp);
00262
00263 }
00264 tmp = (d__1 = b[*n + k * b_dim1], abs(d__1)) + (d__2 = du[*n
00265 - 1] * x[*n - 1 + k * x_dim1], abs(d__2)) + (d__3 =
00266 d__[*n] * x[*n + k * x_dim1], abs(d__3));
00267 axbi = min(axbi,tmp);
00268 }
00269 }
00270
00271 d__1 = axbi, d__2 = nz * unfl;
00272 tmp = berr[k] / (nz * eps + nz * unfl / max(d__1,d__2));
00273 if (k == 1) {
00274 reslts[2] = tmp;
00275 } else {
00276 reslts[2] = max(reslts[2],tmp);
00277 }
00278
00279 }
00280
00281 return 0;
00282
00283
00284
00285 }