cunt03.c
Go to the documentation of this file.
00001 /* cunt03.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 cunt03_(char *rc, integer *mu, integer *mv, integer *n, 
00021         integer *k, complex *u, integer *ldu, complex *v, integer *ldv, 
00022         complex *work, integer *lwork, real *rwork, real *result, integer *
00023         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     real r__1, r__2;
00028     complex q__1, q__2;
00029 
00030     /* Builtin functions */
00031     double c_abs(complex *);
00032     void c_div(complex *, complex *, complex *);
00033 
00034     /* Local variables */
00035     integer i__, j;
00036     complex s, su, sv;
00037     integer irc, lmx;
00038     real ulp, res1, res2;
00039     extern logical lsame_(char *, char *);
00040     extern /* Subroutine */ int cunt01_(char *, integer *, integer *, complex 
00041             *, integer *, complex *, integer *, real *, real *);
00042     extern integer icamax_(integer *, complex *, integer *);
00043     extern doublereal slamch_(char *);
00044     extern /* Subroutine */ int xerbla_(char *, integer *);
00045 
00046 
00047 /*  -- LAPACK test routine (version 3.1) -- */
00048 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00049 /*     November 2006 */
00050 
00051 /*     .. Scalar Arguments .. */
00052 /*     .. */
00053 /*     .. Array Arguments .. */
00054 /*     .. */
00055 
00056 /*  Purpose */
00057 /*  ======= */
00058 
00059 /*  CUNT03 compares two unitary matrices U and V to see if their */
00060 /*  corresponding rows or columns span the same spaces.  The rows are */
00061 /*  checked if RC = 'R', and the columns are checked if RC = 'C'. */
00062 
00063 /*  RESULT is the maximum of */
00064 
00065 /*     | V*V' - I | / ( MV ulp ), if RC = 'R', or */
00066 
00067 /*     | V'*V - I | / ( MV ulp ), if RC = 'C', */
00068 
00069 /*  and the maximum over rows (or columns) 1 to K of */
00070 
00071 /*     | U(i) - S*V(i) |/ ( N ulp ) */
00072 
00073 /*  where abs(S) = 1 (chosen to minimize the expression), U(i) is the */
00074 /*  i-th row (column) of U, and V(i) is the i-th row (column) of V. */
00075 
00076 /*  Arguments */
00077 /*  ========== */
00078 
00079 /*  RC      (input) CHARACTER*1 */
00080 /*          If RC = 'R' the rows of U and V are to be compared. */
00081 /*          If RC = 'C' the columns of U and V are to be compared. */
00082 
00083 /*  MU      (input) INTEGER */
00084 /*          The number of rows of U if RC = 'R', and the number of */
00085 /*          columns if RC = 'C'.  If MU = 0 CUNT03 does nothing. */
00086 /*          MU must be at least zero. */
00087 
00088 /*  MV      (input) INTEGER */
00089 /*          The number of rows of V if RC = 'R', and the number of */
00090 /*          columns if RC = 'C'.  If MV = 0 CUNT03 does nothing. */
00091 /*          MV must be at least zero. */
00092 
00093 /*  N       (input) INTEGER */
00094 /*          If RC = 'R', the number of columns in the matrices U and V, */
00095 /*          and if RC = 'C', the number of rows in U and V.  If N = 0 */
00096 /*          CUNT03 does nothing.  N must be at least zero. */
00097 
00098 /*  K       (input) INTEGER */
00099 /*          The number of rows or columns of U and V to compare. */
00100 /*          0 <= K <= max(MU,MV). */
00101 
00102 /*  U       (input) COMPLEX array, dimension (LDU,N) */
00103 /*          The first matrix to compare.  If RC = 'R', U is MU by N, and */
00104 /*          if RC = 'C', U is N by MU. */
00105 
00106 /*  LDU     (input) INTEGER */
00107 /*          The leading dimension of U.  If RC = 'R', LDU >= max(1,MU), */
00108 /*          and if RC = 'C', LDU >= max(1,N). */
00109 
00110 /*  V       (input) COMPLEX array, dimension (LDV,N) */
00111 /*          The second matrix to compare.  If RC = 'R', V is MV by N, and */
00112 /*          if RC = 'C', V is N by MV. */
00113 
00114 /*  LDV     (input) INTEGER */
00115 /*          The leading dimension of V.  If RC = 'R', LDV >= max(1,MV), */
00116 /*          and if RC = 'C', LDV >= max(1,N). */
00117 
00118 /*  WORK    (workspace) COMPLEX array, dimension (LWORK) */
00119 
00120 /*  LWORK   (input) INTEGER */
00121 /*          The length of the array WORK.  For best performance, LWORK */
00122 /*          should be at least N*N if RC = 'C' or M*M if RC = 'R', but */
00123 /*          the tests will be done even if LWORK is 0. */
00124 
00125 /*  RWORK   (workspace) REAL array, dimension (max(MV,N)) */
00126 
00127 /*  RESULT  (output) REAL */
00128 /*          The value computed by the test described above.  RESULT is */
00129 /*          limited to 1/ulp to avoid overflow. */
00130 
00131 /*  INFO    (output) INTEGER */
00132 /*          0  indicates a successful exit */
00133 /*          -k indicates the k-th parameter had an illegal value */
00134 
00135 /*  ===================================================================== */
00136 
00137 
00138 /*     .. Parameters .. */
00139 /*     .. */
00140 /*     .. Local Scalars .. */
00141 /*     .. */
00142 /*     .. External Functions .. */
00143 /*     .. */
00144 /*     .. Intrinsic Functions .. */
00145 /*     .. */
00146 /*     .. External Subroutines .. */
00147 /*     .. */
00148 /*     .. Executable Statements .. */
00149 
00150 /*     Check inputs */
00151 
00152     /* Parameter adjustments */
00153     u_dim1 = *ldu;
00154     u_offset = 1 + u_dim1;
00155     u -= u_offset;
00156     v_dim1 = *ldv;
00157     v_offset = 1 + v_dim1;
00158     v -= v_offset;
00159     --work;
00160     --rwork;
00161 
00162     /* Function Body */
00163     *info = 0;
00164     if (lsame_(rc, "R")) {
00165         irc = 0;
00166     } else if (lsame_(rc, "C")) {
00167         irc = 1;
00168     } else {
00169         irc = -1;
00170     }
00171     if (irc == -1) {
00172         *info = -1;
00173     } else if (*mu < 0) {
00174         *info = -2;
00175     } else if (*mv < 0) {
00176         *info = -3;
00177     } else if (*n < 0) {
00178         *info = -4;
00179     } else if (*k < 0 || *k > max(*mu,*mv)) {
00180         *info = -5;
00181     } else if (irc == 0 && *ldu < max(1,*mu) || irc == 1 && *ldu < max(1,*n)) 
00182             {
00183         *info = -7;
00184     } else if (irc == 0 && *ldv < max(1,*mv) || irc == 1 && *ldv < max(1,*n)) 
00185             {
00186         *info = -9;
00187     }
00188     if (*info != 0) {
00189         i__1 = -(*info);
00190         xerbla_("CUNT03", &i__1);
00191         return 0;
00192     }
00193 
00194 /*     Initialize result */
00195 
00196     *result = 0.f;
00197     if (*mu == 0 || *mv == 0 || *n == 0) {
00198         return 0;
00199     }
00200 
00201 /*     Machine constants */
00202 
00203     ulp = slamch_("Precision");
00204 
00205     if (irc == 0) {
00206 
00207 /*        Compare rows */
00208 
00209         res1 = 0.f;
00210         i__1 = *k;
00211         for (i__ = 1; i__ <= i__1; ++i__) {
00212             lmx = icamax_(n, &u[i__ + u_dim1], ldu);
00213             i__2 = i__ + lmx * v_dim1;
00214             if (v[i__2].r == 0.f && v[i__2].i == 0.f) {
00215                 sv.r = 1.f, sv.i = 0.f;
00216             } else {
00217                 r__1 = c_abs(&v[i__ + lmx * v_dim1]);
00218                 q__2.r = r__1, q__2.i = 0.f;
00219                 c_div(&q__1, &q__2, &v[i__ + lmx * v_dim1]);
00220                 sv.r = q__1.r, sv.i = q__1.i;
00221             }
00222             i__2 = i__ + lmx * u_dim1;
00223             if (u[i__2].r == 0.f && u[i__2].i == 0.f) {
00224                 su.r = 1.f, su.i = 0.f;
00225             } else {
00226                 r__1 = c_abs(&u[i__ + lmx * u_dim1]);
00227                 q__2.r = r__1, q__2.i = 0.f;
00228                 c_div(&q__1, &q__2, &u[i__ + lmx * u_dim1]);
00229                 su.r = q__1.r, su.i = q__1.i;
00230             }
00231             c_div(&q__1, &sv, &su);
00232             s.r = q__1.r, s.i = q__1.i;
00233             i__2 = *n;
00234             for (j = 1; j <= i__2; ++j) {
00235 /* Computing MAX */
00236                 i__3 = i__ + j * u_dim1;
00237                 i__4 = i__ + j * v_dim1;
00238                 q__2.r = s.r * v[i__4].r - s.i * v[i__4].i, q__2.i = s.r * v[
00239                         i__4].i + s.i * v[i__4].r;
00240                 q__1.r = u[i__3].r - q__2.r, q__1.i = u[i__3].i - q__2.i;
00241                 r__1 = res1, r__2 = c_abs(&q__1);
00242                 res1 = dmax(r__1,r__2);
00243 /* L10: */
00244             }
00245 /* L20: */
00246         }
00247         res1 /= (real) (*n) * ulp;
00248 
00249 /*        Compute orthogonality of rows of V. */
00250 
00251         cunt01_("Rows", mv, n, &v[v_offset], ldv, &work[1], lwork, &rwork[1], 
00252                 &res2);
00253 
00254     } else {
00255 
00256 /*        Compare columns */
00257 
00258         res1 = 0.f;
00259         i__1 = *k;
00260         for (i__ = 1; i__ <= i__1; ++i__) {
00261             lmx = icamax_(n, &u[i__ * u_dim1 + 1], &c__1);
00262             i__2 = lmx + i__ * v_dim1;
00263             if (v[i__2].r == 0.f && v[i__2].i == 0.f) {
00264                 sv.r = 1.f, sv.i = 0.f;
00265             } else {
00266                 r__1 = c_abs(&v[lmx + i__ * v_dim1]);
00267                 q__2.r = r__1, q__2.i = 0.f;
00268                 c_div(&q__1, &q__2, &v[lmx + i__ * v_dim1]);
00269                 sv.r = q__1.r, sv.i = q__1.i;
00270             }
00271             i__2 = lmx + i__ * u_dim1;
00272             if (u[i__2].r == 0.f && u[i__2].i == 0.f) {
00273                 su.r = 1.f, su.i = 0.f;
00274             } else {
00275                 r__1 = c_abs(&u[lmx + i__ * u_dim1]);
00276                 q__2.r = r__1, q__2.i = 0.f;
00277                 c_div(&q__1, &q__2, &u[lmx + i__ * u_dim1]);
00278                 su.r = q__1.r, su.i = q__1.i;
00279             }
00280             c_div(&q__1, &sv, &su);
00281             s.r = q__1.r, s.i = q__1.i;
00282             i__2 = *n;
00283             for (j = 1; j <= i__2; ++j) {
00284 /* Computing MAX */
00285                 i__3 = j + i__ * u_dim1;
00286                 i__4 = j + i__ * v_dim1;
00287                 q__2.r = s.r * v[i__4].r - s.i * v[i__4].i, q__2.i = s.r * v[
00288                         i__4].i + s.i * v[i__4].r;
00289                 q__1.r = u[i__3].r - q__2.r, q__1.i = u[i__3].i - q__2.i;
00290                 r__1 = res1, r__2 = c_abs(&q__1);
00291                 res1 = dmax(r__1,r__2);
00292 /* L30: */
00293             }
00294 /* L40: */
00295         }
00296         res1 /= (real) (*n) * ulp;
00297 
00298 /*        Compute orthogonality of columns of V. */
00299 
00300         cunt01_("Columns", n, mv, &v[v_offset], ldv, &work[1], lwork, &rwork[
00301                 1], &res2);
00302     }
00303 
00304 /* Computing MIN */
00305     r__1 = dmax(res1,res2), r__2 = 1.f / ulp;
00306     *result = dmin(r__1,r__2);
00307     return 0;
00308 
00309 /*     End of CUNT03 */
00310 
00311 } /* cunt03_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:35