cunt01.c
Go to the documentation of this file.
00001 /* cunt01.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 complex c_b7 = {0.f,0.f};
00019 static complex c_b8 = {1.f,0.f};
00020 static real c_b10 = -1.f;
00021 static real c_b11 = 1.f;
00022 static integer c__1 = 1;
00023 
00024 /* Subroutine */ int cunt01_(char *rowcol, integer *m, integer *n, complex *u, 
00025          integer *ldu, complex *work, integer *lwork, real *rwork, real *
00026         resid)
00027 {
00028     /* System generated locals */
00029     integer u_dim1, u_offset, i__1, i__2;
00030     real r__1, r__2, r__3, r__4;
00031     complex q__1, q__2;
00032 
00033     /* Builtin functions */
00034     double r_imag(complex *);
00035 
00036     /* Local variables */
00037     integer i__, j, k;
00038     real eps;
00039     complex tmp;
00040     extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer 
00041             *, complex *, integer *);
00042     extern /* Subroutine */ int cherk_(char *, char *, integer *, integer *, 
00043             real *, complex *, integer *, real *, complex *, integer *);
00044     extern logical lsame_(char *, char *);
00045     integer mnmin;
00046     extern doublereal slamch_(char *);
00047     extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
00048             *, complex *, complex *, integer *);
00049     extern doublereal clansy_(char *, char *, integer *, complex *, integer *, 
00050              real *);
00051     integer ldwork;
00052     char transu[1];
00053 
00054 
00055 /*  -- LAPACK test routine (version 3.1) -- */
00056 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00057 /*     November 2006 */
00058 
00059 /*     .. Scalar Arguments .. */
00060 /*     .. */
00061 /*     .. Array Arguments .. */
00062 /*     .. */
00063 
00064 /*  Purpose */
00065 /*  ======= */
00066 
00067 /*  CUNT01 checks that the matrix U is unitary by computing the ratio */
00068 
00069 /*     RESID = norm( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', */
00070 /*  or */
00071 /*     RESID = norm( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. */
00072 
00073 /*  Alternatively, if there isn't sufficient workspace to form */
00074 /*  I - U*U' or I - U'*U, the ratio is computed as */
00075 
00076 /*     RESID = abs( I - U*U' ) / ( n * EPS ), if ROWCOL = 'R', */
00077 /*  or */
00078 /*     RESID = abs( I - U'*U ) / ( m * EPS ), if ROWCOL = 'C'. */
00079 
00080 /*  where EPS is the machine precision.  ROWCOL is used only if m = n; */
00081 /*  if m > n, ROWCOL is assumed to be 'C', and if m < n, ROWCOL is */
00082 /*  assumed to be 'R'. */
00083 
00084 /*  Arguments */
00085 /*  ========= */
00086 
00087 /*  ROWCOL  (input) CHARACTER */
00088 /*          Specifies whether the rows or columns of U should be checked */
00089 /*          for orthogonality.  Used only if M = N. */
00090 /*          = 'R':  Check for orthogonal rows of U */
00091 /*          = 'C':  Check for orthogonal columns of U */
00092 
00093 /*  M       (input) INTEGER */
00094 /*          The number of rows of the matrix U. */
00095 
00096 /*  N       (input) INTEGER */
00097 /*          The number of columns of the matrix U. */
00098 
00099 /*  U       (input) COMPLEX array, dimension (LDU,N) */
00100 /*          The unitary matrix U.  U is checked for orthogonal columns */
00101 /*          if m > n or if m = n and ROWCOL = 'C'.  U is checked for */
00102 /*          orthogonal rows if m < n or if m = n and ROWCOL = 'R'. */
00103 
00104 /*  LDU     (input) INTEGER */
00105 /*          The leading dimension of the array U.  LDU >= max(1,M). */
00106 
00107 /*  WORK    (workspace) COMPLEX array, dimension (LWORK) */
00108 
00109 /*  LWORK   (input) INTEGER */
00110 /*          The length of the array WORK.  For best performance, LWORK */
00111 /*          should be at least N*N if ROWCOL = 'C' or M*M if */
00112 /*          ROWCOL = 'R', but the test will be done even if LWORK is 0. */
00113 
00114 /*  RWORK   (workspace) REAL array, dimension (min(M,N)) */
00115 /*          Used only if LWORK is large enough to use the Level 3 BLAS */
00116 /*          code. */
00117 
00118 /*  RESID   (output) REAL */
00119 /*          RESID = norm( I - U * U' ) / ( n * EPS ), if ROWCOL = 'R', or */
00120 /*          RESID = norm( I - U' * U ) / ( m * EPS ), if ROWCOL = 'C'. */
00121 
00122 /*  ===================================================================== */
00123 
00124 /*     .. Parameters .. */
00125 /*     .. */
00126 /*     .. Local Scalars .. */
00127 /*     .. */
00128 /*     .. External Functions .. */
00129 /*     .. */
00130 /*     .. External Subroutines .. */
00131 /*     .. */
00132 /*     .. Intrinsic Functions .. */
00133 /*     .. */
00134 /*     .. Statement Functions .. */
00135 /*     .. */
00136 /*     .. Statement Function definitions .. */
00137 /*     .. */
00138 /*     .. Executable Statements .. */
00139 
00140     /* Parameter adjustments */
00141     u_dim1 = *ldu;
00142     u_offset = 1 + u_dim1;
00143     u -= u_offset;
00144     --work;
00145     --rwork;
00146 
00147     /* Function Body */
00148     *resid = 0.f;
00149 
00150 /*     Quick return if possible */
00151 
00152     if (*m <= 0 || *n <= 0) {
00153         return 0;
00154     }
00155 
00156     eps = slamch_("Precision");
00157     if (*m < *n || *m == *n && lsame_(rowcol, "R")) {
00158         *(unsigned char *)transu = 'N';
00159         k = *n;
00160     } else {
00161         *(unsigned char *)transu = 'C';
00162         k = *m;
00163     }
00164     mnmin = min(*m,*n);
00165 
00166     if ((mnmin + 1) * mnmin <= *lwork) {
00167         ldwork = mnmin;
00168     } else {
00169         ldwork = 0;
00170     }
00171     if (ldwork > 0) {
00172 
00173 /*        Compute I - U*U' or I - U'*U. */
00174 
00175         claset_("Upper", &mnmin, &mnmin, &c_b7, &c_b8, &work[1], &ldwork);
00176         cherk_("Upper", transu, &mnmin, &k, &c_b10, &u[u_offset], ldu, &c_b11, 
00177                  &work[1], &ldwork);
00178 
00179 /*        Compute norm( I - U*U' ) / ( K * EPS ) . */
00180 
00181         *resid = clansy_("1", "Upper", &mnmin, &work[1], &ldwork, &rwork[1]);
00182         *resid = *resid / (real) k / eps;
00183     } else if (*(unsigned char *)transu == 'C') {
00184 
00185 /*        Find the maximum element in abs( I - U'*U ) / ( m * EPS ) */
00186 
00187         i__1 = *n;
00188         for (j = 1; j <= i__1; ++j) {
00189             i__2 = j;
00190             for (i__ = 1; i__ <= i__2; ++i__) {
00191                 if (i__ != j) {
00192                     tmp.r = 0.f, tmp.i = 0.f;
00193                 } else {
00194                     tmp.r = 1.f, tmp.i = 0.f;
00195                 }
00196                 cdotc_(&q__2, m, &u[i__ * u_dim1 + 1], &c__1, &u[j * u_dim1 + 
00197                         1], &c__1);
00198                 q__1.r = tmp.r - q__2.r, q__1.i = tmp.i - q__2.i;
00199                 tmp.r = q__1.r, tmp.i = q__1.i;
00200 /* Computing MAX */
00201                 r__3 = *resid, r__4 = (r__1 = tmp.r, dabs(r__1)) + (r__2 = 
00202                         r_imag(&tmp), dabs(r__2));
00203                 *resid = dmax(r__3,r__4);
00204 /* L10: */
00205             }
00206 /* L20: */
00207         }
00208         *resid = *resid / (real) (*m) / eps;
00209     } else {
00210 
00211 /*        Find the maximum element in abs( I - U*U' ) / ( n * EPS ) */
00212 
00213         i__1 = *m;
00214         for (j = 1; j <= i__1; ++j) {
00215             i__2 = j;
00216             for (i__ = 1; i__ <= i__2; ++i__) {
00217                 if (i__ != j) {
00218                     tmp.r = 0.f, tmp.i = 0.f;
00219                 } else {
00220                     tmp.r = 1.f, tmp.i = 0.f;
00221                 }
00222                 cdotc_(&q__2, n, &u[j + u_dim1], ldu, &u[i__ + u_dim1], ldu);
00223                 q__1.r = tmp.r - q__2.r, q__1.i = tmp.i - q__2.i;
00224                 tmp.r = q__1.r, tmp.i = q__1.i;
00225 /* Computing MAX */
00226                 r__3 = *resid, r__4 = (r__1 = tmp.r, dabs(r__1)) + (r__2 = 
00227                         r_imag(&tmp), dabs(r__2));
00228                 *resid = dmax(r__3,r__4);
00229 /* L30: */
00230             }
00231 /* L40: */
00232         }
00233         *resid = *resid / (real) (*n) / eps;
00234     }
00235     return 0;
00236 
00237 /*     End of CUNT01 */
00238 
00239 } /* cunt01_ */


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