text: #include "clapack.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 dgeqrf_(integer* m, integer* n, doublereal* a, integer *
                               lda, doublereal* tau, doublereal* work, integer* lwork, integer* info) {
      /* System generated locals */
      integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
  
      /* Local variables */
      integer i__, k, ib, nb, nx, iws, nbmin, iinfo;
      extern /* Subroutine */ int dgeqr2_(integer*, integer*, doublereal*,
                                          integer*, doublereal*, doublereal*, integer*), dlarfb_(char*,
                                                  char*, char*, char*, integer*, integer*, integer*,
                                                  doublereal*, integer*, doublereal*, integer*, doublereal*,
                                                  integer*, doublereal*, integer*), dlarft_(char*, char*, integer*, integer*, doublereal
                                                          *, integer*, doublereal*, doublereal*, integer*), xerbla_(char*, integer*);
      extern integer ilaenv_(integer*, char*, char*, integer*, integer*,
                             integer*, integer*);
      integer ldwork, lwkopt;
      logical lquery;
  
  
      /*  -- LAPACK routine (version 3.1) -- */
      /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
      /*     November 2006 */
  
      /*     .. Scalar Arguments .. */
      /*     .. */
      /*     .. Array Arguments .. */
      /*     .. */
  
      /*  Purpose */
      /*  ======= */
  
      /*  DGEQRF computes a QR factorization of a real 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) DOUBLE PRECISION 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 orthogonal matrix Q as a */
      /*          product of min(m,n) elementary reflectors (see Further */
      /*          Details). */
  
      /*  LDA     (input) INTEGER */
      /*          The leading dimension of the array A.  LDA >= max(1,M). */
  
      /*  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N)) */
      /*          The scalar factors of the elementary reflectors (see Further */
      /*          Details). */
  
      /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
      /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
  
      /*  LWORK   (input) INTEGER */
      /*          The dimension of the array WORK.  LWORK >= max(1,N). */
      /*          For optimum performance LWORK >= N*NB, where NB is */
      /*          the optimal blocksize. */
  
      /*          If LWORK = -1, then a workspace query is assumed; the routine */
      /*          only calculates the optimal size of the WORK array, returns */
      /*          this value as the first entry of the WORK array, and no error */
      /*          message related to LWORK is issued by XERBLA. */
  
      /*  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 real scalar, and v is a real 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;
      nb = ilaenv_(&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1);
      lwkopt = *n * nb;
      work[1] = (doublereal) lwkopt;
      lquery = *lwork == -1;
      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) && ! lquery) {
          *info = -7;
      }
      if (*info != 0) {
          i__1 = -(*info);
          xerbla_("DGEQRF", &i__1);
          return 0;
      } else if (lquery) {
          return 0;
      }
  
      /*     Quick return if possible */
  
      k = min(*m, *n);
      if (k == 0) {
          work[1] = 1.;
          return 0;
      }
  
      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, "DGEQRF", " ", m, n, &c_n1, &c_n1);
          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, "DGEQRF", " ", m, n, &c_n1, &
                                           c_n1);
                  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;
              dgeqr2_(&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;
                  dlarft_("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;
                  dlarfb_("Left", "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;
          dgeqr2_(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
                  , &iinfo);
      }
  
      work[1] = (doublereal) iws;
      return 0;
  
      /*     End of DGEQRF */
  
  } /* dgeqrf_ */
  