00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016 int zgtt01_(integer *n, doublecomplex *dl, doublecomplex *
00017 d__, doublecomplex *du, doublecomplex *dlf, doublecomplex *df,
00018 doublecomplex *duf, doublecomplex *du2, integer *ipiv, doublecomplex *
00019 work, integer *ldwork, doublereal *rwork, doublereal *resid)
00020 {
00021
00022 integer work_dim1, work_offset, i__1, i__2, i__3, i__4;
00023 doublecomplex z__1;
00024
00025
00026 integer i__, j;
00027 doublecomplex li;
00028 integer ip;
00029 doublereal eps, anorm;
00030 integer lastj;
00031 extern int zswap_(integer *, doublecomplex *, integer *,
00032 doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *,
00033 doublecomplex *, integer *, doublecomplex *, integer *);
00034 extern doublereal dlamch_(char *), zlangt_(char *, integer *,
00035 doublecomplex *, doublecomplex *, doublecomplex *),
00036 zlanhs_(char *, integer *, doublecomplex *, integer *, doublereal
00037 *);
00038
00039
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 --dl;
00120 --d__;
00121 --du;
00122 --dlf;
00123 --df;
00124 --duf;
00125 --du2;
00126 --ipiv;
00127 work_dim1 = *ldwork;
00128 work_offset = 1 + work_dim1;
00129 work -= work_offset;
00130 --rwork;
00131
00132
00133 if (*n <= 0) {
00134 *resid = 0.;
00135 return 0;
00136 }
00137
00138 eps = dlamch_("Epsilon");
00139
00140
00141
00142 i__1 = *n;
00143 for (j = 1; j <= i__1; ++j) {
00144 i__2 = *n;
00145 for (i__ = 1; i__ <= i__2; ++i__) {
00146 i__3 = i__ + j * work_dim1;
00147 work[i__3].r = 0., work[i__3].i = 0.;
00148
00149 }
00150
00151 }
00152 i__1 = *n;
00153 for (i__ = 1; i__ <= i__1; ++i__) {
00154 if (i__ == 1) {
00155 i__2 = i__ + i__ * work_dim1;
00156 i__3 = i__;
00157 work[i__2].r = df[i__3].r, work[i__2].i = df[i__3].i;
00158 if (*n >= 2) {
00159 i__2 = i__ + (i__ + 1) * work_dim1;
00160 i__3 = i__;
00161 work[i__2].r = duf[i__3].r, work[i__2].i = duf[i__3].i;
00162 }
00163 if (*n >= 3) {
00164 i__2 = i__ + (i__ + 2) * work_dim1;
00165 i__3 = i__;
00166 work[i__2].r = du2[i__3].r, work[i__2].i = du2[i__3].i;
00167 }
00168 } else if (i__ == *n) {
00169 i__2 = i__ + i__ * work_dim1;
00170 i__3 = i__;
00171 work[i__2].r = df[i__3].r, work[i__2].i = df[i__3].i;
00172 } else {
00173 i__2 = i__ + i__ * work_dim1;
00174 i__3 = i__;
00175 work[i__2].r = df[i__3].r, work[i__2].i = df[i__3].i;
00176 i__2 = i__ + (i__ + 1) * work_dim1;
00177 i__3 = i__;
00178 work[i__2].r = duf[i__3].r, work[i__2].i = duf[i__3].i;
00179 if (i__ < *n - 1) {
00180 i__2 = i__ + (i__ + 2) * work_dim1;
00181 i__3 = i__;
00182 work[i__2].r = du2[i__3].r, work[i__2].i = du2[i__3].i;
00183 }
00184 }
00185
00186 }
00187
00188
00189
00190 lastj = *n;
00191 for (i__ = *n - 1; i__ >= 1; --i__) {
00192 i__1 = i__;
00193 li.r = dlf[i__1].r, li.i = dlf[i__1].i;
00194 i__1 = lastj - i__ + 1;
00195 zaxpy_(&i__1, &li, &work[i__ + i__ * work_dim1], ldwork, &work[i__ +
00196 1 + i__ * work_dim1], ldwork);
00197 ip = ipiv[i__];
00198 if (ip == i__) {
00199
00200 i__1 = i__ + 2;
00201 lastj = min(i__1,*n);
00202 } else {
00203 i__1 = lastj - i__ + 1;
00204 zswap_(&i__1, &work[i__ + i__ * work_dim1], ldwork, &work[i__ + 1
00205 + i__ * work_dim1], ldwork);
00206 }
00207
00208 }
00209
00210
00211
00212 i__1 = work_dim1 + 1;
00213 i__2 = work_dim1 + 1;
00214 z__1.r = work[i__2].r - d__[1].r, z__1.i = work[i__2].i - d__[1].i;
00215 work[i__1].r = z__1.r, work[i__1].i = z__1.i;
00216 if (*n > 1) {
00217 i__1 = (work_dim1 << 1) + 1;
00218 i__2 = (work_dim1 << 1) + 1;
00219 z__1.r = work[i__2].r - du[1].r, z__1.i = work[i__2].i - du[1].i;
00220 work[i__1].r = z__1.r, work[i__1].i = z__1.i;
00221 i__1 = *n + (*n - 1) * work_dim1;
00222 i__2 = *n + (*n - 1) * work_dim1;
00223 i__3 = *n - 1;
00224 z__1.r = work[i__2].r - dl[i__3].r, z__1.i = work[i__2].i - dl[i__3]
00225 .i;
00226 work[i__1].r = z__1.r, work[i__1].i = z__1.i;
00227 i__1 = *n + *n * work_dim1;
00228 i__2 = *n + *n * work_dim1;
00229 i__3 = *n;
00230 z__1.r = work[i__2].r - d__[i__3].r, z__1.i = work[i__2].i - d__[i__3]
00231 .i;
00232 work[i__1].r = z__1.r, work[i__1].i = z__1.i;
00233 i__1 = *n - 1;
00234 for (i__ = 2; i__ <= i__1; ++i__) {
00235 i__2 = i__ + (i__ - 1) * work_dim1;
00236 i__3 = i__ + (i__ - 1) * work_dim1;
00237 i__4 = i__ - 1;
00238 z__1.r = work[i__3].r - dl[i__4].r, z__1.i = work[i__3].i - dl[
00239 i__4].i;
00240 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00241 i__2 = i__ + i__ * work_dim1;
00242 i__3 = i__ + i__ * work_dim1;
00243 i__4 = i__;
00244 z__1.r = work[i__3].r - d__[i__4].r, z__1.i = work[i__3].i - d__[
00245 i__4].i;
00246 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00247 i__2 = i__ + (i__ + 1) * work_dim1;
00248 i__3 = i__ + (i__ + 1) * work_dim1;
00249 i__4 = i__;
00250 z__1.r = work[i__3].r - du[i__4].r, z__1.i = work[i__3].i - du[
00251 i__4].i;
00252 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00253
00254 }
00255 }
00256
00257
00258
00259 anorm = zlangt_("1", n, &dl[1], &d__[1], &du[1]);
00260
00261
00262
00263
00264 *resid = zlanhs_("1", n, &work[work_offset], ldwork, &rwork[1])
00265 ;
00266
00267
00268
00269 if (anorm <= 0.) {
00270 if (*resid != 0.) {
00271 *resid = 1. / eps;
00272 }
00273 } else {
00274 *resid = *resid / anorm / eps;
00275 }
00276
00277 return 0;
00278
00279
00280
00281 }