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 static doublecomplex c_b7 = {-1.,-0.};
00020 static doublecomplex c_b10 = {1.,0.};
00021
00022 int zbdt01_(integer *m, integer *n, integer *kd,
00023 doublecomplex *a, integer *lda, doublecomplex *q, integer *ldq,
00024 doublereal *d__, doublereal *e, doublecomplex *pt, integer *ldpt,
00025 doublecomplex *work, doublereal *rwork, doublereal *resid)
00026 {
00027
00028 integer a_dim1, a_offset, pt_dim1, pt_offset, q_dim1, q_offset, i__1,
00029 i__2, i__3, i__4, i__5, i__6, i__7;
00030 doublereal d__1, d__2;
00031 doublecomplex z__1, z__2, z__3;
00032
00033
00034 integer i__, j;
00035 doublereal eps, anorm;
00036 extern int zgemv_(char *, integer *, integer *,
00037 doublecomplex *, doublecomplex *, integer *, doublecomplex *,
00038 integer *, doublecomplex *, doublecomplex *, integer *),
00039 zcopy_(integer *, doublecomplex *, integer *, doublecomplex *,
00040 integer *);
00041 extern doublereal dlamch_(char *), zlange_(char *, integer *,
00042 integer *, doublecomplex *, integer *, doublereal *),
00043 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
00131
00132
00133
00134
00135 a_dim1 = *lda;
00136 a_offset = 1 + a_dim1;
00137 a -= a_offset;
00138 q_dim1 = *ldq;
00139 q_offset = 1 + q_dim1;
00140 q -= q_offset;
00141 --d__;
00142 --e;
00143 pt_dim1 = *ldpt;
00144 pt_offset = 1 + pt_dim1;
00145 pt -= pt_offset;
00146 --work;
00147 --rwork;
00148
00149
00150 if (*m <= 0 || *n <= 0) {
00151 *resid = 0.;
00152 return 0;
00153 }
00154
00155
00156
00157 *resid = 0.;
00158 if (*kd != 0) {
00159
00160
00161
00162 if (*kd != 0 && *m >= *n) {
00163
00164
00165
00166 i__1 = *n;
00167 for (j = 1; j <= i__1; ++j) {
00168 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00169 i__2 = *n - 1;
00170 for (i__ = 1; i__ <= i__2; ++i__) {
00171 i__3 = *m + i__;
00172 i__4 = i__;
00173 i__5 = i__ + j * pt_dim1;
00174 z__2.r = d__[i__4] * pt[i__5].r, z__2.i = d__[i__4] * pt[
00175 i__5].i;
00176 i__6 = i__;
00177 i__7 = i__ + 1 + j * pt_dim1;
00178 z__3.r = e[i__6] * pt[i__7].r, z__3.i = e[i__6] * pt[i__7]
00179 .i;
00180 z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
00181 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00182
00183 }
00184 i__2 = *m + *n;
00185 i__3 = *n;
00186 i__4 = *n + j * pt_dim1;
00187 z__1.r = d__[i__3] * pt[i__4].r, z__1.i = d__[i__3] * pt[i__4]
00188 .i;
00189 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00190 zgemv_("No transpose", m, n, &c_b7, &q[q_offset], ldq, &work[*
00191 m + 1], &c__1, &c_b10, &work[1], &c__1);
00192
00193 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00194 *resid = max(d__1,d__2);
00195
00196 }
00197 } else if (*kd < 0) {
00198
00199
00200
00201 i__1 = *n;
00202 for (j = 1; j <= i__1; ++j) {
00203 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00204 i__2 = *m - 1;
00205 for (i__ = 1; i__ <= i__2; ++i__) {
00206 i__3 = *m + i__;
00207 i__4 = i__;
00208 i__5 = i__ + j * pt_dim1;
00209 z__2.r = d__[i__4] * pt[i__5].r, z__2.i = d__[i__4] * pt[
00210 i__5].i;
00211 i__6 = i__;
00212 i__7 = i__ + 1 + j * pt_dim1;
00213 z__3.r = e[i__6] * pt[i__7].r, z__3.i = e[i__6] * pt[i__7]
00214 .i;
00215 z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
00216 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00217
00218 }
00219 i__2 = *m + *m;
00220 i__3 = *m;
00221 i__4 = *m + j * pt_dim1;
00222 z__1.r = d__[i__3] * pt[i__4].r, z__1.i = d__[i__3] * pt[i__4]
00223 .i;
00224 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00225 zgemv_("No transpose", m, m, &c_b7, &q[q_offset], ldq, &work[*
00226 m + 1], &c__1, &c_b10, &work[1], &c__1);
00227
00228 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00229 *resid = max(d__1,d__2);
00230
00231 }
00232 } else {
00233
00234
00235
00236 i__1 = *n;
00237 for (j = 1; j <= i__1; ++j) {
00238 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00239 i__2 = *m + 1;
00240 i__3 = j * pt_dim1 + 1;
00241 z__1.r = d__[1] * pt[i__3].r, z__1.i = d__[1] * pt[i__3].i;
00242 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00243 i__2 = *m;
00244 for (i__ = 2; i__ <= i__2; ++i__) {
00245 i__3 = *m + i__;
00246 i__4 = i__ - 1;
00247 i__5 = i__ - 1 + j * pt_dim1;
00248 z__2.r = e[i__4] * pt[i__5].r, z__2.i = e[i__4] * pt[i__5]
00249 .i;
00250 i__6 = i__;
00251 i__7 = i__ + j * pt_dim1;
00252 z__3.r = d__[i__6] * pt[i__7].r, z__3.i = d__[i__6] * pt[
00253 i__7].i;
00254 z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
00255 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00256
00257 }
00258 zgemv_("No transpose", m, m, &c_b7, &q[q_offset], ldq, &work[*
00259 m + 1], &c__1, &c_b10, &work[1], &c__1);
00260
00261 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00262 *resid = max(d__1,d__2);
00263
00264 }
00265 }
00266 } else {
00267
00268
00269
00270 if (*m >= *n) {
00271 i__1 = *n;
00272 for (j = 1; j <= i__1; ++j) {
00273 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00274 i__2 = *n;
00275 for (i__ = 1; i__ <= i__2; ++i__) {
00276 i__3 = *m + i__;
00277 i__4 = i__;
00278 i__5 = i__ + j * pt_dim1;
00279 z__1.r = d__[i__4] * pt[i__5].r, z__1.i = d__[i__4] * pt[
00280 i__5].i;
00281 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00282
00283 }
00284 zgemv_("No transpose", m, n, &c_b7, &q[q_offset], ldq, &work[*
00285 m + 1], &c__1, &c_b10, &work[1], &c__1);
00286
00287 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00288 *resid = max(d__1,d__2);
00289
00290 }
00291 } else {
00292 i__1 = *n;
00293 for (j = 1; j <= i__1; ++j) {
00294 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00295 i__2 = *m;
00296 for (i__ = 1; i__ <= i__2; ++i__) {
00297 i__3 = *m + i__;
00298 i__4 = i__;
00299 i__5 = i__ + j * pt_dim1;
00300 z__1.r = d__[i__4] * pt[i__5].r, z__1.i = d__[i__4] * pt[
00301 i__5].i;
00302 work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00303
00304 }
00305 zgemv_("No transpose", m, m, &c_b7, &q[q_offset], ldq, &work[*
00306 m + 1], &c__1, &c_b10, &work[1], &c__1);
00307
00308 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00309 *resid = max(d__1,d__2);
00310
00311 }
00312 }
00313 }
00314
00315
00316
00317 anorm = zlange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00318 eps = dlamch_("Precision");
00319
00320 if (anorm <= 0.) {
00321 if (*resid != 0.) {
00322 *resid = 1. / eps;
00323 }
00324 } else {
00325 if (anorm >= *resid) {
00326 *resid = *resid / anorm / ((doublereal) (*n) * eps);
00327 } else {
00328 if (anorm < 1.) {
00329
00330 d__1 = *resid, d__2 = (doublereal) (*n) * anorm;
00331 *resid = min(d__1,d__2) / anorm / ((doublereal) (*n) * eps);
00332 } else {
00333
00334 d__1 = *resid / anorm, d__2 = (doublereal) (*n);
00335 *resid = min(d__1,d__2) / ((doublereal) (*n) * eps);
00336 }
00337 }
00338 }
00339
00340 return 0;
00341
00342
00343
00344 }