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 stbt05_(char *uplo, char *trans, char *diag, integer *n,
00021 integer *kd, integer *nrhs, real *ab, integer *ldab, real *b, integer
00022 *ldb, real *x, integer *ldx, real *xact, integer *ldxact, real *ferr,
00023 real *berr, real *reslts)
00024 {
00025
00026 integer ab_dim1, ab_offset, b_dim1, b_offset, x_dim1, x_offset, xact_dim1,
00027 xact_offset, i__1, i__2, i__3, i__4;
00028 real r__1, r__2, r__3;
00029
00030
00031 integer i__, j, k, nz, ifu;
00032 real eps, tmp, diff, axbi;
00033 integer imax;
00034 real unfl, ovfl;
00035 logical unit;
00036 extern logical lsame_(char *, char *);
00037 logical upper;
00038 real xnorm;
00039 extern doublereal slamch_(char *);
00040 real errbnd;
00041 extern integer isamax_(integer *, real *, integer *);
00042 logical notran;
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
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166 ab_dim1 = *ldab;
00167 ab_offset = 1 + ab_dim1;
00168 ab -= ab_offset;
00169 b_dim1 = *ldb;
00170 b_offset = 1 + b_dim1;
00171 b -= b_offset;
00172 x_dim1 = *ldx;
00173 x_offset = 1 + x_dim1;
00174 x -= x_offset;
00175 xact_dim1 = *ldxact;
00176 xact_offset = 1 + xact_dim1;
00177 xact -= xact_offset;
00178 --ferr;
00179 --berr;
00180 --reslts;
00181
00182
00183 if (*n <= 0 || *nrhs <= 0) {
00184 reslts[1] = 0.f;
00185 reslts[2] = 0.f;
00186 return 0;
00187 }
00188
00189 eps = slamch_("Epsilon");
00190 unfl = slamch_("Safe minimum");
00191 ovfl = 1.f / unfl;
00192 upper = lsame_(uplo, "U");
00193 notran = lsame_(trans, "N");
00194 unit = lsame_(diag, "U");
00195
00196 i__1 = *kd, i__2 = *n - 1;
00197 nz = min(i__1,i__2) + 1;
00198
00199
00200
00201
00202
00203 errbnd = 0.f;
00204 i__1 = *nrhs;
00205 for (j = 1; j <= i__1; ++j) {
00206 imax = isamax_(n, &x[j * x_dim1 + 1], &c__1);
00207
00208 r__2 = (r__1 = x[imax + j * x_dim1], dabs(r__1));
00209 xnorm = dmax(r__2,unfl);
00210 diff = 0.f;
00211 i__2 = *n;
00212 for (i__ = 1; i__ <= i__2; ++i__) {
00213
00214 r__2 = diff, r__3 = (r__1 = x[i__ + j * x_dim1] - xact[i__ + j *
00215 xact_dim1], dabs(r__1));
00216 diff = dmax(r__2,r__3);
00217
00218 }
00219
00220 if (xnorm > 1.f) {
00221 goto L20;
00222 } else if (diff <= ovfl * xnorm) {
00223 goto L20;
00224 } else {
00225 errbnd = 1.f / eps;
00226 goto L30;
00227 }
00228
00229 L20:
00230 if (diff / xnorm <= ferr[j]) {
00231
00232 r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
00233 errbnd = dmax(r__1,r__2);
00234 } else {
00235 errbnd = 1.f / eps;
00236 }
00237 L30:
00238 ;
00239 }
00240 reslts[1] = errbnd;
00241
00242
00243
00244
00245 ifu = 0;
00246 if (unit) {
00247 ifu = 1;
00248 }
00249 i__1 = *nrhs;
00250 for (k = 1; k <= i__1; ++k) {
00251 i__2 = *n;
00252 for (i__ = 1; i__ <= i__2; ++i__) {
00253 tmp = (r__1 = b[i__ + k * b_dim1], dabs(r__1));
00254 if (upper) {
00255 if (! notran) {
00256
00257 i__3 = i__ - *kd;
00258 i__4 = i__ - ifu;
00259 for (j = max(i__3,1); j <= i__4; ++j) {
00260 tmp += (r__1 = ab[*kd + 1 - i__ + j + i__ * ab_dim1],
00261 dabs(r__1)) * (r__2 = x[j + k * x_dim1], dabs(
00262 r__2));
00263
00264 }
00265 if (unit) {
00266 tmp += (r__1 = x[i__ + k * x_dim1], dabs(r__1));
00267 }
00268 } else {
00269 if (unit) {
00270 tmp += (r__1 = x[i__ + k * x_dim1], dabs(r__1));
00271 }
00272
00273 i__3 = i__ + *kd;
00274 i__4 = min(i__3,*n);
00275 for (j = i__ + ifu; j <= i__4; ++j) {
00276 tmp += (r__1 = ab[*kd + 1 + i__ - j + j * ab_dim1],
00277 dabs(r__1)) * (r__2 = x[j + k * x_dim1], dabs(
00278 r__2));
00279
00280 }
00281 }
00282 } else {
00283 if (notran) {
00284
00285 i__4 = i__ - *kd;
00286 i__3 = i__ - ifu;
00287 for (j = max(i__4,1); j <= i__3; ++j) {
00288 tmp += (r__1 = ab[i__ + 1 - j + j * ab_dim1], dabs(
00289 r__1)) * (r__2 = x[j + k * x_dim1], dabs(r__2)
00290 );
00291
00292 }
00293 if (unit) {
00294 tmp += (r__1 = x[i__ + k * x_dim1], dabs(r__1));
00295 }
00296 } else {
00297 if (unit) {
00298 tmp += (r__1 = x[i__ + k * x_dim1], dabs(r__1));
00299 }
00300
00301 i__4 = i__ + *kd;
00302 i__3 = min(i__4,*n);
00303 for (j = i__ + ifu; j <= i__3; ++j) {
00304 tmp += (r__1 = ab[j + 1 - i__ + i__ * ab_dim1], dabs(
00305 r__1)) * (r__2 = x[j + k * x_dim1], dabs(r__2)
00306 );
00307
00308 }
00309 }
00310 }
00311 if (i__ == 1) {
00312 axbi = tmp;
00313 } else {
00314 axbi = dmin(axbi,tmp);
00315 }
00316
00317 }
00318
00319 r__1 = axbi, r__2 = nz * unfl;
00320 tmp = berr[k] / (nz * eps + nz * unfl / dmax(r__1,r__2));
00321 if (k == 1) {
00322 reslts[2] = tmp;
00323 } else {
00324 reslts[2] = dmax(reslts[2],tmp);
00325 }
00326
00327 }
00328
00329 return 0;
00330
00331
00332
00333 }