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 doublecomplex c_b6 = {-1.,-0.};
00019 static integer c__1 = 1;
00020 static doublecomplex c_b9 = {0.,0.};
00021
00022 int zbdt03_(char *uplo, integer *n, integer *kd, doublereal *
00023 d__, doublereal *e, doublecomplex *u, integer *ldu, doublereal *s,
00024 doublecomplex *vt, integer *ldvt, doublecomplex *work, doublereal *
00025 resid)
00026 {
00027
00028 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2, i__3, i__4,
00029 i__5;
00030 doublereal d__1, d__2, d__3, d__4;
00031 doublecomplex z__1;
00032
00033
00034 integer i__, j;
00035 doublereal eps;
00036 extern logical lsame_(char *, char *);
00037 doublereal bnorm;
00038 extern int zgemv_(char *, integer *, integer *,
00039 doublecomplex *, doublecomplex *, integer *, doublecomplex *,
00040 integer *, doublecomplex *, doublecomplex *, integer *);
00041 extern doublereal dlamch_(char *);
00042 extern integer idamax_(integer *, doublereal *, integer *);
00043 extern doublereal dzasum_(integer *, doublecomplex *, integer *);
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 --d__;
00131 --e;
00132 u_dim1 = *ldu;
00133 u_offset = 1 + u_dim1;
00134 u -= u_offset;
00135 --s;
00136 vt_dim1 = *ldvt;
00137 vt_offset = 1 + vt_dim1;
00138 vt -= vt_offset;
00139 --work;
00140
00141
00142 *resid = 0.;
00143 if (*n <= 0) {
00144 return 0;
00145 }
00146
00147
00148
00149 bnorm = 0.;
00150 if (*kd >= 1) {
00151
00152
00153
00154 if (lsame_(uplo, "U")) {
00155
00156
00157
00158 i__1 = *n;
00159 for (j = 1; j <= i__1; ++j) {
00160 i__2 = *n;
00161 for (i__ = 1; i__ <= i__2; ++i__) {
00162 i__3 = *n + i__;
00163 i__4 = i__;
00164 i__5 = i__ + j * vt_dim1;
00165 z__1.r = s[i__4] * vt[i__5].r, z__1.i = s[i__4] * vt[i__5]
00166 .i;
00167 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00168
00169 }
00170 zgemv_("No transpose", n, n, &c_b6, &u[u_offset], ldu, &work[*
00171 n + 1], &c__1, &c_b9, &work[1], &c__1);
00172 i__2 = j;
00173 i__3 = j;
00174 i__4 = j;
00175 z__1.r = work[i__3].r + d__[i__4], z__1.i = work[i__3].i;
00176 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00177 if (j > 1) {
00178 i__2 = j - 1;
00179 i__3 = j - 1;
00180 i__4 = j - 1;
00181 z__1.r = work[i__3].r + e[i__4], z__1.i = work[i__3].i;
00182 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00183
00184 d__3 = bnorm, d__4 = (d__1 = d__[j], abs(d__1)) + (d__2 =
00185 e[j - 1], abs(d__2));
00186 bnorm = max(d__3,d__4);
00187 } else {
00188
00189 d__2 = bnorm, d__3 = (d__1 = d__[j], abs(d__1));
00190 bnorm = max(d__2,d__3);
00191 }
00192
00193 d__1 = *resid, d__2 = dzasum_(n, &work[1], &c__1);
00194 *resid = max(d__1,d__2);
00195
00196 }
00197 } else {
00198
00199
00200
00201 i__1 = *n;
00202 for (j = 1; j <= i__1; ++j) {
00203 i__2 = *n;
00204 for (i__ = 1; i__ <= i__2; ++i__) {
00205 i__3 = *n + i__;
00206 i__4 = i__;
00207 i__5 = i__ + j * vt_dim1;
00208 z__1.r = s[i__4] * vt[i__5].r, z__1.i = s[i__4] * vt[i__5]
00209 .i;
00210 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00211
00212 }
00213 zgemv_("No transpose", n, n, &c_b6, &u[u_offset], ldu, &work[*
00214 n + 1], &c__1, &c_b9, &work[1], &c__1);
00215 i__2 = j;
00216 i__3 = j;
00217 i__4 = j;
00218 z__1.r = work[i__3].r + d__[i__4], z__1.i = work[i__3].i;
00219 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00220 if (j < *n) {
00221 i__2 = j + 1;
00222 i__3 = j + 1;
00223 i__4 = j;
00224 z__1.r = work[i__3].r + e[i__4], z__1.i = work[i__3].i;
00225 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00226
00227 d__3 = bnorm, d__4 = (d__1 = d__[j], abs(d__1)) + (d__2 =
00228 e[j], abs(d__2));
00229 bnorm = max(d__3,d__4);
00230 } else {
00231
00232 d__2 = bnorm, d__3 = (d__1 = d__[j], abs(d__1));
00233 bnorm = max(d__2,d__3);
00234 }
00235
00236 d__1 = *resid, d__2 = dzasum_(n, &work[1], &c__1);
00237 *resid = max(d__1,d__2);
00238
00239 }
00240 }
00241 } else {
00242
00243
00244
00245 i__1 = *n;
00246 for (j = 1; j <= i__1; ++j) {
00247 i__2 = *n;
00248 for (i__ = 1; i__ <= i__2; ++i__) {
00249 i__3 = *n + i__;
00250 i__4 = i__;
00251 i__5 = i__ + j * vt_dim1;
00252 z__1.r = s[i__4] * vt[i__5].r, z__1.i = s[i__4] * vt[i__5].i;
00253 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00254
00255 }
00256 zgemv_("No transpose", n, n, &c_b6, &u[u_offset], ldu, &work[*n +
00257 1], &c__1, &c_b9, &work[1], &c__1);
00258 i__2 = j;
00259 i__3 = j;
00260 i__4 = j;
00261 z__1.r = work[i__3].r + d__[i__4], z__1.i = work[i__3].i;
00262 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00263
00264 d__1 = *resid, d__2 = dzasum_(n, &work[1], &c__1);
00265 *resid = max(d__1,d__2);
00266
00267 }
00268 j = idamax_(n, &d__[1], &c__1);
00269 bnorm = (d__1 = d__[j], abs(d__1));
00270 }
00271
00272
00273
00274 eps = dlamch_("Precision");
00275
00276 if (bnorm <= 0.) {
00277 if (*resid != 0.) {
00278 *resid = 1. / eps;
00279 }
00280 } else {
00281 if (bnorm >= *resid) {
00282 *resid = *resid / bnorm / ((doublereal) (*n) * eps);
00283 } else {
00284 if (bnorm < 1.) {
00285
00286 d__1 = *resid, d__2 = (doublereal) (*n) * bnorm;
00287 *resid = min(d__1,d__2) / bnorm / ((doublereal) (*n) * eps);
00288 } else {
00289
00290 d__1 = *resid / bnorm, d__2 = (doublereal) (*n);
00291 *resid = min(d__1,d__2) / ((doublereal) (*n) * eps);
00292 }
00293 }
00294 }
00295
00296 return 0;
00297
00298
00299
00300 }