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 complex c_b12 = {-1.f,-0.f};
00020
00021 int cgbt01_(integer *m, integer *n, integer *kl, integer *ku,
00022 complex *a, integer *lda, complex *afac, integer *ldafac, integer *
00023 ipiv, complex *work, real *resid)
00024 {
00025
00026 integer a_dim1, a_offset, afac_dim1, afac_offset, i__1, i__2, i__3, i__4;
00027 real r__1, r__2;
00028
00029
00030 integer i__, j;
00031 complex t;
00032 integer i1, i2, kd, il, jl, ip, ju, iw, jua;
00033 real eps;
00034 integer lenj;
00035 real anorm;
00036 extern int ccopy_(integer *, complex *, integer *,
00037 complex *, integer *), caxpy_(integer *, complex *, complex *,
00038 integer *, complex *, integer *);
00039 extern doublereal slamch_(char *), scasum_(integer *, complex *,
00040 integer *);
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
00120
00121
00122 a_dim1 = *lda;
00123 a_offset = 1 + a_dim1;
00124 a -= a_offset;
00125 afac_dim1 = *ldafac;
00126 afac_offset = 1 + afac_dim1;
00127 afac -= afac_offset;
00128 --ipiv;
00129 --work;
00130
00131
00132 *resid = 0.f;
00133 if (*m <= 0 || *n <= 0) {
00134 return 0;
00135 }
00136
00137
00138
00139 eps = slamch_("Epsilon");
00140 kd = *ku + 1;
00141 anorm = 0.f;
00142 i__1 = *n;
00143 for (j = 1; j <= i__1; ++j) {
00144
00145 i__2 = kd + 1 - j;
00146 i1 = max(i__2,1);
00147
00148 i__2 = kd + *m - j, i__3 = *kl + kd;
00149 i2 = min(i__2,i__3);
00150 if (i2 >= i1) {
00151
00152 i__2 = i2 - i1 + 1;
00153 r__1 = anorm, r__2 = scasum_(&i__2, &a[i1 + j * a_dim1], &c__1);
00154 anorm = dmax(r__1,r__2);
00155 }
00156
00157 }
00158
00159
00160
00161 kd = *kl + *ku + 1;
00162 i__1 = *n;
00163 for (j = 1; j <= i__1; ++j) {
00164
00165
00166
00167
00168 i__2 = *kl + *ku, i__3 = j - 1;
00169 ju = min(i__2,i__3);
00170
00171 i__2 = *kl, i__3 = *m - j;
00172 jl = min(i__2,i__3);
00173 lenj = min(*m,j) - j + ju + 1;
00174 if (lenj > 0) {
00175 ccopy_(&lenj, &afac[kd - ju + j * afac_dim1], &c__1, &work[1], &
00176 c__1);
00177 i__2 = ju + jl + 1;
00178 for (i__ = lenj + 1; i__ <= i__2; ++i__) {
00179 i__3 = i__;
00180 work[i__3].r = 0.f, work[i__3].i = 0.f;
00181
00182 }
00183
00184
00185
00186
00187
00188 i__2 = *m - 1;
00189 i__3 = j - ju;
00190 for (i__ = min(i__2,j); i__ >= i__3; --i__) {
00191
00192 i__2 = *kl, i__4 = *m - i__;
00193 il = min(i__2,i__4);
00194 if (il > 0) {
00195 iw = i__ - j + ju + 1;
00196 i__2 = iw;
00197 t.r = work[i__2].r, t.i = work[i__2].i;
00198 caxpy_(&il, &t, &afac[kd + 1 + i__ * afac_dim1], &c__1, &
00199 work[iw + 1], &c__1);
00200 ip = ipiv[i__];
00201 if (i__ != ip) {
00202 ip = ip - j + ju + 1;
00203 i__2 = iw;
00204 i__4 = ip;
00205 work[i__2].r = work[i__4].r, work[i__2].i = work[i__4]
00206 .i;
00207 i__2 = ip;
00208 work[i__2].r = t.r, work[i__2].i = t.i;
00209 }
00210 }
00211
00212 }
00213
00214
00215
00216 jua = min(ju,*ku);
00217 if (jua + jl + 1 > 0) {
00218 i__3 = jua + jl + 1;
00219 caxpy_(&i__3, &c_b12, &a[*ku + 1 - jua + j * a_dim1], &c__1, &
00220 work[ju + 1 - jua], &c__1);
00221 }
00222
00223
00224
00225
00226 i__3 = ju + jl + 1;
00227 r__1 = *resid, r__2 = scasum_(&i__3, &work[1], &c__1);
00228 *resid = dmax(r__1,r__2);
00229 }
00230
00231 }
00232
00233
00234
00235 if (anorm <= 0.f) {
00236 if (*resid != 0.f) {
00237 *resid = 1.f / eps;
00238 }
00239 } else {
00240 *resid = *resid / (real) (*n) / anorm / eps;
00241 }
00242
00243 return 0;
00244
00245
00246
00247 }