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 sgtt05_(char *trans, integer *n, integer *nrhs, real *dl,
00021 real *d__, real *du, real *b, integer *ldb, real *x, integer *ldx,
00022 real *xact, integer *ldxact, real *ferr, real *berr, real *reslts)
00023 {
00024
00025 integer b_dim1, b_offset, x_dim1, x_offset, xact_dim1, xact_offset, i__1,
00026 i__2;
00027 real r__1, r__2, r__3, r__4;
00028
00029
00030 integer i__, j, k, nz;
00031 real eps, tmp, diff, axbi;
00032 integer imax;
00033 real unfl, ovfl;
00034 extern logical lsame_(char *, char *);
00035 real xnorm;
00036 extern doublereal slamch_(char *);
00037 real errbnd;
00038 extern integer isamax_(integer *, real *, integer *);
00039 logical notran;
00040
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 --dl;
00146 --d__;
00147 --du;
00148 b_dim1 = *ldb;
00149 b_offset = 1 + b_dim1;
00150 b -= b_offset;
00151 x_dim1 = *ldx;
00152 x_offset = 1 + x_dim1;
00153 x -= x_offset;
00154 xact_dim1 = *ldxact;
00155 xact_offset = 1 + xact_dim1;
00156 xact -= xact_offset;
00157 --ferr;
00158 --berr;
00159 --reslts;
00160
00161
00162 if (*n <= 0 || *nrhs <= 0) {
00163 reslts[1] = 0.f;
00164 reslts[2] = 0.f;
00165 return 0;
00166 }
00167
00168 eps = slamch_("Epsilon");
00169 unfl = slamch_("Safe minimum");
00170 ovfl = 1.f / unfl;
00171 notran = lsame_(trans, "N");
00172 nz = 4;
00173
00174
00175
00176
00177
00178 errbnd = 0.f;
00179 i__1 = *nrhs;
00180 for (j = 1; j <= i__1; ++j) {
00181 imax = isamax_(n, &x[j * x_dim1 + 1], &c__1);
00182
00183 r__2 = (r__1 = x[imax + j * x_dim1], dabs(r__1));
00184 xnorm = dmax(r__2,unfl);
00185 diff = 0.f;
00186 i__2 = *n;
00187 for (i__ = 1; i__ <= i__2; ++i__) {
00188
00189 r__2 = diff, r__3 = (r__1 = x[i__ + j * x_dim1] - xact[i__ + j *
00190 xact_dim1], dabs(r__1));
00191 diff = dmax(r__2,r__3);
00192
00193 }
00194
00195 if (xnorm > 1.f) {
00196 goto L20;
00197 } else if (diff <= ovfl * xnorm) {
00198 goto L20;
00199 } else {
00200 errbnd = 1.f / eps;
00201 goto L30;
00202 }
00203
00204 L20:
00205 if (diff / xnorm <= ferr[j]) {
00206
00207 r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
00208 errbnd = dmax(r__1,r__2);
00209 } else {
00210 errbnd = 1.f / eps;
00211 }
00212 L30:
00213 ;
00214 }
00215 reslts[1] = errbnd;
00216
00217
00218
00219
00220 i__1 = *nrhs;
00221 for (k = 1; k <= i__1; ++k) {
00222 if (notran) {
00223 if (*n == 1) {
00224 axbi = (r__1 = b[k * b_dim1 + 1], dabs(r__1)) + (r__2 = d__[1]
00225 * x[k * x_dim1 + 1], dabs(r__2));
00226 } else {
00227 axbi = (r__1 = b[k * b_dim1 + 1], dabs(r__1)) + (r__2 = d__[1]
00228 * x[k * x_dim1 + 1], dabs(r__2)) + (r__3 = du[1] * x[
00229 k * x_dim1 + 2], dabs(r__3));
00230 i__2 = *n - 1;
00231 for (i__ = 2; i__ <= i__2; ++i__) {
00232 tmp = (r__1 = b[i__ + k * b_dim1], dabs(r__1)) + (r__2 =
00233 dl[i__ - 1] * x[i__ - 1 + k * x_dim1], dabs(r__2))
00234 + (r__3 = d__[i__] * x[i__ + k * x_dim1], dabs(
00235 r__3)) + (r__4 = du[i__] * x[i__ + 1 + k * x_dim1]
00236 , dabs(r__4));
00237 axbi = dmin(axbi,tmp);
00238
00239 }
00240 tmp = (r__1 = b[*n + k * b_dim1], dabs(r__1)) + (r__2 = dl[*n
00241 - 1] * x[*n - 1 + k * x_dim1], dabs(r__2)) + (r__3 =
00242 d__[*n] * x[*n + k * x_dim1], dabs(r__3));
00243 axbi = dmin(axbi,tmp);
00244 }
00245 } else {
00246 if (*n == 1) {
00247 axbi = (r__1 = b[k * b_dim1 + 1], dabs(r__1)) + (r__2 = d__[1]
00248 * x[k * x_dim1 + 1], dabs(r__2));
00249 } else {
00250 axbi = (r__1 = b[k * b_dim1 + 1], dabs(r__1)) + (r__2 = d__[1]
00251 * x[k * x_dim1 + 1], dabs(r__2)) + (r__3 = dl[1] * x[
00252 k * x_dim1 + 2], dabs(r__3));
00253 i__2 = *n - 1;
00254 for (i__ = 2; i__ <= i__2; ++i__) {
00255 tmp = (r__1 = b[i__ + k * b_dim1], dabs(r__1)) + (r__2 =
00256 du[i__ - 1] * x[i__ - 1 + k * x_dim1], dabs(r__2))
00257 + (r__3 = d__[i__] * x[i__ + k * x_dim1], dabs(
00258 r__3)) + (r__4 = dl[i__] * x[i__ + 1 + k * x_dim1]
00259 , dabs(r__4));
00260 axbi = dmin(axbi,tmp);
00261
00262 }
00263 tmp = (r__1 = b[*n + k * b_dim1], dabs(r__1)) + (r__2 = du[*n
00264 - 1] * x[*n - 1 + k * x_dim1], dabs(r__2)) + (r__3 =
00265 d__[*n] * x[*n + k * x_dim1], dabs(r__3));
00266 axbi = dmin(axbi,tmp);
00267 }
00268 }
00269
00270 r__1 = axbi, r__2 = nz * unfl;
00271 tmp = berr[k] / (nz * eps + nz * unfl / dmax(r__1,r__2));
00272 if (k == 1) {
00273 reslts[2] = tmp;
00274 } else {
00275 reslts[2] = dmax(reslts[2],tmp);
00276 }
00277
00278 }
00279
00280 return 0;
00281
00282
00283
00284 }