zunt03.c
Go to the documentation of this file.
00001 /* zunt03.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 integer c__1 = 1;
00019 
00020 /* Subroutine */ 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     /* System generated locals */
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     /* Builtin functions */
00031     double z_abs(doublecomplex *);
00032     void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00033 
00034     /* Local variables */
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 /* Subroutine */ int zunt01_(char *, integer *, integer *, 
00041             doublecomplex *, integer *, doublecomplex *, integer *, 
00042             doublereal *, doublereal *);
00043     extern doublereal dlamch_(char *);
00044     extern /* Subroutine */ int xerbla_(char *, integer *);
00045     extern integer izamax_(integer *, doublecomplex *, integer *);
00046 
00047 
00048 /*  -- LAPACK test routine (version 3.1) -- */
00049 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00050 /*     November 2006 */
00051 
00052 /*     .. Scalar Arguments .. */
00053 /*     .. */
00054 /*     .. Array Arguments .. */
00055 /*     .. */
00056 
00057 /*  Purpose */
00058 /*  ======= */
00059 
00060 /*  ZUNT03 compares two unitary matrices U and V to see if their */
00061 /*  corresponding rows or columns span the same spaces.  The rows are */
00062 /*  checked if RC = 'R', and the columns are checked if RC = 'C'. */
00063 
00064 /*  RESULT is the maximum of */
00065 
00066 /*     | V*V' - I | / ( MV ulp ), if RC = 'R', or */
00067 
00068 /*     | V'*V - I | / ( MV ulp ), if RC = 'C', */
00069 
00070 /*  and the maximum over rows (or columns) 1 to K of */
00071 
00072 /*     | U(i) - S*V(i) |/ ( N ulp ) */
00073 
00074 /*  where abs(S) = 1 (chosen to minimize the expression), U(i) is the */
00075 /*  i-th row (column) of U, and V(i) is the i-th row (column) of V. */
00076 
00077 /*  Arguments */
00078 /*  ========== */
00079 
00080 /*  RC      (input) CHARACTER*1 */
00081 /*          If RC = 'R' the rows of U and V are to be compared. */
00082 /*          If RC = 'C' the columns of U and V are to be compared. */
00083 
00084 /*  MU      (input) INTEGER */
00085 /*          The number of rows of U if RC = 'R', and the number of */
00086 /*          columns if RC = 'C'.  If MU = 0 ZUNT03 does nothing. */
00087 /*          MU must be at least zero. */
00088 
00089 /*  MV      (input) INTEGER */
00090 /*          The number of rows of V if RC = 'R', and the number of */
00091 /*          columns if RC = 'C'.  If MV = 0 ZUNT03 does nothing. */
00092 /*          MV must be at least zero. */
00093 
00094 /*  N       (input) INTEGER */
00095 /*          If RC = 'R', the number of columns in the matrices U and V, */
00096 /*          and if RC = 'C', the number of rows in U and V.  If N = 0 */
00097 /*          ZUNT03 does nothing.  N must be at least zero. */
00098 
00099 /*  K       (input) INTEGER */
00100 /*          The number of rows or columns of U and V to compare. */
00101 /*          0 <= K <= max(MU,MV). */
00102 
00103 /*  U       (input) COMPLEX*16 array, dimension (LDU,N) */
00104 /*          The first matrix to compare.  If RC = 'R', U is MU by N, and */
00105 /*          if RC = 'C', U is N by MU. */
00106 
00107 /*  LDU     (input) INTEGER */
00108 /*          The leading dimension of U.  If RC = 'R', LDU >= max(1,MU), */
00109 /*          and if RC = 'C', LDU >= max(1,N). */
00110 
00111 /*  V       (input) COMPLEX*16 array, dimension (LDV,N) */
00112 /*          The second matrix to compare.  If RC = 'R', V is MV by N, and */
00113 /*          if RC = 'C', V is N by MV. */
00114 
00115 /*  LDV     (input) INTEGER */
00116 /*          The leading dimension of V.  If RC = 'R', LDV >= max(1,MV), */
00117 /*          and if RC = 'C', LDV >= max(1,N). */
00118 
00119 /*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) */
00120 
00121 /*  LWORK   (input) INTEGER */
00122 /*          The length of the array WORK.  For best performance, LWORK */
00123 /*          should be at least N*N if RC = 'C' or M*M if RC = 'R', but */
00124 /*          the tests will be done even if LWORK is 0. */
00125 
00126 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(MV,N)) */
00127 
00128 /*  RESULT  (output) DOUBLE PRECISION */
00129 /*          The value computed by the test described above.  RESULT is */
00130 /*          limited to 1/ulp to avoid overflow. */
00131 
00132 /*  INFO    (output) INTEGER */
00133 /*          0  indicates a successful exit */
00134 /*          -k indicates the k-th parameter had an illegal value */
00135 
00136 /*  ===================================================================== */
00137 
00138 
00139 /*     .. Parameters .. */
00140 /*     .. */
00141 /*     .. Local Scalars .. */
00142 /*     .. */
00143 /*     .. External Functions .. */
00144 /*     .. */
00145 /*     .. Intrinsic Functions .. */
00146 /*     .. */
00147 /*     .. External Subroutines .. */
00148 /*     .. */
00149 /*     .. Executable Statements .. */
00150 
00151 /*     Check inputs */
00152 
00153     /* Parameter adjustments */
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     /* Function Body */
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 /*     Initialize result */
00196 
00197     *result = 0.;
00198     if (*mu == 0 || *mv == 0 || *n == 0) {
00199         return 0;
00200     }
00201 
00202 /*     Machine constants */
00203 
00204     ulp = dlamch_("Precision");
00205 
00206     if (irc == 0) {
00207 
00208 /*        Compare rows */
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 /* Computing MAX */
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 /* L10: */
00245             }
00246 /* L20: */
00247         }
00248         res1 /= (doublereal) (*n) * ulp;
00249 
00250 /*        Compute orthogonality of rows of V. */
00251 
00252         zunt01_("Rows", mv, n, &v[v_offset], ldv, &work[1], lwork, &rwork[1], 
00253                 &res2);
00254 
00255     } else {
00256 
00257 /*        Compare columns */
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 /* Computing MAX */
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 /* L30: */
00294             }
00295 /* L40: */
00296         }
00297         res1 /= (doublereal) (*n) * ulp;
00298 
00299 /*        Compute orthogonality of columns of V. */
00300 
00301         zunt01_("Columns", n, mv, &v[v_offset], ldv, &work[1], lwork, &rwork[
00302                 1], &res2);
00303     }
00304 
00305 /* Computing MIN */
00306     d__1 = max(res1,res2), d__2 = 1. / ulp;
00307     *result = min(d__1,d__2);
00308     return 0;
00309 
00310 /*     End of ZUNT03 */
00311 
00312 } /* zunt03_ */

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