zqlt02.c
Go to the documentation of this file.
00001 /* zqlt02.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 /* Common Block Declarations */
00017 
00018 struct {
00019     char srnamt[32];
00020 } srnamc_;
00021 
00022 #define srnamc_1 srnamc_
00023 
00024 /* Table of constant values */
00025 
00026 static doublecomplex c_b1 = {-1e10,-1e10};
00027 static doublecomplex c_b9 = {0.,0.};
00028 static doublecomplex c_b14 = {-1.,0.};
00029 static doublecomplex c_b15 = {1.,0.};
00030 static doublereal c_b23 = -1.;
00031 static doublereal c_b24 = 1.;
00032 
00033 /* Subroutine */ int zqlt02_(integer *m, integer *n, integer *k, 
00034         doublecomplex *a, doublecomplex *af, doublecomplex *q, doublecomplex *
00035         l, integer *lda, doublecomplex *tau, doublecomplex *work, integer *
00036         lwork, doublereal *rwork, doublereal *result)
00037 {
00038     /* System generated locals */
00039     integer a_dim1, a_offset, af_dim1, af_offset, l_dim1, l_offset, q_dim1, 
00040             q_offset, i__1, i__2;
00041 
00042     /* Builtin functions */
00043     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00044 
00045     /* Local variables */
00046     doublereal eps;
00047     integer info;
00048     doublereal resid, anorm;
00049     extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 
00050             integer *, doublecomplex *, doublecomplex *, integer *, 
00051             doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
00052             integer *), zherk_(char *, char *, integer *, 
00053             integer *, doublereal *, doublecomplex *, integer *, doublereal *, 
00054              doublecomplex *, integer *);
00055     extern doublereal dlamch_(char *), zlange_(char *, integer *, 
00056             integer *, doublecomplex *, integer *, doublereal *);
00057     extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 
00058             doublecomplex *, integer *, doublecomplex *, integer *), 
00059             zlaset_(char *, integer *, integer *, doublecomplex *, 
00060             doublecomplex *, doublecomplex *, integer *);
00061     extern doublereal zlansy_(char *, char *, integer *, doublecomplex *, 
00062             integer *, doublereal *);
00063     extern /* Subroutine */ int zungql_(integer *, integer *, integer *, 
00064             doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
00065             integer *, integer *);
00066 
00067 
00068 /*  -- LAPACK test routine (version 3.1) -- */
00069 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00070 /*     November 2006 */
00071 
00072 /*     .. Scalar Arguments .. */
00073 /*     .. */
00074 /*     .. Array Arguments .. */
00075 /*     .. */
00076 
00077 /*  Purpose */
00078 /*  ======= */
00079 
00080 /*  ZQLT02 tests ZUNGQL, which generates an m-by-n matrix Q with */
00081 /*  orthonornmal columns that is defined as the product of k elementary */
00082 /*  reflectors. */
00083 
00084 /*  Given the QL factorization of an m-by-n matrix A, ZQLT02 generates */
00085 /*  the orthogonal matrix Q defined by the factorization of the last k */
00086 /*  columns of A; it compares L(m-n+1:m,n-k+1:n) with */
00087 /*  Q(1:m,m-n+1:m)'*A(1:m,n-k+1:n), and checks that the columns of Q are */
00088 /*  orthonormal. */
00089 
00090 /*  Arguments */
00091 /*  ========= */
00092 
00093 /*  M       (input) INTEGER */
00094 /*          The number of rows of the matrix Q to be generated.  M >= 0. */
00095 
00096 /*  N       (input) INTEGER */
00097 /*          The number of columns of the matrix Q to be generated. */
00098 /*          M >= N >= 0. */
00099 
00100 /*  K       (input) INTEGER */
00101 /*          The number of elementary reflectors whose product defines the */
00102 /*          matrix Q. N >= K >= 0. */
00103 
00104 /*  A       (input) COMPLEX*16 array, dimension (LDA,N) */
00105 /*          The m-by-n matrix A which was factorized by ZQLT01. */
00106 
00107 /*  AF      (input) COMPLEX*16 array, dimension (LDA,N) */
00108 /*          Details of the QL factorization of A, as returned by ZGEQLF. */
00109 /*          See ZGEQLF for further details. */
00110 
00111 /*  Q       (workspace) COMPLEX*16 array, dimension (LDA,N) */
00112 
00113 /*  L       (workspace) COMPLEX*16 array, dimension (LDA,N) */
00114 
00115 /*  LDA     (input) INTEGER */
00116 /*          The leading dimension of the arrays A, AF, Q and L. LDA >= M. */
00117 
00118 /*  TAU     (input) COMPLEX*16 array, dimension (N) */
00119 /*          The scalar factors of the elementary reflectors corresponding */
00120 /*          to the QL factorization in AF. */
00121 
00122 /*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) */
00123 
00124 /*  LWORK   (input) INTEGER */
00125 /*          The dimension of the array WORK. */
00126 
00127 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (M) */
00128 
00129 /*  RESULT  (output) DOUBLE PRECISION array, dimension (2) */
00130 /*          The test ratios: */
00131 /*          RESULT(1) = norm( L - Q'*A ) / ( M * norm(A) * EPS ) */
00132 /*          RESULT(2) = norm( I - Q'*Q ) / ( M * EPS ) */
00133 
00134 /*  ===================================================================== */
00135 
00136 /*     .. Parameters .. */
00137 /*     .. */
00138 /*     .. Local Scalars .. */
00139 /*     .. */
00140 /*     .. External Functions .. */
00141 /*     .. */
00142 /*     .. External Subroutines .. */
00143 /*     .. */
00144 /*     .. Intrinsic Functions .. */
00145 /*     .. */
00146 /*     .. Scalars in Common .. */
00147 /*     .. */
00148 /*     .. Common blocks .. */
00149 /*     .. */
00150 /*     .. Executable Statements .. */
00151 
00152 /*     Quick return if possible */
00153 
00154     /* Parameter adjustments */
00155     l_dim1 = *lda;
00156     l_offset = 1 + l_dim1;
00157     l -= l_offset;
00158     q_dim1 = *lda;
00159     q_offset = 1 + q_dim1;
00160     q -= q_offset;
00161     af_dim1 = *lda;
00162     af_offset = 1 + af_dim1;
00163     af -= af_offset;
00164     a_dim1 = *lda;
00165     a_offset = 1 + a_dim1;
00166     a -= a_offset;
00167     --tau;
00168     --work;
00169     --rwork;
00170     --result;
00171 
00172     /* Function Body */
00173     if (*m == 0 || *n == 0 || *k == 0) {
00174         result[1] = 0.;
00175         result[2] = 0.;
00176         return 0;
00177     }
00178 
00179     eps = dlamch_("Epsilon");
00180 
00181 /*     Copy the last k columns of the factorization to the array Q */
00182 
00183     zlaset_("Full", m, n, &c_b1, &c_b1, &q[q_offset], lda);
00184     if (*k < *m) {
00185         i__1 = *m - *k;
00186         zlacpy_("Full", &i__1, k, &af[(*n - *k + 1) * af_dim1 + 1], lda, &q[(*
00187                 n - *k + 1) * q_dim1 + 1], lda);
00188     }
00189     if (*k > 1) {
00190         i__1 = *k - 1;
00191         i__2 = *k - 1;
00192         zlacpy_("Upper", &i__1, &i__2, &af[*m - *k + 1 + (*n - *k + 2) * 
00193                 af_dim1], lda, &q[*m - *k + 1 + (*n - *k + 2) * q_dim1], lda);
00194     }
00195 
00196 /*     Generate the last n columns of the matrix Q */
00197 
00198     s_copy(srnamc_1.srnamt, "ZUNGQL", (ftnlen)32, (ftnlen)6);
00199     zungql_(m, n, k, &q[q_offset], lda, &tau[*n - *k + 1], &work[1], lwork, &
00200             info);
00201 
00202 /*     Copy L(m-n+1:m,n-k+1:n) */
00203 
00204     zlaset_("Full", n, k, &c_b9, &c_b9, &l[*m - *n + 1 + (*n - *k + 1) * 
00205             l_dim1], lda);
00206     zlacpy_("Lower", k, k, &af[*m - *k + 1 + (*n - *k + 1) * af_dim1], lda, &
00207             l[*m - *k + 1 + (*n - *k + 1) * l_dim1], lda);
00208 
00209 /*     Compute L(m-n+1:m,n-k+1:n) - Q(1:m,m-n+1:m)' * A(1:m,n-k+1:n) */
00210 
00211     zgemm_("Conjugate transpose", "No transpose", n, k, m, &c_b14, &q[
00212             q_offset], lda, &a[(*n - *k + 1) * a_dim1 + 1], lda, &c_b15, &l[*
00213             m - *n + 1 + (*n - *k + 1) * l_dim1], lda)
00214             ;
00215 
00216 /*     Compute norm( L - Q'*A ) / ( M * norm(A) * EPS ) . */
00217 
00218     anorm = zlange_("1", m, k, &a[(*n - *k + 1) * a_dim1 + 1], lda, &rwork[1]);
00219     resid = zlange_("1", n, k, &l[*m - *n + 1 + (*n - *k + 1) * l_dim1], lda, 
00220             &rwork[1]);
00221     if (anorm > 0.) {
00222         result[1] = resid / (doublereal) max(1,*m) / anorm / eps;
00223     } else {
00224         result[1] = 0.;
00225     }
00226 
00227 /*     Compute I - Q'*Q */
00228 
00229     zlaset_("Full", n, n, &c_b9, &c_b15, &l[l_offset], lda);
00230     zherk_("Upper", "Conjugate transpose", n, m, &c_b23, &q[q_offset], lda, &
00231             c_b24, &l[l_offset], lda);
00232 
00233 /*     Compute norm( I - Q'*Q ) / ( M * EPS ) . */
00234 
00235     resid = zlansy_("1", "Upper", n, &l[l_offset], lda, &rwork[1]);
00236 
00237     result[2] = resid / (doublereal) max(1,*m) / eps;
00238 
00239     return 0;
00240 
00241 /*     End of ZQLT02 */
00242 
00243 } /* zqlt02_ */


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