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


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