3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static integer c__3 = 3;
8 static integer c__2 = 2;
10 /* Subroutine */ int dgelqf_(integer *m, integer *n, doublereal *a, integer *
11 lda, doublereal *tau, doublereal *work, integer *lwork, integer *info)
13 /* System generated locals */
14 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
17 integer i__, k, ib, nb, nx, iws, nbmin, iinfo;
18 extern /* Subroutine */ int dgelq2_(integer *, integer *, doublereal *,
19 integer *, doublereal *, doublereal *, integer *), dlarfb_(char *,
20 char *, char *, char *, integer *, integer *, integer *,
21 doublereal *, integer *, doublereal *, integer *, doublereal *,
22 integer *, doublereal *, integer *), dlarft_(char *, char *, integer *, integer *, doublereal
23 *, integer *, doublereal *, doublereal *, integer *), xerbla_(char *, integer *);
24 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
25 integer *, integer *);
26 integer ldwork, lwkopt;
30 /* -- LAPACK routine (version 3.1) -- */
31 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
34 /* .. Scalar Arguments .. */
36 /* .. Array Arguments .. */
42 /* DGELQF computes an LQ factorization of a real M-by-N matrix A: */
48 /* M (input) INTEGER */
49 /* The number of rows of the matrix A. M >= 0. */
51 /* N (input) INTEGER */
52 /* The number of columns of the matrix A. N >= 0. */
54 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
55 /* On entry, the M-by-N matrix A. */
56 /* On exit, the elements on and below the diagonal of the array */
57 /* contain the m-by-min(m,n) lower trapezoidal matrix L (L is */
58 /* lower triangular if m <= n); the elements above the diagonal, */
59 /* with the array TAU, represent the orthogonal matrix Q as a */
60 /* product of elementary reflectors (see Further Details). */
62 /* LDA (input) INTEGER */
63 /* The leading dimension of the array A. LDA >= max(1,M). */
65 /* TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) */
66 /* The scalar factors of the elementary reflectors (see Further */
69 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
70 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
72 /* LWORK (input) INTEGER */
73 /* The dimension of the array WORK. LWORK >= max(1,M). */
74 /* For optimum performance LWORK >= M*NB, where NB is the */
75 /* optimal blocksize. */
77 /* If LWORK = -1, then a workspace query is assumed; the routine */
78 /* only calculates the optimal size of the WORK array, returns */
79 /* this value as the first entry of the WORK array, and no error */
80 /* message related to LWORK is issued by XERBLA. */
82 /* INFO (output) INTEGER */
83 /* = 0: successful exit */
84 /* < 0: if INFO = -i, the i-th argument had an illegal value */
89 /* The matrix Q is represented as a product of elementary reflectors */
91 /* Q = H(k) . . . H(2) H(1), where k = min(m,n). */
93 /* Each H(i) has the form */
95 /* H(i) = I - tau * v * v' */
97 /* where tau is a real scalar, and v is a real vector with */
98 /* v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), */
99 /* and tau in TAU(i). */
101 /* ===================================================================== */
103 /* .. Local Scalars .. */
105 /* .. External Subroutines .. */
107 /* .. Intrinsic Functions .. */
109 /* .. External Functions .. */
111 /* .. Executable Statements .. */
113 /* Test the input arguments */
115 /* Parameter adjustments */
117 a_offset = 1 + a_dim1;
124 nb = ilaenv_(&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1);
126 work[1] = (doublereal) lwkopt;
127 lquery = *lwork == -1;
132 } else if (*lda < max(1,*m)) {
134 } else if (*lwork < max(1,*m) && ! lquery) {
139 xerbla_("DGELQF", &i__1);
145 /* Quick return if possible */
156 if (nb > 1 && nb < k) {
158 /* Determine when to cross over from blocked to unblocked code. */
161 i__1 = 0, i__2 = ilaenv_(&c__3, "DGELQF", " ", m, n, &c_n1, &c_n1);
165 /* Determine if workspace is large enough for blocked code. */
171 /* Not enough workspace to use optimal NB: reduce NB and */
172 /* determine the minimum value of NB. */
174 nb = *lwork / ldwork;
176 i__1 = 2, i__2 = ilaenv_(&c__2, "DGELQF", " ", m, n, &c_n1, &
178 nbmin = max(i__1,i__2);
183 if (nb >= nbmin && nb < k && nx < k) {
185 /* Use blocked code initially */
189 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
194 /* Compute the LQ factorization of the current block */
195 /* A(i:i+ib-1,i:n) */
198 dgelq2_(&ib, &i__3, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
200 if (i__ + ib <= *m) {
202 /* Form the triangular factor of the block reflector */
203 /* H = H(i) H(i+1) . . . H(i+ib-1) */
206 dlarft_("Forward", "Rowwise", &i__3, &ib, &a[i__ + i__ *
207 a_dim1], lda, &tau[i__], &work[1], &ldwork);
209 /* Apply H to A(i+ib:m,i:n) from the right */
211 i__3 = *m - i__ - ib + 1;
213 dlarfb_("Right", "No transpose", "Forward", "Rowwise", &i__3,
214 &i__4, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
215 ldwork, &a[i__ + ib + i__ * a_dim1], lda, &work[ib +
224 /* Use unblocked code to factor the last or only block. */
229 dgelq2_(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
233 work[1] = (doublereal) iws;