zgsvts.c
Go to the documentation of this file.
00001 /* zgsvts.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
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 /* Subroutine */ 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     /* System generated locals */
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     /* Local variables */
00040     integer i__, j, k, l;
00041     doublereal ulp;
00042     integer info;
00043     doublereal unfl, temp, resid, anorm, bnorm;
00044     extern /* Subroutine */ 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 /* Subroutine */ 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 /*  -- LAPACK test routine (version 3.1) -- */
00068 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00069 /*     November 2006 */
00070 
00071 /*     .. Scalar Arguments .. */
00072 /*     .. */
00073 /*     .. Array Arguments .. */
00074 /*     .. */
00075 
00076 /*  Purpose */
00077 /*  ======= */
00078 
00079 /*  ZGSVTS tests ZGGSVD, which computes the GSVD of an M-by-N matrix A */
00080 /*  and a P-by-N matrix B: */
00081 /*               U'*A*Q = D1*R and V'*B*Q = D2*R. */
00082 
00083 /*  Arguments */
00084 /*  ========= */
00085 
00086 /*  M       (input) INTEGER */
00087 /*          The number of rows of the matrix A.  M >= 0. */
00088 
00089 /*  P       (input) INTEGER */
00090 /*          The number of rows of the matrix B.  P >= 0. */
00091 
00092 /*  N       (input) INTEGER */
00093 /*          The number of columns of the matrices A and B.  N >= 0. */
00094 
00095 /*  A       (input) COMPLEX*16 array, dimension (LDA,M) */
00096 /*          The M-by-N matrix A. */
00097 
00098 /*  AF      (output) COMPLEX*16 array, dimension (LDA,N) */
00099 /*          Details of the GSVD of A and B, as returned by ZGGSVD, */
00100 /*          see ZGGSVD for further details. */
00101 
00102 /*  LDA     (input) INTEGER */
00103 /*          The leading dimension of the arrays A and AF. */
00104 /*          LDA >= max( 1,M ). */
00105 
00106 /*  B       (input) COMPLEX*16 array, dimension (LDB,P) */
00107 /*          On entry, the P-by-N matrix B. */
00108 
00109 /*  BF      (output) COMPLEX*16 array, dimension (LDB,N) */
00110 /*          Details of the GSVD of A and B, as returned by ZGGSVD, */
00111 /*          see ZGGSVD for further details. */
00112 
00113 /*  LDB     (input) INTEGER */
00114 /*          The leading dimension of the arrays B and BF. */
00115 /*          LDB >= max(1,P). */
00116 
00117 /*  U       (output) COMPLEX*16 array, dimension(LDU,M) */
00118 /*          The M by M unitary matrix U. */
00119 
00120 /*  LDU     (input) INTEGER */
00121 /*          The leading dimension of the array U. LDU >= max(1,M). */
00122 
00123 /*  V       (output) COMPLEX*16 array, dimension(LDV,M) */
00124 /*          The P by P unitary matrix V. */
00125 
00126 /*  LDV     (input) INTEGER */
00127 /*          The leading dimension of the array V. LDV >= max(1,P). */
00128 
00129 /*  Q       (output) COMPLEX*16 array, dimension(LDQ,N) */
00130 /*          The N by N unitary matrix Q. */
00131 
00132 /*  LDQ     (input) INTEGER */
00133 /*          The leading dimension of the array Q. LDQ >= max(1,N). */
00134 
00135 /*  ALPHA   (output) DOUBLE PRECISION array, dimension (N) */
00136 /*  BETA    (output) DOUBLE PRECISION array, dimension (N) */
00137 /*          The generalized singular value pairs of A and B, the */
00138 /*          ``diagonal'' matrices D1 and D2 are constructed from */
00139 /*          ALPHA and BETA, see subroutine ZGGSVD for details. */
00140 
00141 /*  R       (output) COMPLEX*16 array, dimension(LDQ,N) */
00142 /*          The upper triangular matrix R. */
00143 
00144 /*  LDR     (input) INTEGER */
00145 /*          The leading dimension of the array R. LDR >= max(1,N). */
00146 
00147 /*  IWORK   (workspace) INTEGER array, dimension (N) */
00148 
00149 /*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) */
00150 
00151 /*  LWORK   (input) INTEGER */
00152 /*          The dimension of the array WORK, */
00153 /*          LWORK >= max(M,P,N)*max(M,P,N). */
00154 
00155 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(M,P,N)) */
00156 
00157 /*  RESULT  (output) DOUBLE PRECISION array, dimension (5) */
00158 /*          The test ratios: */
00159 /*          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP) */
00160 /*          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP) */
00161 /*          RESULT(3) = norm( I - U'*U ) / ( M*ULP ) */
00162 /*          RESULT(4) = norm( I - V'*V ) / ( P*ULP ) */
00163 /*          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP ) */
00164 /*          RESULT(6) = 0        if ALPHA is in decreasing order; */
00165 /*                    = ULPINV   otherwise. */
00166 
00167 /*  ===================================================================== */
00168 
00169 /*     .. Parameters .. */
00170 /*     .. */
00171 /*     .. Local Scalars .. */
00172 /*     .. */
00173 /*     .. External Functions .. */
00174 /*     .. */
00175 /*     .. External Subroutines .. */
00176 /*     .. */
00177 /*     .. Intrinsic Functions .. */
00178 /*     .. */
00179 /*     .. Executable Statements .. */
00180 
00181     /* Parameter adjustments */
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     /* Function Body */
00214     ulp = dlamch_("Precision");
00215     ulpinv = 1. / ulp;
00216     unfl = dlamch_("Safe minimum");
00217 
00218 /*     Copy the matrix A to the array AF. */
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 /* Computing MAX */
00224     d__1 = zlange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00225     anorm = max(d__1,unfl);
00226 /* Computing MAX */
00227     d__1 = zlange_("1", p, n, &b[b_offset], ldb, &rwork[1]);
00228     bnorm = max(d__1,unfl);
00229 
00230 /*     Factorize the matrices A and B in the arrays AF and BF. */
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 /*     Copy R */
00238 
00239 /* Computing MIN */
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 /* L10: */
00249         }
00250 /* L20: */
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 /* L30: */
00262             }
00263 /* L40: */
00264         }
00265     }
00266 
00267 /*     Compute A:= U'*A*Q - D1*R */
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 /* L50: */
00286         }
00287 /* L60: */
00288     }
00289 
00290 /* Computing MIN */
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 /* L70: */
00305         }
00306 /* L80: */
00307     }
00308 
00309 /*     Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) . */
00310 
00311     resid = zlange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00312     if (anorm > 0.) {
00313 /* Computing MAX */
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 /*     Compute B := V'*B*Q - D2*R */
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 /* L90: */
00341         }
00342 /* L100: */
00343     }
00344 
00345 /*     Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) . */
00346 
00347     resid = zlange_("1", p, n, &b[b_offset], ldb, &rwork[1]);
00348     if (bnorm > 0.) {
00349 /* Computing MAX */
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 /*     Compute I - U'*U */
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 /*     Compute norm( I - U'*U ) / ( M * ULP ) . */
00363 
00364     resid = zlanhe_("1", "Upper", m, &work[1], ldu, &rwork[1]);
00365     result[3] = resid / (doublereal) max(1,*m) / ulp;
00366 
00367 /*     Compute I - V'*V */
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 /*     Compute norm( I - V'*V ) / ( P * ULP ) . */
00374 
00375     resid = zlanhe_("1", "Upper", p, &work[1], ldv, &rwork[1]);
00376     result[4] = resid / (doublereal) max(1,*p) / ulp;
00377 
00378 /*     Compute I - Q'*Q */
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 /*     Compute norm( I - Q'*Q ) / ( N * ULP ) . */
00385 
00386     resid = zlanhe_("1", "Upper", n, &work[1], ldq, &rwork[1]);
00387     result[5] = resid / (doublereal) max(1,*n) / ulp;
00388 
00389 /*     Check sorting */
00390 
00391     dcopy_(n, &alpha[1], &c__1, &rwork[1], &c__1);
00392 /* Computing MIN */
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 /* L110: */
00403     }
00404 
00405     result[6] = 0.;
00406 /* Computing MIN */
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 /* L120: */
00414     }
00415 
00416     return 0;
00417 
00418 /*     End of ZGSVTS */
00419 
00420 } /* zgsvts_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:36