cqpt01.c
Go to the documentation of this file.
00001 /* cqpt01.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__10 = 10;
00019 static integer c__1 = 1;
00020 static complex c_b16 = {-1.f,0.f};
00021 
00022 doublereal cqpt01_(integer *m, integer *n, integer *k, complex *a, complex *
00023         af, integer *lda, complex *tau, integer *jpvt, complex *work, integer 
00024         *lwork)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, af_dim1, af_offset, i__1, i__2, i__3, i__4;
00028     real ret_val;
00029 
00030     /* Local variables */
00031     integer i__, j, info;
00032     real norma;
00033     extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
00034             complex *, integer *), caxpy_(integer *, complex *, complex *, 
00035             integer *, complex *, integer *);
00036     real rwork[1];
00037     extern doublereal clange_(char *, integer *, integer *, complex *, 
00038             integer *, real *), slamch_(char *);
00039     extern /* Subroutine */ int xerbla_(char *, integer *), cunmqr_(
00040             char *, char *, integer *, integer *, integer *, complex *, 
00041             integer *, complex *, complex *, integer *, complex *, integer *, 
00042             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 /*  CQPT01 tests the QR-factorization with pivoting of a matrix A.  The */
00058 /*  array AF contains the (possibly partial) QR-factorization of A, where */
00059 /*  the upper triangle of AF(1:k,1:k) is a partial triangular factor, */
00060 /*  the entries below the diagonal in the first k columns are the */
00061 /*  Householder vectors, and the rest of AF contains a partially updated */
00062 /*  matrix. */
00063 
00064 /*  This function returns ||A*P - Q*R||/(||norm(A)||*eps*M) */
00065 
00066 /*  Arguments */
00067 /*  ========= */
00068 
00069 /*  M       (input) INTEGER */
00070 /*          The number of rows of the matrices A and AF. */
00071 
00072 /*  N       (input) INTEGER */
00073 /*          The number of columns of the matrices A and AF. */
00074 
00075 /*  K       (input) INTEGER */
00076 /*          The number of columns of AF that have been reduced */
00077 /*          to upper triangular form. */
00078 
00079 /*  A       (input) COMPLEX array, dimension (LDA, N) */
00080 /*          The original matrix A. */
00081 
00082 /*  AF      (input) COMPLEX array, dimension (LDA,N) */
00083 /*          The (possibly partial) output of CGEQPF.  The upper triangle */
00084 /*          of AF(1:k,1:k) is a partial triangular factor, the entries */
00085 /*          below the diagonal in the first k columns are the Householder */
00086 /*          vectors, and the rest of AF contains a partially updated */
00087 /*          matrix. */
00088 
00089 /*  LDA     (input) INTEGER */
00090 /*          The leading dimension of the arrays A and AF. */
00091 
00092 /*  TAU     (input) COMPLEX array, dimension (K) */
00093 /*          Details of the Householder transformations as returned by */
00094 /*          CGEQPF. */
00095 
00096 /*  JPVT    (input) INTEGER array, dimension (N) */
00097 /*          Pivot information as returned by CGEQPF. */
00098 
00099 /*  WORK    (workspace) COMPLEX array, dimension (LWORK) */
00100 
00101 /*  LWORK   (input) INTEGER */
00102 /*          The length of the array WORK.  LWORK >= M*N+N. */
00103 
00104 /*  ===================================================================== */
00105 
00106 /*     .. Parameters .. */
00107 /*     .. */
00108 /*     .. Local Scalars .. */
00109 /*     .. */
00110 /*     .. Local Arrays .. */
00111 /*     .. */
00112 /*     .. External Functions .. */
00113 /*     .. */
00114 /*     .. External Subroutines .. */
00115 /*     .. */
00116 /*     .. Intrinsic Functions .. */
00117 /*     .. */
00118 /*     .. Executable Statements .. */
00119 
00120     /* Parameter adjustments */
00121     af_dim1 = *lda;
00122     af_offset = 1 + af_dim1;
00123     af -= af_offset;
00124     a_dim1 = *lda;
00125     a_offset = 1 + a_dim1;
00126     a -= a_offset;
00127     --tau;
00128     --jpvt;
00129     --work;
00130 
00131     /* Function Body */
00132     ret_val = 0.f;
00133 
00134 /*     Test if there is enough workspace */
00135 
00136     if (*lwork < *m * *n + *n) {
00137         xerbla_("CQPT01", &c__10);
00138         return ret_val;
00139     }
00140 
00141 /*     Quick return if possible */
00142 
00143     if (*m <= 0 || *n <= 0) {
00144         return ret_val;
00145     }
00146 
00147     norma = clange_("One-norm", m, n, &a[a_offset], lda, rwork);
00148 
00149     i__1 = *k;
00150     for (j = 1; j <= i__1; ++j) {
00151         i__2 = min(j,*m);
00152         for (i__ = 1; i__ <= i__2; ++i__) {
00153             i__3 = (j - 1) * *m + i__;
00154             i__4 = i__ + j * af_dim1;
00155             work[i__3].r = af[i__4].r, work[i__3].i = af[i__4].i;
00156 /* L10: */
00157         }
00158         i__2 = *m;
00159         for (i__ = j + 1; i__ <= i__2; ++i__) {
00160             i__3 = (j - 1) * *m + i__;
00161             work[i__3].r = 0.f, work[i__3].i = 0.f;
00162 /* L20: */
00163         }
00164 /* L30: */
00165     }
00166     i__1 = *n;
00167     for (j = *k + 1; j <= i__1; ++j) {
00168         ccopy_(m, &af[j * af_dim1 + 1], &c__1, &work[(j - 1) * *m + 1], &c__1)
00169                 ;
00170 /* L40: */
00171     }
00172 
00173     i__1 = *lwork - *m * *n;
00174     cunmqr_("Left", "No transpose", m, n, k, &af[af_offset], lda, &tau[1], &
00175             work[1], m, &work[*m * *n + 1], &i__1, &info);
00176 
00177     i__1 = *n;
00178     for (j = 1; j <= i__1; ++j) {
00179 
00180 /*        Compare i-th column of QR and jpvt(i)-th column of A */
00181 
00182         caxpy_(m, &c_b16, &a[jpvt[j] * a_dim1 + 1], &c__1, &work[(j - 1) * *m 
00183                 + 1], &c__1);
00184 /* L50: */
00185     }
00186 
00187     ret_val = clange_("One-norm", m, n, &work[1], m, rwork) / ((
00188             real) max(*m,*n) * slamch_("Epsilon"));
00189     if (norma != 0.f) {
00190         ret_val /= norma;
00191     }
00192 
00193     return ret_val;
00194 
00195 /*     End of CQPT01 */
00196 
00197 } /* cqpt01_ */


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