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 complex c_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020 static integer c__1 = 1;
00021
00022 int chbt21_(char *uplo, integer *n, integer *ka, integer *ks,
00023 complex *a, integer *lda, real *d__, real *e, complex *u, integer *
00024 ldu, complex *work, real *rwork, real *result)
00025 {
00026
00027 integer a_dim1, a_offset, u_dim1, u_offset, i__1, i__2, i__3, i__4;
00028 real r__1, r__2;
00029 complex q__1, q__2;
00030
00031
00032 integer j, jc, jr, ika;
00033 real ulp;
00034 extern int chpr_(char *, integer *, real *, complex *,
00035 integer *, complex *);
00036 real unfl;
00037 extern int chpr2_(char *, integer *, complex *, complex *
00038 , integer *, complex *, integer *, complex *), cgemm_(
00039 char *, char *, integer *, integer *, integer *, complex *,
00040 complex *, integer *, complex *, integer *, complex *, complex *,
00041 integer *);
00042 extern logical lsame_(char *, char *);
00043 real anorm;
00044 char cuplo[1];
00045 logical lower;
00046 real wnorm;
00047 extern doublereal clanhb_(char *, char *, integer *, integer *, complex *,
00048 integer *, real *), clange_(char *, integer *,
00049 integer *, complex *, integer *, real *), clanhp_(char *,
00050 char *, integer *, complex *, real *), slamch_(
00051 char *);
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 a_dim1 = *lda;
00154 a_offset = 1 + a_dim1;
00155 a -= a_offset;
00156 --d__;
00157 --e;
00158 u_dim1 = *ldu;
00159 u_offset = 1 + u_dim1;
00160 u -= u_offset;
00161 --work;
00162 --rwork;
00163 --result;
00164
00165
00166 result[1] = 0.f;
00167 result[2] = 0.f;
00168 if (*n <= 0) {
00169 return 0;
00170 }
00171
00172
00173
00174 i__3 = *n - 1;
00175 i__1 = 0, i__2 = min(i__3,*ka);
00176 ika = max(i__1,i__2);
00177
00178 if (lsame_(uplo, "U")) {
00179 lower = FALSE_;
00180 *(unsigned char *)cuplo = 'U';
00181 } else {
00182 lower = TRUE_;
00183 *(unsigned char *)cuplo = 'L';
00184 }
00185
00186 unfl = slamch_("Safe minimum");
00187 ulp = slamch_("Epsilon") * slamch_("Base");
00188
00189
00190
00191
00192
00193
00194
00195
00196 r__1 = clanhb_("1", cuplo, n, &ika, &a[a_offset], lda, &rwork[1]);
00197 anorm = dmax(r__1,unfl);
00198
00199
00200
00201
00202
00203 j = 0;
00204 i__1 = *n;
00205 for (jc = 1; jc <= i__1; ++jc) {
00206 if (lower) {
00207
00208 i__3 = ika + 1, i__4 = *n + 1 - jc;
00209 i__2 = min(i__3,i__4);
00210 for (jr = 1; jr <= i__2; ++jr) {
00211 ++j;
00212 i__3 = j;
00213 i__4 = jr + jc * a_dim1;
00214 work[i__3].r = a[i__4].r, work[i__3].i = a[i__4].i;
00215
00216 }
00217 i__2 = *n + 1 - jc;
00218 for (jr = ika + 2; jr <= i__2; ++jr) {
00219 ++j;
00220 i__3 = j;
00221 work[i__3].r = 0.f, work[i__3].i = 0.f;
00222
00223 }
00224 } else {
00225 i__2 = jc;
00226 for (jr = ika + 2; jr <= i__2; ++jr) {
00227 ++j;
00228 i__3 = j;
00229 work[i__3].r = 0.f, work[i__3].i = 0.f;
00230
00231 }
00232
00233 i__2 = ika, i__3 = jc - 1;
00234 for (jr = min(i__2,i__3); jr >= 0; --jr) {
00235 ++j;
00236 i__2 = j;
00237 i__3 = ika + 1 - jr + jc * a_dim1;
00238 work[i__2].r = a[i__3].r, work[i__2].i = a[i__3].i;
00239
00240 }
00241 }
00242
00243 }
00244
00245 i__1 = *n;
00246 for (j = 1; j <= i__1; ++j) {
00247 r__1 = -d__[j];
00248 chpr_(cuplo, n, &r__1, &u[j * u_dim1 + 1], &c__1, &work[1])
00249 ;
00250
00251 }
00252
00253 if (*n > 1 && *ks == 1) {
00254 i__1 = *n - 1;
00255 for (j = 1; j <= i__1; ++j) {
00256 i__2 = j;
00257 q__2.r = e[i__2], q__2.i = 0.f;
00258 q__1.r = -q__2.r, q__1.i = -q__2.i;
00259 chpr2_(cuplo, n, &q__1, &u[j * u_dim1 + 1], &c__1, &u[(j + 1) *
00260 u_dim1 + 1], &c__1, &work[1]);
00261
00262 }
00263 }
00264 wnorm = clanhp_("1", cuplo, n, &work[1], &rwork[1]);
00265
00266 if (anorm > wnorm) {
00267 result[1] = wnorm / anorm / (*n * ulp);
00268 } else {
00269 if (anorm < 1.f) {
00270
00271 r__1 = wnorm, r__2 = *n * anorm;
00272 result[1] = dmin(r__1,r__2) / anorm / (*n * ulp);
00273 } else {
00274
00275 r__1 = wnorm / anorm, r__2 = (real) (*n);
00276 result[1] = dmin(r__1,r__2) / (*n * ulp);
00277 }
00278 }
00279
00280
00281
00282
00283
00284 cgemm_("N", "C", n, n, n, &c_b2, &u[u_offset], ldu, &u[u_offset], ldu, &
00285 c_b1, &work[1], n);
00286
00287 i__1 = *n;
00288 for (j = 1; j <= i__1; ++j) {
00289 i__2 = (*n + 1) * (j - 1) + 1;
00290 i__3 = (*n + 1) * (j - 1) + 1;
00291 q__1.r = work[i__3].r - 1.f, q__1.i = work[i__3].i - 0.f;
00292 work[i__2].r = q__1.r, work[i__2].i = q__1.i;
00293
00294 }
00295
00296
00297 r__1 = clange_("1", n, n, &work[1], n, &rwork[1]), r__2 = (
00298 real) (*n);
00299 result[2] = dmin(r__1,r__2) / (*n * ulp);
00300
00301 return 0;
00302
00303
00304
00305 }