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
00020 int zunt03_(char *rc, integer *mu, integer *mv, integer *n,
00021 integer *k, doublecomplex *u, integer *ldu, doublecomplex *v, integer
00022 *ldv, doublecomplex *work, integer *lwork, doublereal *rwork,
00023 doublereal *result, integer *info)
00024 {
00025
00026 integer u_dim1, u_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;
00027 doublereal d__1, d__2;
00028 doublecomplex z__1, z__2;
00029
00030
00031 double z_abs(doublecomplex *);
00032 void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00033
00034
00035 integer i__, j;
00036 doublecomplex s, su, sv;
00037 integer irc, lmx;
00038 doublereal ulp, res1, res2;
00039 extern logical lsame_(char *, char *);
00040 extern int zunt01_(char *, integer *, integer *,
00041 doublecomplex *, integer *, doublecomplex *, integer *,
00042 doublereal *, doublereal *);
00043 extern doublereal dlamch_(char *);
00044 extern int xerbla_(char *, integer *);
00045 extern integer izamax_(integer *, doublecomplex *, integer *);
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
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 u_dim1 = *ldu;
00155 u_offset = 1 + u_dim1;
00156 u -= u_offset;
00157 v_dim1 = *ldv;
00158 v_offset = 1 + v_dim1;
00159 v -= v_offset;
00160 --work;
00161 --rwork;
00162
00163
00164 *info = 0;
00165 if (lsame_(rc, "R")) {
00166 irc = 0;
00167 } else if (lsame_(rc, "C")) {
00168 irc = 1;
00169 } else {
00170 irc = -1;
00171 }
00172 if (irc == -1) {
00173 *info = -1;
00174 } else if (*mu < 0) {
00175 *info = -2;
00176 } else if (*mv < 0) {
00177 *info = -3;
00178 } else if (*n < 0) {
00179 *info = -4;
00180 } else if (*k < 0 || *k > max(*mu,*mv)) {
00181 *info = -5;
00182 } else if (irc == 0 && *ldu < max(1,*mu) || irc == 1 && *ldu < max(1,*n))
00183 {
00184 *info = -7;
00185 } else if (irc == 0 && *ldv < max(1,*mv) || irc == 1 && *ldv < max(1,*n))
00186 {
00187 *info = -9;
00188 }
00189 if (*info != 0) {
00190 i__1 = -(*info);
00191 xerbla_("ZUNT03", &i__1);
00192 return 0;
00193 }
00194
00195
00196
00197 *result = 0.;
00198 if (*mu == 0 || *mv == 0 || *n == 0) {
00199 return 0;
00200 }
00201
00202
00203
00204 ulp = dlamch_("Precision");
00205
00206 if (irc == 0) {
00207
00208
00209
00210 res1 = 0.;
00211 i__1 = *k;
00212 for (i__ = 1; i__ <= i__1; ++i__) {
00213 lmx = izamax_(n, &u[i__ + u_dim1], ldu);
00214 i__2 = i__ + lmx * v_dim1;
00215 if (v[i__2].r == 0. && v[i__2].i == 0.) {
00216 sv.r = 1., sv.i = 0.;
00217 } else {
00218 d__1 = z_abs(&v[i__ + lmx * v_dim1]);
00219 z__2.r = d__1, z__2.i = 0.;
00220 z_div(&z__1, &z__2, &v[i__ + lmx * v_dim1]);
00221 sv.r = z__1.r, sv.i = z__1.i;
00222 }
00223 i__2 = i__ + lmx * u_dim1;
00224 if (u[i__2].r == 0. && u[i__2].i == 0.) {
00225 su.r = 1., su.i = 0.;
00226 } else {
00227 d__1 = z_abs(&u[i__ + lmx * u_dim1]);
00228 z__2.r = d__1, z__2.i = 0.;
00229 z_div(&z__1, &z__2, &u[i__ + lmx * u_dim1]);
00230 su.r = z__1.r, su.i = z__1.i;
00231 }
00232 z_div(&z__1, &sv, &su);
00233 s.r = z__1.r, s.i = z__1.i;
00234 i__2 = *n;
00235 for (j = 1; j <= i__2; ++j) {
00236
00237 i__3 = i__ + j * u_dim1;
00238 i__4 = i__ + j * v_dim1;
00239 z__2.r = s.r * v[i__4].r - s.i * v[i__4].i, z__2.i = s.r * v[
00240 i__4].i + s.i * v[i__4].r;
00241 z__1.r = u[i__3].r - z__2.r, z__1.i = u[i__3].i - z__2.i;
00242 d__1 = res1, d__2 = z_abs(&z__1);
00243 res1 = max(d__1,d__2);
00244
00245 }
00246
00247 }
00248 res1 /= (doublereal) (*n) * ulp;
00249
00250
00251
00252 zunt01_("Rows", mv, n, &v[v_offset], ldv, &work[1], lwork, &rwork[1],
00253 &res2);
00254
00255 } else {
00256
00257
00258
00259 res1 = 0.;
00260 i__1 = *k;
00261 for (i__ = 1; i__ <= i__1; ++i__) {
00262 lmx = izamax_(n, &u[i__ * u_dim1 + 1], &c__1);
00263 i__2 = lmx + i__ * v_dim1;
00264 if (v[i__2].r == 0. && v[i__2].i == 0.) {
00265 sv.r = 1., sv.i = 0.;
00266 } else {
00267 d__1 = z_abs(&v[lmx + i__ * v_dim1]);
00268 z__2.r = d__1, z__2.i = 0.;
00269 z_div(&z__1, &z__2, &v[lmx + i__ * v_dim1]);
00270 sv.r = z__1.r, sv.i = z__1.i;
00271 }
00272 i__2 = lmx + i__ * u_dim1;
00273 if (u[i__2].r == 0. && u[i__2].i == 0.) {
00274 su.r = 1., su.i = 0.;
00275 } else {
00276 d__1 = z_abs(&u[lmx + i__ * u_dim1]);
00277 z__2.r = d__1, z__2.i = 0.;
00278 z_div(&z__1, &z__2, &u[lmx + i__ * u_dim1]);
00279 su.r = z__1.r, su.i = z__1.i;
00280 }
00281 z_div(&z__1, &sv, &su);
00282 s.r = z__1.r, s.i = z__1.i;
00283 i__2 = *n;
00284 for (j = 1; j <= i__2; ++j) {
00285
00286 i__3 = j + i__ * u_dim1;
00287 i__4 = j + i__ * v_dim1;
00288 z__2.r = s.r * v[i__4].r - s.i * v[i__4].i, z__2.i = s.r * v[
00289 i__4].i + s.i * v[i__4].r;
00290 z__1.r = u[i__3].r - z__2.r, z__1.i = u[i__3].i - z__2.i;
00291 d__1 = res1, d__2 = z_abs(&z__1);
00292 res1 = max(d__1,d__2);
00293
00294 }
00295
00296 }
00297 res1 /= (doublereal) (*n) * ulp;
00298
00299
00300
00301 zunt01_("Columns", n, mv, &v[v_offset], ldv, &work[1], lwork, &rwork[
00302 1], &res2);
00303 }
00304
00305
00306 d__1 = max(res1,res2), d__2 = 1. / ulp;
00307 *result = min(d__1,d__2);
00308 return 0;
00309
00310
00311
00312 }