text: /* cmylap.f -- translated by f2c (version 20041007).
     You must link the resulting object file with libf2c:
  	on Microsoft Windows system, link with libf2c.lib;
  	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
  	or, if you install libf2c.a in a standard place, with -lf2c -lm
  	-- in that order, at the end of the command line, as in
  		cc *.o -lf2c -lm
  	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
  
  		http://www.netlib.org/f2c/libf2c.zip
  */
  
  #include "rrqr.h"
  
  /* Table of constant values */
  
  static integer c__1 = 1;
  static integer c_n1 = -1;
  static integer c__3 = 3;
  static integer c__2 = 2;
  
  /* ******************************************************************** */
  /* Subroutine */ int cgeqpf(integer *m, integer *n, complex *a, integer *lda,
  	 integer *jpvt, complex *tau, complex *work, real *rwork, integer *
  	info)
  {
      /* System generated locals */
      integer a_dim1, a_offset, i__1, i__2, i__3;
      real r__1;
      complex q__1;
  
      /* Builtin functions */
      void r_cnjg(complex *, complex *);
      double c_abs(complex *), sqrt(doublereal);
  
      /* Local variables */
      static integer i__, j, ma, mn;
      static complex aii;
      static integer pvt;
      static real temp, temp2;
      static integer itemp;
  
  
  /*  -- LAPACK auxiliary routine (version 1.0b) -- */
  /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
  /*     Courant Institute, Argonne National Lab, and Rice University */
  /*     February 29, 1992 */
  
  /*     .. Scalar Arguments .. */
  /*     .. */
  /*     .. Array Arguments .. */
  /*     .. */
  
  /*  Purpose */
  /*  ======= */
  
  /*  CGEQPF computes a QR factorization with column pivoting of a */
  /*  complex m by n matrix A: A*P = Q*R */
  
  /*  Arguments */
  /*  ========= */
  
  /*  M       (input) INTEGER */
  /*          The number of rows of the matrix A. M >= 0. */
  
  /*  N       (input) INTEGER */
  /*          The number of columns of the matrix A. N >= 0 */
  
  /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
  /*          On entry, the m by n matrix A. */
  /*          On exit, the upper triangle of the array contains the */
  /*          min(m,n) by n upper triangular matrix R; the elements */
  /*          below the diagonal, together with the array TAU, */
  /*          represent the orthogonal matrix Q as a product of */
  /*          min(m,n) elementary reflectors. */
  
  /*  LDA     (input) INTEGER */
  /*          The leading dimension of the array A. LDA >= max(1,M). */
  
  /*  JPVT    (input/output) INTEGER array, dimension (N) */
  /*          on entry: If JPVT(I) <> 0, column I of A is permuted */
  /*          to the front of AP (a leading column) */
  /*          IF JPVT(I) == 0, column I of A is a free column. */
  /*          on exit: If JPVT(I) = K, then the Ith column of AP */
  /*          was the Kth column of A. */
  
  /*  TAU     (output) COMPLEX array, dimension (min(M,N)) */
  /*          Stores further details of */
  /*          the orthogonal matrix Q (see A). */
  
  /*  WORK    (workspace) COMPLEX array, dimension (N) */
  
  /*  RWORK   (workspace) REAL array, dimension (2*N) */
  
  /*  INFO    (output) INTEGER */
  /*          = 0: successful exit */
  /*          < 0: if INFO = -i, the i-th argument had an illegal value */
  
  /*  Further Details */
  /*  =============== */
  
  /*  The matrix Q is represented as a product of elementary reflectors */
  
  /*     Q = H(1) H(2) . . . H(n) */
  
  /*  Each H(i) has the form */
  
  /*     H = I - tau * v * v' */
  
  /*  where tau is a complex scalar, and v is a complex vector with */
  /*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i). */
  
  /*  The matrix P is represented in jpvt as follows: If */
  /*     jpvt(j) = i */
  /*  then the jth column of P is the ith canonical unit vector. */
  
  /*  ===================================================================== */
  
  /*     .. Parameters .. */
  /*     .. */
  /*     .. Local Scalars .. */
  /*     .. */
  /*     .. External Subroutines .. */
  /*     .. */
  /*     .. Intrinsic Functions .. */
  /*     .. */
  /*     .. External Functions .. */
  /*     .. */
  /*     .. Executable Statements .. */
  
  /*     Test the input arguments */
  
      /* Parameter adjustments */
      a_dim1 = *lda;
      a_offset = 1 + a_dim1;
      a -= a_offset;
      --jpvt;
      --tau;
      --work;
      --rwork;
  
      /* Function Body */
      *info = 0;
      if (*m < 0) {
  	*info = -1;
      } else if (*n < 0) {
  	*info = -2;
      } else if (*lda < max(1,*m)) {
  	*info = -4;
      }
      if (*info != 0) {
  	i__1 = -(*info);
  	xerbla("CGEQPF", &i__1);
  	return 0;
      }
  
  /*     Quick return if possible */
  
      if (*m == 0 || *n == 0) {
  	return 0;
      }
  
      mn = min(*m,*n);
  
  /*     Move initial columns up front */
  
      itemp = 1;
      i__1 = *n;
      for (i__ = 1; i__ <= i__1; ++i__) {
  	if (jpvt[i__] != 0) {
  	    if (i__ != itemp) {
  		cswap(m, &a[i__ * a_dim1 + 1], &c__1, &a[itemp * a_dim1 + 1],
  			 &c__1);
  		jpvt[i__] = jpvt[itemp];
  		jpvt[itemp] = i__;
  	    } else {
  		jpvt[i__] = i__;
  	    }
  	    ++itemp;
  	} else {
  	    jpvt[i__] = i__;
  	}
  /* L10: */
      }
      --itemp;
  
  /*     Compute the QR factorization and update remaining columns */
  
      if (itemp > 0) {
  	ma = min(itemp,*m);
  	cgeqr2(m, &ma, &a[a_offset], lda, &tau[1], &work[1], info);
  	if (ma < *n) {
  	    i__1 = *n - ma;
  	    cunm2r("Left", "Conjugate transpose", m, &i__1, &ma, &a[a_offset]
  		    , lda, &tau[1], &a[(ma + 1) * a_dim1 + 1], lda, &work[1],
  		    info, (ftnlen)4, (ftnlen)19);
  	}
      }
  
      if (itemp < mn) {
  
  /*        Initialize partial column norms. The first n entries of */
  /*        work store the exact column norms. */
  
  	i__1 = *n;
  	for (i__ = itemp + 1; i__ <= i__1; ++i__) {
  	    i__2 = *m - itemp;
  	    rwork[i__] = scnrm2(&i__2, &a[itemp + 1 + i__ * a_dim1], &c__1);
  	    rwork[*n + i__] = rwork[i__];
  /* L20: */
  	}
  
  /*        Compute factorization */
  
  	i__1 = mn;
  	for (i__ = itemp + 1; i__ <= i__1; ++i__) {
  
  /*           Determine ith pivot column and swap if necessary */
  
  	    i__2 = *n - i__ + 1;
  	    pvt = i__ - 1 + isamax(&i__2, &rwork[i__], &c__1);
  
  	    if (pvt != i__) {
  		cswap(m, &a[pvt * a_dim1 + 1], &c__1, &a[i__ * a_dim1 + 1], &
  			c__1);
  		itemp = jpvt[pvt];
  		jpvt[pvt] = jpvt[i__];
  		jpvt[i__] = itemp;
  		rwork[pvt] = rwork[i__];
  		rwork[*n + pvt] = rwork[*n + i__];
  	    }
  
  /*           Generate elementary reflector H(i) */
  
  	    i__2 = i__ + i__ * a_dim1;
  	    aii.r = a[i__2].r, aii.i = a[i__2].i;
  	    i__2 = *m - i__ + 1;
  /* Computing MIN */
  	    i__3 = i__ + 1;
  	    clarfg(&i__2, &aii, &a[min(i__3,*m) + i__ * a_dim1], &c__1, &tau[
  		    i__]);
  	    i__2 = i__ + i__ * a_dim1;
  	    a[i__2].r = aii.r, a[i__2].i = aii.i;
  
  	    if (i__ < *n) {
  
  /*              Apply H(i) to A(i:m,i+1:n) from the left */
  
  		i__2 = i__ + i__ * a_dim1;
  		aii.r = a[i__2].r, aii.i = a[i__2].i;
  		i__2 = i__ + i__ * a_dim1;
  		a[i__2].r = 1.f, a[i__2].i = 0.f;
  		i__2 = *m - i__ + 1;
  		i__3 = *n - i__;
  		r_cnjg(&q__1, &tau[i__]);
  		clarf("Left", &i__2, &i__3, &a[i__ + i__ * a_dim1], &c__1, &
  			q__1, &a[i__ + (i__ + 1) * a_dim1], lda, &work[1]);
  		i__2 = i__ + i__ * a_dim1;
  		a[i__2].r = aii.r, a[i__2].i = aii.i;
  	    }
  
  /*           Update partial column norms */
  
  	    i__2 = *n;
  	    for (j = i__ + 1; j <= i__2; ++j) {
  		if (rwork[j] != 0.f) {
  /* Computing 2nd power */
  		    r__1 = c_abs(&a[i__ + j * a_dim1]) / rwork[j];
  		    temp = 1.f - r__1 * r__1;
  		    temp = dmax(temp,0.f);
  /* Computing 2nd power */
  		    r__1 = rwork[j] / rwork[*n + j];
  		    temp2 = temp * .05f * (r__1 * r__1) + 1.f;
  		    if (temp2 == 1.f) {
  			if (*m - i__ > 0) {
  			    i__3 = *m - i__;
  			    rwork[j] = scnrm2(&i__3, &a[i__ + 1 + j * a_dim1]
  				    , &c__1);
  			    rwork[*n + j] = rwork[j];
  			} else {
  			    rwork[j] = 0.f;
  			    rwork[*n + j] = 0.f;
  			}
  		    } else {
  			rwork[j] *= sqrt(temp);
  		    }
  		}
  /* L30: */
  	    }
  
  /* L40: */
  	}
      }
      return 0;
  
  /*     End of CGEQPF */
  
  } /* cgeqpf_ */
  
  /* ******************************************************************** */
  /* Subroutine */ int cgeqrf(integer *m, integer *n, complex *a, integer *lda,
  	 complex *tau, complex *work, integer *lwork, integer *info)
  {
      /* System generated locals */
      integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
  
      /* Local variables */
      static integer i__, k, ib, nb, nx, iws, nbmin, iinfo;
      static integer ldwork;
  
  
  /*  -- LAPACK routine (version 1.0b) -- */
  /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
  /*     Courant Institute, Argonne National Lab, and Rice University */
  /*     February 29, 1992 */
  
  /*     .. Scalar Arguments .. */
  /*     .. */
  /*     .. Array Arguments .. */
  /*     .. */
  
  /*  Purpose */
  /*  ======= */
  
  /*  CGEQRF computes a QR factorization of a complex m by n matrix A: */
  /*  A = Q * R. */
  
  /*  Arguments */
  /*  ========= */
  
  /*  M       (input) INTEGER */
  /*          The number of rows of the matrix A.  M >= 0. */
  
  /*  N       (input) INTEGER */
  /*          The number of columns of the matrix A.  N >= 0. */
  
  /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
  /*          On entry, the m by n matrix A. */
  /*          On exit, the elements on and above the diagonal of the array */
  /*          contain the min(m,n) by n upper trapezoidal matrix R (R is */
  /*          upper triangular if m >= n); the elements below the diagonal, */
  /*          with the array TAU, represent the unitary matrix Q as a */
  /*          product of elementary reflectors (see Further Details). */
  
  /*  LDA     (input) INTEGER */
  /*          The leading dimension of the array A.  LDA >= max(1,M). */
  
  /*  TAU     (output) COMPLEX array, dimension (min(M,N)) */
  /*          The scalar factors of the elementary reflectors (see Further */
  /*          Details). */
  
  /*  WORK    (workspace) COMPLEX array, dimension (LWORK) */
  /*          On exit, if INFO = 0, WORK(1) returns the minimum value of */
  /*          LWORK required to use the optimal blocksize. */
  
  /*  LWORK   (input) INTEGER */
  /*          The dimension of the array WORK.  LWORK >= max(1,N). */
  /*          For optimum performance LWORK should be at least N*NB, */
  /*          where NB is the optimal blocksize. */
  
  /*  INFO    (output) INTEGER */
  /*          = 0: successful exit */
  /*          < 0: if INFO = -i, the i-th argument had an illegal value */
  
  /*  Further Details */
  /*  =============== */
  
  /*  The matrix Q is represented as a product of elementary reflectors */
  
  /*     Q = H(1) H(2) . . . H(k), where k = min(m,n). */
  
  /*  Each H(i) has the form */
  
  /*     H(i) = I - tau * v * v' */
  
  /*  where tau is a complex scalar, and v is a complex vector with */
  /*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), */
  /*  and tau in TAU(i). */
  
  /*  ===================================================================== */
  
  /*     .. Local Scalars .. */
  /*     .. */
  /*     .. External Subroutines .. */
  /*     .. */
  /*     .. Intrinsic Functions .. */
  /*     .. */
  /*     .. External Functions .. */
  /*     .. */
  /*     .. Executable Statements .. */
  
  /*     Test the input arguments */
  
      /* Parameter adjustments */
      a_dim1 = *lda;
      a_offset = 1 + a_dim1;
      a -= a_offset;
      --tau;
      --work;
  
      /* Function Body */
      *info = 0;
      if (*m < 0) {
  	*info = -1;
      } else if (*n < 0) {
  	*info = -2;
      } else if (*lda < max(1,*m)) {
  	*info = -4;
      } else if (*lwork < max(1,*n)) {
  	*info = -7;
      }
      if (*info != 0) {
  	i__1 = -(*info);
  	xerbla("CGEQRF", &i__1);
  	return 0;
      }
  
  /*     Quick return if possible */
  
      k = min(*m,*n);
      if (k == 0) {
  	work[1].r = 1.f, work[1].i = 0.f;
  	return 0;
      }
  
  /*     Determine the block size. */
  
      nb = ilaenv(&c__1, "CGEQRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)
  	    1);
      nbmin = 2;
      nx = 0;
      iws = *n;
      if (nb > 1 && nb < k) {
  
  /*        Determine when to cross over from blocked to unblocked code. */
  
  /* Computing MAX */
  	i__1 = 0, i__2 = ilaenv(&c__3, "CGEQRF", " ", m, n, &c_n1, &c_n1, (
  		ftnlen)6, (ftnlen)1);
  	nx = max(i__1,i__2);
  	if (nx < k) {
  
  /*           Determine if workspace is large enough for blocked code. */
  
  	    ldwork = *n;
  	    iws = ldwork * nb;
  	    if (*lwork < iws) {
  
  /*              Not enough workspace to use optimal NB:  reduce NB and */
  /*              determine the minimum value of NB. */
  
  		nb = *lwork / ldwork;
  /* Computing MAX */
  		i__1 = 2, i__2 = ilaenv(&c__2, "CGEQRF", " ", m, n, &c_n1, &
  			c_n1, (ftnlen)6, (ftnlen)1);
  		nbmin = max(i__1,i__2);
  	    }
  	}
      }
  
      if (nb >= nbmin && nb < k && nx < k) {
  
  /*        Use blocked code initially */
  
  	i__1 = k - nx;
  	i__2 = nb;
  	for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
  /* Computing MIN */
  	    i__3 = k - i__ + 1;
  	    ib = min(i__3,nb);
  
  /*           Compute the QR factorization of the current block */
  /*           A(i:m,i:i+ib-1) */
  
  	    i__3 = *m - i__ + 1;
  	    cgeqr2(&i__3, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
  		    1], &iinfo);
  	    if (i__ + ib <= *n) {
  
  /*              Form the triangular factor of the block reflector */
  /*              H = H(i) H(i+1) . . . H(i+ib-1) */
  
  		i__3 = *m - i__ + 1;
  		clarft("Forward", "Columnwise", &i__3, &ib, &a[i__ + i__ *
  			a_dim1], lda, &tau[i__], &work[1], &ldwork);
  
  /*              Apply H' to A(i:m,i+ib:n) from the left */
  
  		i__3 = *m - i__ + 1;
  		i__4 = *n - i__ - ib + 1;
  		clarfb("Left", "Conjugate transpose", "Forward", "Columnwise"
  			, &i__3, &i__4, &ib, &a[i__ + i__ * a_dim1], lda, &
  			work[1], &ldwork, &a[i__ + (i__ + ib) * a_dim1], lda,
  			&work[ib + 1], &ldwork);
  	    }
  /* L10: */
  	}
      } else {
  	i__ = 1;
      }
  
  /*     Use unblocked code to factor the last or only block. */
  
      if (i__ <= k) {
  	i__2 = *m - i__ + 1;
  	i__1 = *n - i__ + 1;
  	cgeqr2(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
  		, &iinfo);
      }
  
      work[1].r = (real) iws, work[1].i = 0.f;
      return 0;
  
  /*     End of CGEQRF */
  
  } /* cgeqrf_ */
  
  /* ******************************************************************** */
  /* Subroutine */ int cgeqr2(integer *m, integer *n, complex *a, integer *lda,
  	 complex *tau, complex *work, integer *info)
  {
      /* System generated locals */
      integer a_dim1, a_offset, i__1, i__2, i__3;
      complex q__1;
  
      /* Builtin functions */
      void r_cnjg(complex *, complex *);
  
      /* Local variables */
      static integer i__, k;
      static complex alpha;
  
  
  /*  -- LAPACK routine (version 1.0b) -- */
  /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
  /*     Courant Institute, Argonne National Lab, and Rice University */
  /*     February 29, 1992 */
  
  /*     .. Scalar Arguments .. */
  /*     .. */
  /*     .. Array Arguments .. */
  /*     .. */
  
  /*  Purpose */
  /*  ======= */
  
  /*  CGEQR2 computes a QR factorization of a complex m by n matrix A: */
  /*  A = Q * R. */
  
  /*  Arguments */
  /*  ========= */
  
  /*  M       (input) INTEGER */
  /*          The number of rows of the matrix A.  M >= 0. */
  
  /*  N       (input) INTEGER */
  /*          The number of columns of the matrix A.  N >= 0. */
  
  /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
  /*          On entry, the m by n matrix A. */
  /*          On exit, the elements on and above the diagonal of the array */
  /*          contain the min(m,n) by n upper trapezoidal matrix R (R is */
  /*          upper triangular if m >= n); the elements below the diagonal, */
  /*          with the array TAU, represent the unitary matrix Q as a */
  /*          product of elementary reflectors (see Further Details). */
  
  /*  LDA     (input) INTEGER */
  /*          The leading dimension of the array A.  LDA >= max(1,M). */
  
  /*  TAU     (output) COMPLEX array, dimension (min(M,N)) */
  /*          The scalar factors of the elementary reflectors (see Further */
  /*          Details). */
  
  /*  WORK    (workspace) COMPLEX array, dimension (N) */
  
  /*  INFO    (output) INTEGER */
  /*          = 0: successful exit */
  /*          < 0: if INFO = -i, the i-th argument had an illegal value */
  
  /*  Further Details */
  /*  =============== */
  
  /*  The matrix Q is represented as a product of elementary reflectors */
  
  /*     Q = H(1) H(2) . . . H(k), where k = min(m,n). */
  
  /*  Each H(i) has the form */
  
  /*     H(i) = I - tau * v * v' */
  
  /*  where tau is a complex scalar, and v is a complex vector with */
  /*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), */
  /*  and tau in TAU(i). */
  
  /*  ===================================================================== */
  
  /*     .. Parameters .. */
  /*     .. */
  /*     .. Local Scalars .. */
  /*     .. */
  /*     .. External Subroutines .. */
  /*     .. */
  /*     .. Intrinsic Functions .. */
  /*     .. */
  /*     .. Executable Statements .. */
  
  /*     Test the input arguments */
  
      /* Parameter adjustments */
      a_dim1 = *lda;
      a_offset = 1 + a_dim1;
      a -= a_offset;
      --tau;
      --work;
  
      /* Function Body */
      *info = 0;
      if (*m < 0) {
  	*info = -1;
      } else if (*n < 0) {
  	*info = -2;
      } else if (*lda < max(1,*m)) {
  	*info = -4;
      }
      if (*info != 0) {
  	i__1 = -(*info);
  	xerbla("CGEQR2", &i__1);
  	return 0;
      }
  
      k = min(*m,*n);
  
      i__1 = k;
      for (i__ = 1; i__ <= i__1; ++i__) {
  
  /*        Generate elementary reflector H(i) to annihilate A(i+1:m,i) */
  
  	i__2 = *m - i__ + 1;
  /* Computing MIN */
  	i__3 = i__ + 1;
  	clarfg(&i__2, &a[i__ + i__ * a_dim1], &a[min(i__3,*m) + i__ * a_dim1]
  		, &c__1, &tau[i__]);
  	if (i__ < *n) {
  
  /*           Apply H(i)' to A(i:m,i+1:n) from the left */
  
  	    i__2 = i__ + i__ * a_dim1;
  	    alpha.r = a[i__2].r, alpha.i = a[i__2].i;
  	    i__2 = i__ + i__ * a_dim1;
  	    a[i__2].r = 1.f, a[i__2].i = 0.f;
  	    i__2 = *m - i__ + 1;
  	    i__3 = *n - i__;
  	    r_cnjg(&q__1, &tau[i__]);
  	    clarf("Left", &i__2, &i__3, &a[i__ + i__ * a_dim1], &c__1, &q__1,
  		     &a[i__ + (i__ + 1) * a_dim1], lda, &work[1]);
  	    i__2 = i__ + i__ * a_dim1;
  	    a[i__2].r = alpha.r, a[i__2].i = alpha.i;
  	}
  /* L10: */
      }
      return 0;
  
  /*     End of CGEQR2 */
  
  } /* cgeqr2_ */
  
  