cget54.c
Go to the documentation of this file.
00001 /* cget54.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_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020 
00021 /* Subroutine */ int cget54_(integer *n, complex *a, integer *lda, complex *b, 
00022          integer *ldb, complex *s, integer *lds, complex *t, integer *ldt, 
00023         complex *u, integer *ldu, complex *v, integer *ldv, complex *work, 
00024         real *result)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, b_dim1, b_offset, s_dim1, s_offset, t_dim1, 
00028             t_offset, u_dim1, u_offset, v_dim1, v_offset, i__1;
00029     real r__1, r__2;
00030     complex q__1;
00031 
00032     /* Local variables */
00033     real dum[1], ulp, unfl;
00034     extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, 
00035             integer *, complex *, complex *, integer *, complex *, integer *, 
00036             complex *, complex *, integer *);
00037     real wnorm;
00038     extern doublereal clange_(char *, integer *, integer *, complex *, 
00039             integer *, real *), slamch_(char *);
00040     extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
00041             *, integer *, complex *, integer *);
00042     real abnorm;
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 /*  CGET54 checks a generalized decomposition of the form */
00058 
00059 /*           A = U*S*V'  and B = U*T* V' */
00060 
00061 /*  where ' means conjugate transpose and U and V are unitary. */
00062 
00063 /*  Specifically, */
00064 
00065 /*    RESULT = ||( A - U*S*V', B - U*T*V' )|| / (||( A, B )||*n*ulp ) */
00066 
00067 /*  Arguments */
00068 /*  ========= */
00069 
00070 /*  N       (input) INTEGER */
00071 /*          The size of the matrix.  If it is zero, SGET54 does nothing. */
00072 /*          It must be at least zero. */
00073 
00074 /*  A       (input) COMPLEX array, dimension (LDA, N) */
00075 /*          The original (unfactored) matrix A. */
00076 
00077 /*  LDA     (input) INTEGER */
00078 /*          The leading dimension of A.  It must be at least 1 */
00079 /*          and at least N. */
00080 
00081 /*  B       (input) COMPLEX array, dimension (LDB, N) */
00082 /*          The original (unfactored) matrix B. */
00083 
00084 /*  LDB     (input) INTEGER */
00085 /*          The leading dimension of B.  It must be at least 1 */
00086 /*          and at least N. */
00087 
00088 /*  S       (input) COMPLEX array, dimension (LDS, N) */
00089 /*          The factored matrix S. */
00090 
00091 /*  LDS     (input) INTEGER */
00092 /*          The leading dimension of S.  It must be at least 1 */
00093 /*          and at least N. */
00094 
00095 /*  T       (input) COMPLEX array, dimension (LDT, N) */
00096 /*          The factored matrix T. */
00097 
00098 /*  LDT     (input) INTEGER */
00099 /*          The leading dimension of T.  It must be at least 1 */
00100 /*          and at least N. */
00101 
00102 /*  U       (input) COMPLEX array, dimension (LDU, N) */
00103 /*          The orthogonal matrix on the left-hand side in the */
00104 /*          decomposition. */
00105 
00106 /*  LDU     (input) INTEGER */
00107 /*          The leading dimension of U.  LDU must be at least N and */
00108 /*          at least 1. */
00109 
00110 /*  V       (input) COMPLEX array, dimension (LDV, N) */
00111 /*          The orthogonal matrix on the left-hand side in the */
00112 /*          decomposition. */
00113 
00114 /*  LDV     (input) INTEGER */
00115 /*          The leading dimension of V.  LDV must be at least N and */
00116 /*          at least 1. */
00117 
00118 /*  WORK    (workspace) COMPLEX array, dimension (3*N**2) */
00119 
00120 /*  RESULT  (output) REAL */
00121 /*          The value RESULT, It is currently limited to 1/ulp, to */
00122 /*          avoid overflow. Errors are flagged by RESULT=10/ulp. */
00123 
00124 /*  ===================================================================== */
00125 
00126 /*     .. Parameters .. */
00127 /*     .. */
00128 /*     .. Local Scalars .. */
00129 /*     .. */
00130 /*     .. Local Arrays .. */
00131 /*     .. */
00132 /*     .. External Functions .. */
00133 /*     .. */
00134 /*     .. External Subroutines .. */
00135 /*     .. */
00136 /*     .. Intrinsic Functions .. */
00137 /*     .. */
00138 /*     .. Executable Statements .. */
00139 
00140     /* Parameter adjustments */
00141     a_dim1 = *lda;
00142     a_offset = 1 + a_dim1;
00143     a -= a_offset;
00144     b_dim1 = *ldb;
00145     b_offset = 1 + b_dim1;
00146     b -= b_offset;
00147     s_dim1 = *lds;
00148     s_offset = 1 + s_dim1;
00149     s -= s_offset;
00150     t_dim1 = *ldt;
00151     t_offset = 1 + t_dim1;
00152     t -= t_offset;
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 
00161     /* Function Body */
00162     *result = 0.f;
00163     if (*n <= 0) {
00164         return 0;
00165     }
00166 
00167 /*     Constants */
00168 
00169     unfl = slamch_("Safe minimum");
00170     ulp = slamch_("Epsilon") * slamch_("Base");
00171 
00172 /*     compute the norm of (A,B) */
00173 
00174     clacpy_("Full", n, n, &a[a_offset], lda, &work[1], n);
00175     clacpy_("Full", n, n, &b[b_offset], ldb, &work[*n * *n + 1], n)
00176             ;
00177 /* Computing MAX */
00178     i__1 = *n << 1;
00179     r__1 = clange_("1", n, &i__1, &work[1], n, dum);
00180     abnorm = dmax(r__1,unfl);
00181 
00182 /*     Compute W1 = A - U*S*V', and put in the array WORK(1:N*N) */
00183 
00184     clacpy_(" ", n, n, &a[a_offset], lda, &work[1], n);
00185     cgemm_("N", "N", n, n, n, &c_b2, &u[u_offset], ldu, &s[s_offset], lds, &
00186             c_b1, &work[*n * *n + 1], n);
00187 
00188     q__1.r = -1.f, q__1.i = -0.f;
00189     cgemm_("N", "C", n, n, n, &q__1, &work[*n * *n + 1], n, &v[v_offset], ldv, 
00190              &c_b2, &work[1], n);
00191 
00192 /*     Compute W2 = B - U*T*V', and put in the workarray W(N*N+1:2*N*N) */
00193 
00194     clacpy_(" ", n, n, &b[b_offset], ldb, &work[*n * *n + 1], n);
00195     cgemm_("N", "N", n, n, n, &c_b2, &u[u_offset], ldu, &t[t_offset], ldt, &
00196             c_b1, &work[(*n << 1) * *n + 1], n);
00197 
00198     q__1.r = -1.f, q__1.i = -0.f;
00199     cgemm_("N", "C", n, n, n, &q__1, &work[(*n << 1) * *n + 1], n, &v[
00200             v_offset], ldv, &c_b2, &work[*n * *n + 1], n);
00201 
00202 /*     Compute norm(W)/ ( ulp*norm((A,B)) ) */
00203 
00204     i__1 = *n << 1;
00205     wnorm = clange_("1", n, &i__1, &work[1], n, dum);
00206 
00207     if (abnorm > wnorm) {
00208         *result = wnorm / abnorm / ((*n << 1) * ulp);
00209     } else {
00210         if (abnorm < 1.f) {
00211 /* Computing MIN */
00212             r__1 = wnorm, r__2 = (*n << 1) * abnorm;
00213             *result = dmin(r__1,r__2) / abnorm / ((*n << 1) * ulp);
00214         } else {
00215 /* Computing MIN */
00216             r__1 = wnorm / abnorm, r__2 = (real) (*n << 1);
00217             *result = dmin(r__1,r__2) / ((*n << 1) * ulp);
00218         }
00219     }
00220 
00221     return 0;
00222 
00223 /*     End of CGET54 */
00224 
00225 } /* cget54_ */


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