cgeqp3.c
Go to the documentation of this file.
00001 /* cgeqp3.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__1 = 1;
00019 static integer c_n1 = -1;
00020 static integer c__3 = 3;
00021 static integer c__2 = 2;
00022 
00023 /* Subroutine */ int cgeqp3_(integer *m, integer *n, complex *a, integer *lda, 
00024          integer *jpvt, complex *tau, complex *work, integer *lwork, real *
00025         rwork, integer *info)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, i__1, i__2, i__3;
00029 
00030     /* Local variables */
00031     integer j, jb, na, nb, sm, sn, nx, fjb, iws, nfxd, nbmin;
00032     extern /* Subroutine */ int cswap_(integer *, complex *, integer *, 
00033             complex *, integer *);
00034     integer minmn, minws;
00035     extern /* Subroutine */ int claqp2_(integer *, integer *, integer *, 
00036             complex *, integer *, integer *, complex *, real *, real *, 
00037             complex *);
00038     extern doublereal scnrm2_(integer *, complex *, integer *);
00039     extern /* Subroutine */ int cgeqrf_(integer *, integer *, complex *, 
00040             integer *, complex *, complex *, integer *, integer *), xerbla_(
00041             char *, integer *);
00042     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00043             integer *, integer *);
00044     extern /* Subroutine */ int claqps_(integer *, integer *, integer *, 
00045             integer *, integer *, complex *, integer *, integer *, complex *, 
00046             real *, real *, complex *, complex *, integer *);
00047     integer topbmn, sminmn;
00048     extern /* Subroutine */ int cunmqr_(char *, char *, integer *, integer *, 
00049             integer *, complex *, integer *, complex *, complex *, integer *, 
00050             complex *, integer *, integer *);
00051     integer lwkopt;
00052     logical lquery;
00053 
00054 
00055 /*  -- LAPACK routine (version 3.2) -- */
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 /*  CGEQP3 computes a QR factorization with column pivoting of a */
00068 /*  matrix A:  A*P = Q*R  using Level 3 BLAS. */
00069 
00070 /*  Arguments */
00071 /*  ========= */
00072 
00073 /*  M       (input) INTEGER */
00074 /*          The number of rows of the matrix A. M >= 0. */
00075 
00076 /*  N       (input) INTEGER */
00077 /*          The number of columns of the matrix A.  N >= 0. */
00078 
00079 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00080 /*          On entry, the M-by-N matrix A. */
00081 /*          On exit, the upper triangle of the array contains the */
00082 /*          min(M,N)-by-N upper trapezoidal matrix R; the elements below */
00083 /*          the diagonal, together with the array TAU, represent the */
00084 /*          unitary matrix Q as a product of min(M,N) elementary */
00085 /*          reflectors. */
00086 
00087 /*  LDA     (input) INTEGER */
00088 /*          The leading dimension of the array A. LDA >= max(1,M). */
00089 
00090 /*  JPVT    (input/output) INTEGER array, dimension (N) */
00091 /*          On entry, if JPVT(J).ne.0, the J-th column of A is permuted */
00092 /*          to the front of A*P (a leading column); if JPVT(J)=0, */
00093 /*          the J-th column of A is a free column. */
00094 /*          On exit, if JPVT(J)=K, then the J-th column of A*P was the */
00095 /*          the K-th column of A. */
00096 
00097 /*  TAU     (output) COMPLEX array, dimension (min(M,N)) */
00098 /*          The scalar factors of the elementary reflectors. */
00099 
00100 /*  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */
00101 /*          On exit, if INFO=0, WORK(1) returns the optimal LWORK. */
00102 
00103 /*  LWORK   (input) INTEGER */
00104 /*          The dimension of the array WORK. LWORK >= N+1. */
00105 /*          For optimal performance LWORK >= ( N+1 )*NB, where NB */
00106 /*          is the optimal blocksize. */
00107 
00108 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00109 /*          only calculates the optimal size of the WORK array, returns */
00110 /*          this value as the first entry of the WORK array, and no error */
00111 /*          message related to LWORK is issued by XERBLA. */
00112 
00113 /*  RWORK   (workspace) REAL array, dimension (2*N) */
00114 
00115 /*  INFO    (output) INTEGER */
00116 /*          = 0: successful exit. */
00117 /*          < 0: if INFO = -i, the i-th argument had an illegal value. */
00118 
00119 /*  Further Details */
00120 /*  =============== */
00121 
00122 /*  The matrix Q is represented as a product of elementary reflectors */
00123 
00124 /*     Q = H(1) H(2) . . . H(k), where k = min(m,n). */
00125 
00126 /*  Each H(i) has the form */
00127 
00128 /*     H(i) = I - tau * v * v' */
00129 
00130 /*  where tau is a real/complex scalar, and v is a real/complex vector */
00131 /*  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in */
00132 /*  A(i+1:m,i), and tau in TAU(i). */
00133 
00134 /*  Based on contributions by */
00135 /*    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain */
00136 /*    X. Sun, Computer Science Dept., Duke University, USA */
00137 
00138 /*  ===================================================================== */
00139 
00140 /*     .. Parameters .. */
00141 /*     .. */
00142 /*     .. Local Scalars .. */
00143 /*     .. */
00144 /*     .. External Subroutines .. */
00145 /*     .. */
00146 /*     .. External Functions .. */
00147 /*     .. */
00148 /*     .. Intrinsic Functions .. */
00149 /*     .. */
00150 /*     .. Executable Statements .. */
00151 
00152 /*     Test input arguments */
00153 /*     ==================== */
00154 
00155     /* Parameter adjustments */
00156     a_dim1 = *lda;
00157     a_offset = 1 + a_dim1;
00158     a -= a_offset;
00159     --jpvt;
00160     --tau;
00161     --work;
00162     --rwork;
00163 
00164     /* Function Body */
00165     *info = 0;
00166     lquery = *lwork == -1;
00167     if (*m < 0) {
00168         *info = -1;
00169     } else if (*n < 0) {
00170         *info = -2;
00171     } else if (*lda < max(1,*m)) {
00172         *info = -4;
00173     }
00174 
00175     if (*info == 0) {
00176         minmn = min(*m,*n);
00177         if (minmn == 0) {
00178             iws = 1;
00179             lwkopt = 1;
00180         } else {
00181             iws = *n + 1;
00182             nb = ilaenv_(&c__1, "CGEQRF", " ", m, n, &c_n1, &c_n1);
00183             lwkopt = (*n + 1) * nb;
00184         }
00185         work[1].r = (real) lwkopt, work[1].i = 0.f;
00186 
00187         if (*lwork < iws && ! lquery) {
00188             *info = -8;
00189         }
00190     }
00191 
00192     if (*info != 0) {
00193         i__1 = -(*info);
00194         xerbla_("CGEQP3", &i__1);
00195         return 0;
00196     } else if (lquery) {
00197         return 0;
00198     }
00199 
00200 /*     Quick return if possible. */
00201 
00202     if (minmn == 0) {
00203         return 0;
00204     }
00205 
00206 /*     Move initial columns up front. */
00207 
00208     nfxd = 1;
00209     i__1 = *n;
00210     for (j = 1; j <= i__1; ++j) {
00211         if (jpvt[j] != 0) {
00212             if (j != nfxd) {
00213                 cswap_(m, &a[j * a_dim1 + 1], &c__1, &a[nfxd * a_dim1 + 1], &
00214                         c__1);
00215                 jpvt[j] = jpvt[nfxd];
00216                 jpvt[nfxd] = j;
00217             } else {
00218                 jpvt[j] = j;
00219             }
00220             ++nfxd;
00221         } else {
00222             jpvt[j] = j;
00223         }
00224 /* L10: */
00225     }
00226     --nfxd;
00227 
00228 /*     Factorize fixed columns */
00229 /*     ======================= */
00230 
00231 /*     Compute the QR factorization of fixed columns and update */
00232 /*     remaining columns. */
00233 
00234     if (nfxd > 0) {
00235         na = min(*m,nfxd);
00236 /* CC      CALL CGEQR2( M, NA, A, LDA, TAU, WORK, INFO ) */
00237         cgeqrf_(m, &na, &a[a_offset], lda, &tau[1], &work[1], lwork, info);
00238 /* Computing MAX */
00239         i__1 = iws, i__2 = (integer) work[1].r;
00240         iws = max(i__1,i__2);
00241         if (na < *n) {
00242 /* CC         CALL CUNM2R( 'Left', 'Conjugate Transpose', M, N-NA, */
00243 /* CC  $                   NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK, */
00244 /* CC  $                   INFO ) */
00245             i__1 = *n - na;
00246             cunmqr_("Left", "Conjugate Transpose", m, &i__1, &na, &a[a_offset]
00247 , lda, &tau[1], &a[(na + 1) * a_dim1 + 1], lda, &work[1], 
00248                     lwork, info);
00249 /* Computing MAX */
00250             i__1 = iws, i__2 = (integer) work[1].r;
00251             iws = max(i__1,i__2);
00252         }
00253     }
00254 
00255 /*     Factorize free columns */
00256 /*     ====================== */
00257 
00258     if (nfxd < minmn) {
00259 
00260         sm = *m - nfxd;
00261         sn = *n - nfxd;
00262         sminmn = minmn - nfxd;
00263 
00264 /*        Determine the block size. */
00265 
00266         nb = ilaenv_(&c__1, "CGEQRF", " ", &sm, &sn, &c_n1, &c_n1);
00267         nbmin = 2;
00268         nx = 0;
00269 
00270         if (nb > 1 && nb < sminmn) {
00271 
00272 /*           Determine when to cross over from blocked to unblocked code. */
00273 
00274 /* Computing MAX */
00275             i__1 = 0, i__2 = ilaenv_(&c__3, "CGEQRF", " ", &sm, &sn, &c_n1, &
00276                     c_n1);
00277             nx = max(i__1,i__2);
00278 
00279 
00280             if (nx < sminmn) {
00281 
00282 /*              Determine if workspace is large enough for blocked code. */
00283 
00284                 minws = (sn + 1) * nb;
00285                 iws = max(iws,minws);
00286                 if (*lwork < minws) {
00287 
00288 /*                 Not enough workspace to use optimal NB: Reduce NB and */
00289 /*                 determine the minimum value of NB. */
00290 
00291                     nb = *lwork / (sn + 1);
00292 /* Computing MAX */
00293                     i__1 = 2, i__2 = ilaenv_(&c__2, "CGEQRF", " ", &sm, &sn, &
00294                             c_n1, &c_n1);
00295                     nbmin = max(i__1,i__2);
00296 
00297 
00298                 }
00299             }
00300         }
00301 
00302 /*        Initialize partial column norms. The first N elements of work */
00303 /*        store the exact column norms. */
00304 
00305         i__1 = *n;
00306         for (j = nfxd + 1; j <= i__1; ++j) {
00307             rwork[j] = scnrm2_(&sm, &a[nfxd + 1 + j * a_dim1], &c__1);
00308             rwork[*n + j] = rwork[j];
00309 /* L20: */
00310         }
00311 
00312         if (nb >= nbmin && nb < sminmn && nx < sminmn) {
00313 
00314 /*           Use blocked code initially. */
00315 
00316             j = nfxd + 1;
00317 
00318 /*           Compute factorization: while loop. */
00319 
00320 
00321             topbmn = minmn - nx;
00322 L30:
00323             if (j <= topbmn) {
00324 /* Computing MIN */
00325                 i__1 = nb, i__2 = topbmn - j + 1;
00326                 jb = min(i__1,i__2);
00327 
00328 /*              Factorize JB columns among columns J:N. */
00329 
00330                 i__1 = *n - j + 1;
00331                 i__2 = j - 1;
00332                 i__3 = *n - j + 1;
00333                 claqps_(m, &i__1, &i__2, &jb, &fjb, &a[j * a_dim1 + 1], lda, &
00334                         jpvt[j], &tau[j], &rwork[j], &rwork[*n + j], &work[1], 
00335                          &work[jb + 1], &i__3);
00336 
00337                 j += fjb;
00338                 goto L30;
00339             }
00340         } else {
00341             j = nfxd + 1;
00342         }
00343 
00344 /*        Use unblocked code to factor the last or only block. */
00345 
00346 
00347         if (j <= minmn) {
00348             i__1 = *n - j + 1;
00349             i__2 = j - 1;
00350             claqp2_(m, &i__1, &i__2, &a[j * a_dim1 + 1], lda, &jpvt[j], &tau[
00351                     j], &rwork[j], &rwork[*n + j], &work[1]);
00352         }
00353 
00354     }
00355 
00356     work[1].r = (real) iws, work[1].i = 0.f;
00357     return 0;
00358 
00359 /*     End of CGEQP3 */
00360 
00361 } /* cgeqp3_ */


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