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 cgtt05_(char *trans, integer *n, integer *nrhs, complex *
00021 dl, complex *d__, complex *du, complex *b, integer *ldb, complex *x,
00022 integer *ldx, complex *xact, integer *ldxact, real *ferr, real *berr,
00023 real *reslts)
00024 {
00025
00026 integer b_dim1, b_offset, x_dim1, x_offset, xact_dim1, xact_offset, i__1,
00027 i__2, i__3, i__4, i__5, i__6, i__7, i__8, i__9;
00028 real r__1, r__2, r__3, r__4, r__5, r__6, r__7, r__8, r__9, r__10, r__11,
00029 r__12, r__13, r__14;
00030 complex q__1, q__2;
00031
00032
00033 double r_imag(complex *);
00034
00035
00036 integer i__, j, k, nz;
00037 real eps, tmp, diff, axbi;
00038 integer imax;
00039 real unfl, ovfl;
00040 extern logical lsame_(char *, char *);
00041 real xnorm;
00042 extern integer icamax_(integer *, complex *, integer *);
00043 extern doublereal slamch_(char *);
00044 real errbnd;
00045 logical notran;
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 --dl;
00156 --d__;
00157 --du;
00158 b_dim1 = *ldb;
00159 b_offset = 1 + b_dim1;
00160 b -= b_offset;
00161 x_dim1 = *ldx;
00162 x_offset = 1 + x_dim1;
00163 x -= x_offset;
00164 xact_dim1 = *ldxact;
00165 xact_offset = 1 + xact_dim1;
00166 xact -= xact_offset;
00167 --ferr;
00168 --berr;
00169 --reslts;
00170
00171
00172 if (*n <= 0 || *nrhs <= 0) {
00173 reslts[1] = 0.f;
00174 reslts[2] = 0.f;
00175 return 0;
00176 }
00177
00178 eps = slamch_("Epsilon");
00179 unfl = slamch_("Safe minimum");
00180 ovfl = 1.f / unfl;
00181 notran = lsame_(trans, "N");
00182 nz = 4;
00183
00184
00185
00186
00187
00188 errbnd = 0.f;
00189 i__1 = *nrhs;
00190 for (j = 1; j <= i__1; ++j) {
00191 imax = icamax_(n, &x[j * x_dim1 + 1], &c__1);
00192
00193 i__2 = imax + j * x_dim1;
00194 r__3 = (r__1 = x[i__2].r, dabs(r__1)) + (r__2 = r_imag(&x[imax + j *
00195 x_dim1]), dabs(r__2));
00196 xnorm = dmax(r__3,unfl);
00197 diff = 0.f;
00198 i__2 = *n;
00199 for (i__ = 1; i__ <= i__2; ++i__) {
00200 i__3 = i__ + j * x_dim1;
00201 i__4 = i__ + j * xact_dim1;
00202 q__2.r = x[i__3].r - xact[i__4].r, q__2.i = x[i__3].i - xact[i__4]
00203 .i;
00204 q__1.r = q__2.r, q__1.i = q__2.i;
00205
00206 r__3 = diff, r__4 = (r__1 = q__1.r, dabs(r__1)) + (r__2 = r_imag(&
00207 q__1), dabs(r__2));
00208 diff = dmax(r__3,r__4);
00209
00210 }
00211
00212 if (xnorm > 1.f) {
00213 goto L20;
00214 } else if (diff <= ovfl * xnorm) {
00215 goto L20;
00216 } else {
00217 errbnd = 1.f / eps;
00218 goto L30;
00219 }
00220
00221 L20:
00222 if (diff / xnorm <= ferr[j]) {
00223
00224 r__1 = errbnd, r__2 = diff / xnorm / ferr[j];
00225 errbnd = dmax(r__1,r__2);
00226 } else {
00227 errbnd = 1.f / eps;
00228 }
00229 L30:
00230 ;
00231 }
00232 reslts[1] = errbnd;
00233
00234
00235
00236
00237 i__1 = *nrhs;
00238 for (k = 1; k <= i__1; ++k) {
00239 if (notran) {
00240 if (*n == 1) {
00241 i__2 = k * b_dim1 + 1;
00242 i__3 = k * x_dim1 + 1;
00243 axbi = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[k *
00244 b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, dabs(
00245 r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * ((
00246 r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[k *
00247 x_dim1 + 1]), dabs(r__6)));
00248 } else {
00249 i__2 = k * b_dim1 + 1;
00250 i__3 = k * x_dim1 + 1;
00251 i__4 = k * x_dim1 + 2;
00252 axbi = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[k *
00253 b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, dabs(
00254 r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * ((
00255 r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[k *
00256 x_dim1 + 1]), dabs(r__6))) + ((r__7 = du[1].r, dabs(
00257 r__7)) + (r__8 = r_imag(&du[1]), dabs(r__8))) * ((
00258 r__9 = x[i__4].r, dabs(r__9)) + (r__10 = r_imag(&x[k *
00259 x_dim1 + 2]), dabs(r__10)));
00260 i__2 = *n - 1;
00261 for (i__ = 2; i__ <= i__2; ++i__) {
00262 i__3 = i__ + k * b_dim1;
00263 i__4 = i__ - 1;
00264 i__5 = i__ - 1 + k * x_dim1;
00265 i__6 = i__;
00266 i__7 = i__ + k * x_dim1;
00267 i__8 = i__;
00268 i__9 = i__ + 1 + k * x_dim1;
00269 tmp = (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = r_imag(&b[
00270 i__ + k * b_dim1]), dabs(r__2)) + ((r__3 = dl[
00271 i__4].r, dabs(r__3)) + (r__4 = r_imag(&dl[i__ - 1]
00272 ), dabs(r__4))) * ((r__5 = x[i__5].r, dabs(r__5))
00273 + (r__6 = r_imag(&x[i__ - 1 + k * x_dim1]), dabs(
00274 r__6))) + ((r__7 = d__[i__6].r, dabs(r__7)) + (
00275 r__8 = r_imag(&d__[i__]), dabs(r__8))) * ((r__9 =
00276 x[i__7].r, dabs(r__9)) + (r__10 = r_imag(&x[i__ +
00277 k * x_dim1]), dabs(r__10))) + ((r__11 = du[i__8]
00278 .r, dabs(r__11)) + (r__12 = r_imag(&du[i__]),
00279 dabs(r__12))) * ((r__13 = x[i__9].r, dabs(r__13))
00280 + (r__14 = r_imag(&x[i__ + 1 + k * x_dim1]), dabs(
00281 r__14)));
00282 axbi = dmin(axbi,tmp);
00283
00284 }
00285 i__2 = *n + k * b_dim1;
00286 i__3 = *n - 1;
00287 i__4 = *n - 1 + k * x_dim1;
00288 i__5 = *n;
00289 i__6 = *n + k * x_dim1;
00290 tmp = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[*n +
00291 k * b_dim1]), dabs(r__2)) + ((r__3 = dl[i__3].r, dabs(
00292 r__3)) + (r__4 = r_imag(&dl[*n - 1]), dabs(r__4))) * (
00293 (r__5 = x[i__4].r, dabs(r__5)) + (r__6 = r_imag(&x[*n
00294 - 1 + k * x_dim1]), dabs(r__6))) + ((r__7 = d__[i__5]
00295 .r, dabs(r__7)) + (r__8 = r_imag(&d__[*n]), dabs(r__8)
00296 )) * ((r__9 = x[i__6].r, dabs(r__9)) + (r__10 =
00297 r_imag(&x[*n + k * x_dim1]), dabs(r__10)));
00298 axbi = dmin(axbi,tmp);
00299 }
00300 } else {
00301 if (*n == 1) {
00302 i__2 = k * b_dim1 + 1;
00303 i__3 = k * x_dim1 + 1;
00304 axbi = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[k *
00305 b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, dabs(
00306 r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * ((
00307 r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[k *
00308 x_dim1 + 1]), dabs(r__6)));
00309 } else {
00310 i__2 = k * b_dim1 + 1;
00311 i__3 = k * x_dim1 + 1;
00312 i__4 = k * x_dim1 + 2;
00313 axbi = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[k *
00314 b_dim1 + 1]), dabs(r__2)) + ((r__3 = d__[1].r, dabs(
00315 r__3)) + (r__4 = r_imag(&d__[1]), dabs(r__4))) * ((
00316 r__5 = x[i__3].r, dabs(r__5)) + (r__6 = r_imag(&x[k *
00317 x_dim1 + 1]), dabs(r__6))) + ((r__7 = dl[1].r, dabs(
00318 r__7)) + (r__8 = r_imag(&dl[1]), dabs(r__8))) * ((
00319 r__9 = x[i__4].r, dabs(r__9)) + (r__10 = r_imag(&x[k *
00320 x_dim1 + 2]), dabs(r__10)));
00321 i__2 = *n - 1;
00322 for (i__ = 2; i__ <= i__2; ++i__) {
00323 i__3 = i__ + k * b_dim1;
00324 i__4 = i__ - 1;
00325 i__5 = i__ - 1 + k * x_dim1;
00326 i__6 = i__;
00327 i__7 = i__ + k * x_dim1;
00328 i__8 = i__;
00329 i__9 = i__ + 1 + k * x_dim1;
00330 tmp = (r__1 = b[i__3].r, dabs(r__1)) + (r__2 = r_imag(&b[
00331 i__ + k * b_dim1]), dabs(r__2)) + ((r__3 = du[
00332 i__4].r, dabs(r__3)) + (r__4 = r_imag(&du[i__ - 1]
00333 ), dabs(r__4))) * ((r__5 = x[i__5].r, dabs(r__5))
00334 + (r__6 = r_imag(&x[i__ - 1 + k * x_dim1]), dabs(
00335 r__6))) + ((r__7 = d__[i__6].r, dabs(r__7)) + (
00336 r__8 = r_imag(&d__[i__]), dabs(r__8))) * ((r__9 =
00337 x[i__7].r, dabs(r__9)) + (r__10 = r_imag(&x[i__ +
00338 k * x_dim1]), dabs(r__10))) + ((r__11 = dl[i__8]
00339 .r, dabs(r__11)) + (r__12 = r_imag(&dl[i__]),
00340 dabs(r__12))) * ((r__13 = x[i__9].r, dabs(r__13))
00341 + (r__14 = r_imag(&x[i__ + 1 + k * x_dim1]), dabs(
00342 r__14)));
00343 axbi = dmin(axbi,tmp);
00344
00345 }
00346 i__2 = *n + k * b_dim1;
00347 i__3 = *n - 1;
00348 i__4 = *n - 1 + k * x_dim1;
00349 i__5 = *n;
00350 i__6 = *n + k * x_dim1;
00351 tmp = (r__1 = b[i__2].r, dabs(r__1)) + (r__2 = r_imag(&b[*n +
00352 k * b_dim1]), dabs(r__2)) + ((r__3 = du[i__3].r, dabs(
00353 r__3)) + (r__4 = r_imag(&du[*n - 1]), dabs(r__4))) * (
00354 (r__5 = x[i__4].r, dabs(r__5)) + (r__6 = r_imag(&x[*n
00355 - 1 + k * x_dim1]), dabs(r__6))) + ((r__7 = d__[i__5]
00356 .r, dabs(r__7)) + (r__8 = r_imag(&d__[*n]), dabs(r__8)
00357 )) * ((r__9 = x[i__6].r, dabs(r__9)) + (r__10 =
00358 r_imag(&x[*n + k * x_dim1]), dabs(r__10)));
00359 axbi = dmin(axbi,tmp);
00360 }
00361 }
00362
00363 r__1 = axbi, r__2 = nz * unfl;
00364 tmp = berr[k] / (nz * eps + nz * unfl / dmax(r__1,r__2));
00365 if (k == 1) {
00366 reslts[2] = tmp;
00367 } else {
00368 reslts[2] = dmax(reslts[2],tmp);
00369 }
00370
00371 }
00372
00373 return 0;
00374
00375
00376
00377 }