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