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 real c_b36 = -1.f;
00021 static real c_b37 = 1.f;
00022 static integer c__1 = 1;
00023
00024 int cgsvts_(integer *m, integer *p, integer *n, complex *a,
00025 complex *af, integer *lda, complex *b, complex *bf, integer *ldb,
00026 complex *u, integer *ldu, complex *v, integer *ldv, complex *q,
00027 integer *ldq, real *alpha, real *beta, complex *r__, integer *ldr,
00028 integer *iwork, complex *work, integer *lwork, real *rwork, real *
00029 result)
00030 {
00031
00032 integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, bf_dim1,
00033 bf_offset, q_dim1, q_offset, r_dim1, r_offset, u_dim1, u_offset,
00034 v_dim1, v_offset, i__1, i__2, i__3, i__4, i__5, i__6;
00035 real r__1;
00036 complex q__1, q__2;
00037
00038
00039 integer i__, j, k, l;
00040 real ulp;
00041 integer info;
00042 real unfl, temp;
00043 extern int cgemm_(char *, char *, integer *, integer *,
00044 integer *, complex *, complex *, integer *, complex *, integer *,
00045 complex *, complex *, integer *), cherk_(char *,
00046 char *, integer *, integer *, real *, complex *, integer *, real *
00047 , complex *, integer *);
00048 real resid, anorm, bnorm;
00049 extern int scopy_(integer *, real *, integer *, real *,
00050 integer *);
00051 extern doublereal clange_(char *, integer *, integer *, complex *,
00052 integer *, real *), clanhe_(char *, char *, integer *,
00053 complex *, integer *, real *), slamch_(char *);
00054 extern int clacpy_(char *, integer *, integer *, complex
00055 *, integer *, complex *, integer *), claset_(char *,
00056 integer *, integer *, complex *, complex *, complex *, integer *), cggsvd_(char *, char *, char *, integer *, integer *,
00057 integer *, integer *, integer *, complex *, integer *, complex *,
00058 integer *, real *, real *, complex *, integer *, complex *,
00059 integer *, complex *, integer *, complex *, real *, integer *,
00060 integer *);
00061 real ulpinv;
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
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179 af_dim1 = *lda;
00180 af_offset = 1 + af_dim1;
00181 af -= af_offset;
00182 a_dim1 = *lda;
00183 a_offset = 1 + a_dim1;
00184 a -= a_offset;
00185 bf_dim1 = *ldb;
00186 bf_offset = 1 + bf_dim1;
00187 bf -= bf_offset;
00188 b_dim1 = *ldb;
00189 b_offset = 1 + b_dim1;
00190 b -= b_offset;
00191 u_dim1 = *ldu;
00192 u_offset = 1 + u_dim1;
00193 u -= u_offset;
00194 v_dim1 = *ldv;
00195 v_offset = 1 + v_dim1;
00196 v -= v_offset;
00197 q_dim1 = *ldq;
00198 q_offset = 1 + q_dim1;
00199 q -= q_offset;
00200 --alpha;
00201 --beta;
00202 r_dim1 = *ldr;
00203 r_offset = 1 + r_dim1;
00204 r__ -= r_offset;
00205 --iwork;
00206 --work;
00207 --rwork;
00208 --result;
00209
00210
00211 ulp = slamch_("Precision");
00212 ulpinv = 1.f / ulp;
00213 unfl = slamch_("Safe minimum");
00214
00215
00216
00217 clacpy_("Full", m, n, &a[a_offset], lda, &af[af_offset], lda);
00218 clacpy_("Full", p, n, &b[b_offset], ldb, &bf[bf_offset], ldb);
00219
00220
00221 r__1 = clange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00222 anorm = dmax(r__1,unfl);
00223
00224 r__1 = clange_("1", p, n, &b[b_offset], ldb, &rwork[1]);
00225 bnorm = dmax(r__1,unfl);
00226
00227
00228
00229 cggsvd_("U", "V", "Q", m, n, p, &k, &l, &af[af_offset], lda, &bf[
00230 bf_offset], ldb, &alpha[1], &beta[1], &u[u_offset], ldu, &v[
00231 v_offset], ldv, &q[q_offset], ldq, &work[1], &rwork[1], &iwork[1],
00232 &info);
00233
00234
00235
00236
00237 i__2 = k + l;
00238 i__1 = min(i__2,*m);
00239 for (i__ = 1; i__ <= i__1; ++i__) {
00240 i__2 = k + l;
00241 for (j = i__; j <= i__2; ++j) {
00242 i__3 = i__ + j * r_dim1;
00243 i__4 = i__ + (*n - k - l + j) * af_dim1;
00244 r__[i__3].r = af[i__4].r, r__[i__3].i = af[i__4].i;
00245
00246 }
00247
00248 }
00249
00250 if (*m - k - l < 0) {
00251 i__1 = k + l;
00252 for (i__ = *m + 1; i__ <= i__1; ++i__) {
00253 i__2 = k + l;
00254 for (j = i__; j <= i__2; ++j) {
00255 i__3 = i__ + j * r_dim1;
00256 i__4 = i__ - k + (*n - k - l + j) * bf_dim1;
00257 r__[i__3].r = bf[i__4].r, r__[i__3].i = bf[i__4].i;
00258
00259 }
00260
00261 }
00262 }
00263
00264
00265
00266 cgemm_("No transpose", "No transpose", m, n, n, &c_b2, &a[a_offset], lda,
00267 &q[q_offset], ldq, &c_b1, &work[1], lda);
00268
00269 cgemm_("Conjugate transpose", "No transpose", m, n, m, &c_b2, &u[u_offset]
00270 , ldu, &work[1], lda, &c_b1, &a[a_offset], lda);
00271
00272 i__1 = k;
00273 for (i__ = 1; i__ <= i__1; ++i__) {
00274 i__2 = k + l;
00275 for (j = i__; j <= i__2; ++j) {
00276 i__3 = i__ + (*n - k - l + j) * a_dim1;
00277 i__4 = i__ + (*n - k - l + j) * a_dim1;
00278 i__5 = i__ + j * r_dim1;
00279 q__1.r = a[i__4].r - r__[i__5].r, q__1.i = a[i__4].i - r__[i__5]
00280 .i;
00281 a[i__3].r = q__1.r, a[i__3].i = q__1.i;
00282
00283 }
00284
00285 }
00286
00287
00288 i__2 = k + l;
00289 i__1 = min(i__2,*m);
00290 for (i__ = k + 1; i__ <= i__1; ++i__) {
00291 i__2 = k + l;
00292 for (j = i__; j <= i__2; ++j) {
00293 i__3 = i__ + (*n - k - l + j) * a_dim1;
00294 i__4 = i__ + (*n - k - l + j) * a_dim1;
00295 i__5 = i__;
00296 i__6 = i__ + j * r_dim1;
00297 q__2.r = alpha[i__5] * r__[i__6].r, q__2.i = alpha[i__5] * r__[
00298 i__6].i;
00299 q__1.r = a[i__4].r - q__2.r, q__1.i = a[i__4].i - q__2.i;
00300 a[i__3].r = q__1.r, a[i__3].i = q__1.i;
00301
00302 }
00303
00304 }
00305
00306
00307
00308 resid = clange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00309 if (anorm > 0.f) {
00310
00311 i__1 = max(1,*m);
00312 result[1] = resid / (real) max(i__1,*n) / anorm / ulp;
00313 } else {
00314 result[1] = 0.f;
00315 }
00316
00317
00318
00319 cgemm_("No transpose", "No transpose", p, n, n, &c_b2, &b[b_offset], ldb,
00320 &q[q_offset], ldq, &c_b1, &work[1], ldb);
00321
00322 cgemm_("Conjugate transpose", "No transpose", p, n, p, &c_b2, &v[v_offset]
00323 , ldv, &work[1], ldb, &c_b1, &b[b_offset], ldb);
00324
00325 i__1 = l;
00326 for (i__ = 1; i__ <= i__1; ++i__) {
00327 i__2 = l;
00328 for (j = i__; j <= i__2; ++j) {
00329 i__3 = i__ + (*n - l + j) * b_dim1;
00330 i__4 = i__ + (*n - l + j) * b_dim1;
00331 i__5 = k + i__;
00332 i__6 = k + i__ + (k + j) * r_dim1;
00333 q__2.r = beta[i__5] * r__[i__6].r, q__2.i = beta[i__5] * r__[i__6]
00334 .i;
00335 q__1.r = b[i__4].r - q__2.r, q__1.i = b[i__4].i - q__2.i;
00336 b[i__3].r = q__1.r, b[i__3].i = q__1.i;
00337
00338 }
00339
00340 }
00341
00342
00343
00344 resid = clange_("1", p, n, &b[b_offset], ldb, &rwork[1]);
00345 if (bnorm > 0.f) {
00346
00347 i__1 = max(1,*p);
00348 result[2] = resid / (real) max(i__1,*n) / bnorm / ulp;
00349 } else {
00350 result[2] = 0.f;
00351 }
00352
00353
00354
00355 claset_("Full", m, m, &c_b1, &c_b2, &work[1], ldq);
00356 cherk_("Upper", "Conjugate transpose", m, m, &c_b36, &u[u_offset], ldu, &
00357 c_b37, &work[1], ldu);
00358
00359
00360
00361 resid = clanhe_("1", "Upper", m, &work[1], ldu, &rwork[1]);
00362 result[3] = resid / (real) max(1,*m) / ulp;
00363
00364
00365
00366 claset_("Full", p, p, &c_b1, &c_b2, &work[1], ldv);
00367 cherk_("Upper", "Conjugate transpose", p, p, &c_b36, &v[v_offset], ldv, &
00368 c_b37, &work[1], ldv);
00369
00370
00371
00372 resid = clanhe_("1", "Upper", p, &work[1], ldv, &rwork[1]);
00373 result[4] = resid / (real) max(1,*p) / ulp;
00374
00375
00376
00377 claset_("Full", n, n, &c_b1, &c_b2, &work[1], ldq);
00378 cherk_("Upper", "Conjugate transpose", n, n, &c_b36, &q[q_offset], ldq, &
00379 c_b37, &work[1], ldq);
00380
00381
00382
00383 resid = clanhe_("1", "Upper", n, &work[1], ldq, &rwork[1]);
00384 result[5] = resid / (real) max(1,*n) / ulp;
00385
00386
00387
00388 scopy_(n, &alpha[1], &c__1, &rwork[1], &c__1);
00389
00390 i__2 = k + l;
00391 i__1 = min(i__2,*m);
00392 for (i__ = k + 1; i__ <= i__1; ++i__) {
00393 j = iwork[i__];
00394 if (i__ != j) {
00395 temp = rwork[i__];
00396 rwork[i__] = rwork[j];
00397 rwork[j] = temp;
00398 }
00399
00400 }
00401
00402 result[6] = 0.f;
00403
00404 i__2 = k + l;
00405 i__1 = min(i__2,*m) - 1;
00406 for (i__ = k + 1; i__ <= i__1; ++i__) {
00407 if (rwork[i__] < rwork[i__ + 1]) {
00408 result[6] = ulpinv;
00409 }
00410
00411 }
00412
00413 return 0;
00414
00415
00416
00417 }