1 /* ./src_f77/sgeqlf.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static integer c_n1 = -1;
12 static integer c__3 = 3;
13 static integer c__2 = 2;
14
sgeqlf_(integer * m,integer * n,real * a,integer * lda,real * tau,real * work,integer * lwork,integer * info)15 /* Subroutine */ int sgeqlf_(integer *m, integer *n, real *a, integer *lda,
16 real *tau, real *work, integer *lwork, integer *info)
17 {
18 /* System generated locals */
19 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
20
21 /* Local variables */
22 static integer i__, k, ib, nb, ki, kk, mu, nu, nx, iws, nbmin, iinfo;
23 extern /* Subroutine */ int sgeql2_(integer *, integer *, real *, integer
24 *, real *, real *, integer *), slarfb_(char *, char *, char *,
25 char *, integer *, integer *, integer *, real *, integer *, real *
26 , integer *, real *, integer *, real *, integer *, ftnlen, ftnlen,
27 ftnlen, ftnlen), xerbla_(char *, integer *, ftnlen);
28 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
29 integer *, integer *, ftnlen, ftnlen);
30 extern /* Subroutine */ int slarft_(char *, char *, integer *, integer *,
31 real *, integer *, real *, real *, integer *, ftnlen, ftnlen);
32 static integer ldwork, lwkopt;
33 static logical lquery;
34
35
36 /* -- LAPACK routine (version 3.0) -- */
37 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
38 /* Courant Institute, Argonne National Lab, and Rice University */
39 /* June 30, 1999 */
40
41 /* .. Scalar Arguments .. */
42 /* .. */
43 /* .. Array Arguments .. */
44 /* .. */
45
46 /* Purpose */
47 /* ======= */
48
49 /* SGEQLF computes a QL factorization of a real M-by-N matrix A: */
50 /* A = Q * L. */
51
52 /* Arguments */
53 /* ========= */
54
55 /* M (input) INTEGER */
56 /* The number of rows of the matrix A. M >= 0. */
57
58 /* N (input) INTEGER */
59 /* The number of columns of the matrix A. N >= 0. */
60
61 /* A (input/output) REAL array, dimension (LDA,N) */
62 /* On entry, the M-by-N matrix A. */
63 /* On exit, */
64 /* if m >= n, the lower triangle of the subarray */
65 /* A(m-n+1:m,1:n) contains the N-by-N lower triangular matrix L; */
66 /* if m <= n, the elements on and below the (n-m)-th */
67 /* superdiagonal contain the M-by-N lower trapezoidal matrix L; */
68 /* the remaining elements, with the array TAU, represent the */
69 /* orthogonal matrix Q as a product of elementary reflectors */
70 /* (see Further Details). */
71
72 /* LDA (input) INTEGER */
73 /* The leading dimension of the array A. LDA >= max(1,M). */
74
75 /* TAU (output) REAL array, dimension (min(M,N)) */
76 /* The scalar factors of the elementary reflectors (see Further */
77 /* Details). */
78
79 /* WORK (workspace/output) REAL array, dimension (LWORK) */
80 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
81
82 /* LWORK (input) INTEGER */
83 /* The dimension of the array WORK. LWORK >= max(1,N). */
84 /* For optimum performance LWORK >= N*NB, where NB is the */
85 /* optimal blocksize. */
86
87 /* If LWORK = -1, then a workspace query is assumed; the routine */
88 /* only calculates the optimal size of the WORK array, returns */
89 /* this value as the first entry of the WORK array, and no error */
90 /* message related to LWORK is issued by XERBLA. */
91
92 /* INFO (output) INTEGER */
93 /* = 0: successful exit */
94 /* < 0: if INFO = -i, the i-th argument had an illegal value */
95
96 /* Further Details */
97 /* =============== */
98
99 /* The matrix Q is represented as a product of elementary reflectors */
100
101 /* Q = H(k) . . . H(2) H(1), where k = min(m,n). */
102
103 /* Each H(i) has the form */
104
105 /* H(i) = I - tau * v * v' */
106
107 /* where tau is a real scalar, and v is a real vector with */
108 /* v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in */
109 /* A(1:m-k+i-1,n-k+i), and tau in TAU(i). */
110
111 /* ===================================================================== */
112
113 /* .. Local Scalars .. */
114 /* .. */
115 /* .. External Subroutines .. */
116 /* .. */
117 /* .. Intrinsic Functions .. */
118 /* .. */
119 /* .. External Functions .. */
120 /* .. */
121 /* .. Executable Statements .. */
122
123 /* Test the input arguments */
124
125 /* Parameter adjustments */
126 a_dim1 = *lda;
127 a_offset = 1 + a_dim1;
128 a -= a_offset;
129 --tau;
130 --work;
131
132 /* Function Body */
133 *info = 0;
134 nb = ilaenv_(&c__1, "SGEQLF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)
135 1);
136 lwkopt = *n * nb;
137 work[1] = (real) lwkopt;
138 lquery = *lwork == -1;
139 if (*m < 0) {
140 *info = -1;
141 } else if (*n < 0) {
142 *info = -2;
143 } else if (*lda < max(1,*m)) {
144 *info = -4;
145 } else if (*lwork < max(1,*n) && ! lquery) {
146 *info = -7;
147 }
148 if (*info != 0) {
149 i__1 = -(*info);
150 xerbla_("SGEQLF", &i__1, (ftnlen)6);
151 return 0;
152 } else if (lquery) {
153 return 0;
154 }
155
156 /* Quick return if possible */
157
158 k = min(*m,*n);
159 if (k == 0) {
160 work[1] = 1.f;
161 return 0;
162 }
163
164 nbmin = 2;
165 nx = 1;
166 iws = *n;
167 if (nb > 1 && nb < k) {
168
169 /* Determine when to cross over from blocked to unblocked code. */
170
171 /* Computing MAX */
172 i__1 = 0, i__2 = ilaenv_(&c__3, "SGEQLF", " ", m, n, &c_n1, &c_n1, (
173 ftnlen)6, (ftnlen)1);
174 nx = max(i__1,i__2);
175 if (nx < k) {
176
177 /* Determine if workspace is large enough for blocked code. */
178
179 ldwork = *n;
180 iws = ldwork * nb;
181 if (*lwork < iws) {
182
183 /* Not enough workspace to use optimal NB: reduce NB and */
184 /* determine the minimum value of NB. */
185
186 nb = *lwork / ldwork;
187 /* Computing MAX */
188 i__1 = 2, i__2 = ilaenv_(&c__2, "SGEQLF", " ", m, n, &c_n1, &
189 c_n1, (ftnlen)6, (ftnlen)1);
190 nbmin = max(i__1,i__2);
191 }
192 }
193 }
194
195 if (nb >= nbmin && nb < k && nx < k) {
196
197 /* Use blocked code initially. */
198 /* The last kk columns are handled by the block method. */
199
200 ki = (k - nx - 1) / nb * nb;
201 /* Computing MIN */
202 i__1 = k, i__2 = ki + nb;
203 kk = min(i__1,i__2);
204
205 i__1 = k - kk + 1;
206 i__2 = -nb;
207 for (i__ = k - kk + ki + 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__
208 += i__2) {
209 /* Computing MIN */
210 i__3 = k - i__ + 1;
211 ib = min(i__3,nb);
212
213 /* Compute the QL factorization of the current block */
214 /* A(1:m-k+i+ib-1,n-k+i:n-k+i+ib-1) */
215
216 i__3 = *m - k + i__ + ib - 1;
217 sgeql2_(&i__3, &ib, &a[(*n - k + i__) * a_dim1 + 1], lda, &tau[
218 i__], &work[1], &iinfo);
219 if (*n - k + i__ > 1) {
220
221 /* Form the triangular factor of the block reflector */
222 /* H = H(i+ib-1) . . . H(i+1) H(i) */
223
224 i__3 = *m - k + i__ + ib - 1;
225 slarft_("Backward", "Columnwise", &i__3, &ib, &a[(*n - k +
226 i__) * a_dim1 + 1], lda, &tau[i__], &work[1], &ldwork,
227 (ftnlen)8, (ftnlen)10);
228
229 /* Apply H' to A(1:m-k+i+ib-1,1:n-k+i-1) from the left */
230
231 i__3 = *m - k + i__ + ib - 1;
232 i__4 = *n - k + i__ - 1;
233 slarfb_("Left", "Transpose", "Backward", "Columnwise", &i__3,
234 &i__4, &ib, &a[(*n - k + i__) * a_dim1 + 1], lda, &
235 work[1], &ldwork, &a[a_offset], lda, &work[ib + 1], &
236 ldwork, (ftnlen)4, (ftnlen)9, (ftnlen)8, (ftnlen)10);
237 }
238 /* L10: */
239 }
240 mu = *m - k + i__ + nb - 1;
241 nu = *n - k + i__ + nb - 1;
242 } else {
243 mu = *m;
244 nu = *n;
245 }
246
247 /* Use unblocked code to factor the last or only block */
248
249 if (mu > 0 && nu > 0) {
250 sgeql2_(&mu, &nu, &a[a_offset], lda, &tau[1], &work[1], &iinfo);
251 }
252
253 work[1] = (real) iws;
254 return 0;
255
256 /* End of SGEQLF */
257
258 } /* sgeqlf_ */
259
260