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