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