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_b7 = {0.f,0.f};
00019 static complex c_b8 = {1.f,0.f};
00020 static real c_b10 = -1.f;
00021 static real c_b11 = 1.f;
00022 static integer c__1 = 1;
00023
00024 int cunt01_(char *rowcol, integer *m, integer *n, complex *u,
00025 integer *ldu, complex *work, integer *lwork, real *rwork, real *
00026 resid)
00027 {
00028
00029 integer u_dim1, u_offset, i__1, i__2;
00030 real r__1, r__2, r__3, r__4;
00031 complex q__1, q__2;
00032
00033
00034 double r_imag(complex *);
00035
00036
00037 integer i__, j, k;
00038 real eps;
00039 complex tmp;
00040 extern VOID cdotc_(complex *, integer *, complex *, integer
00041 *, complex *, integer *);
00042 extern int cherk_(char *, char *, integer *, integer *,
00043 real *, complex *, integer *, real *, complex *, integer *);
00044 extern logical lsame_(char *, char *);
00045 integer mnmin;
00046 extern doublereal slamch_(char *);
00047 extern int claset_(char *, integer *, integer *, complex
00048 *, complex *, complex *, integer *);
00049 extern doublereal clansy_(char *, char *, integer *, complex *, integer *,
00050 real *);
00051 integer ldwork;
00052 char transu[1];
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 u_dim1 = *ldu;
00142 u_offset = 1 + u_dim1;
00143 u -= u_offset;
00144 --work;
00145 --rwork;
00146
00147
00148 *resid = 0.f;
00149
00150
00151
00152 if (*m <= 0 || *n <= 0) {
00153 return 0;
00154 }
00155
00156 eps = slamch_("Precision");
00157 if (*m < *n || *m == *n && lsame_(rowcol, "R")) {
00158 *(unsigned char *)transu = 'N';
00159 k = *n;
00160 } else {
00161 *(unsigned char *)transu = 'C';
00162 k = *m;
00163 }
00164 mnmin = min(*m,*n);
00165
00166 if ((mnmin + 1) * mnmin <= *lwork) {
00167 ldwork = mnmin;
00168 } else {
00169 ldwork = 0;
00170 }
00171 if (ldwork > 0) {
00172
00173
00174
00175 claset_("Upper", &mnmin, &mnmin, &c_b7, &c_b8, &work[1], &ldwork);
00176 cherk_("Upper", transu, &mnmin, &k, &c_b10, &u[u_offset], ldu, &c_b11,
00177 &work[1], &ldwork);
00178
00179
00180
00181 *resid = clansy_("1", "Upper", &mnmin, &work[1], &ldwork, &rwork[1]);
00182 *resid = *resid / (real) k / eps;
00183 } else if (*(unsigned char *)transu == 'C') {
00184
00185
00186
00187 i__1 = *n;
00188 for (j = 1; j <= i__1; ++j) {
00189 i__2 = j;
00190 for (i__ = 1; i__ <= i__2; ++i__) {
00191 if (i__ != j) {
00192 tmp.r = 0.f, tmp.i = 0.f;
00193 } else {
00194 tmp.r = 1.f, tmp.i = 0.f;
00195 }
00196 cdotc_(&q__2, m, &u[i__ * u_dim1 + 1], &c__1, &u[j * u_dim1 +
00197 1], &c__1);
00198 q__1.r = tmp.r - q__2.r, q__1.i = tmp.i - q__2.i;
00199 tmp.r = q__1.r, tmp.i = q__1.i;
00200
00201 r__3 = *resid, r__4 = (r__1 = tmp.r, dabs(r__1)) + (r__2 =
00202 r_imag(&tmp), dabs(r__2));
00203 *resid = dmax(r__3,r__4);
00204
00205 }
00206
00207 }
00208 *resid = *resid / (real) (*m) / eps;
00209 } else {
00210
00211
00212
00213 i__1 = *m;
00214 for (j = 1; j <= i__1; ++j) {
00215 i__2 = j;
00216 for (i__ = 1; i__ <= i__2; ++i__) {
00217 if (i__ != j) {
00218 tmp.r = 0.f, tmp.i = 0.f;
00219 } else {
00220 tmp.r = 1.f, tmp.i = 0.f;
00221 }
00222 cdotc_(&q__2, n, &u[j + u_dim1], ldu, &u[i__ + u_dim1], ldu);
00223 q__1.r = tmp.r - q__2.r, q__1.i = tmp.i - q__2.i;
00224 tmp.r = q__1.r, tmp.i = q__1.i;
00225
00226 r__3 = *resid, r__4 = (r__1 = tmp.r, dabs(r__1)) + (r__2 =
00227 r_imag(&tmp), dabs(r__2));
00228 *resid = dmax(r__3,r__4);
00229
00230 }
00231
00232 }
00233 *resid = *resid / (real) (*n) / eps;
00234 }
00235 return 0;
00236
00237
00238
00239 }