1 /* ./src_f77/sgelqf.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
sgelqf_(integer * m,integer * n,real * a,integer * lda,real * tau,real * work,integer * lwork,integer * info)15 /* Subroutine */ int sgelqf_(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, nx, iws, nbmin, iinfo;
23 extern /* Subroutine */ int sgelq2_(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 /* SGELQF computes an LQ factorization of a real M-by-N matrix A: */
50 /* A = L * Q. */
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, the elements on and below the diagonal of the array */
64 /* contain the m-by-min(m,n) lower trapezoidal matrix L (L is */
65 /* lower triangular if m <= n); the elements above the diagonal, */
66 /* with the array TAU, represent the orthogonal matrix Q as a */
67 /* product of elementary reflectors (see Further Details). */
68
69 /* LDA (input) INTEGER */
70 /* The leading dimension of the array A. LDA >= max(1,M). */
71
72 /* TAU (output) REAL array, dimension (min(M,N)) */
73 /* The scalar factors of the elementary reflectors (see Further */
74 /* Details). */
75
76 /* WORK (workspace/output) REAL array, dimension (LWORK) */
77 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
78
79 /* LWORK (input) INTEGER */
80 /* The dimension of the array WORK. LWORK >= max(1,M). */
81 /* For optimum performance LWORK >= M*NB, where NB is the */
82 /* optimal blocksize. */
83
84 /* If LWORK = -1, then a workspace query is assumed; the routine */
85 /* only calculates the optimal size of the WORK array, returns */
86 /* this value as the first entry of the WORK array, and no error */
87 /* message related to LWORK is issued by XERBLA. */
88
89 /* INFO (output) INTEGER */
90 /* = 0: successful exit */
91 /* < 0: if INFO = -i, the i-th argument had an illegal value */
92
93 /* Further Details */
94 /* =============== */
95
96 /* The matrix Q is represented as a product of elementary reflectors */
97
98 /* Q = H(k) . . . H(2) H(1), where k = min(m,n). */
99
100 /* Each H(i) has the form */
101
102 /* H(i) = I - tau * v * v' */
103
104 /* where tau is a real scalar, and v is a real vector with */
105 /* v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i,i+1:n), */
106 /* and tau in TAU(i). */
107
108 /* ===================================================================== */
109
110 /* .. Local Scalars .. */
111 /* .. */
112 /* .. External Subroutines .. */
113 /* .. */
114 /* .. Intrinsic Functions .. */
115 /* .. */
116 /* .. External Functions .. */
117 /* .. */
118 /* .. Executable Statements .. */
119
120 /* Test the input arguments */
121
122 /* Parameter adjustments */
123 a_dim1 = *lda;
124 a_offset = 1 + a_dim1;
125 a -= a_offset;
126 --tau;
127 --work;
128
129 /* Function Body */
130 *info = 0;
131 nb = ilaenv_(&c__1, "SGELQF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)
132 1);
133 lwkopt = *m * nb;
134 work[1] = (real) lwkopt;
135 lquery = *lwork == -1;
136 if (*m < 0) {
137 *info = -1;
138 } else if (*n < 0) {
139 *info = -2;
140 } else if (*lda < max(1,*m)) {
141 *info = -4;
142 } else if (*lwork < max(1,*m) && ! lquery) {
143 *info = -7;
144 }
145 if (*info != 0) {
146 i__1 = -(*info);
147 xerbla_("SGELQF", &i__1, (ftnlen)6);
148 return 0;
149 } else if (lquery) {
150 return 0;
151 }
152
153 /* Quick return if possible */
154
155 k = min(*m,*n);
156 if (k == 0) {
157 work[1] = 1.f;
158 return 0;
159 }
160
161 nbmin = 2;
162 nx = 0;
163 iws = *m;
164 if (nb > 1 && nb < k) {
165
166 /* Determine when to cross over from blocked to unblocked code. */
167
168 /* Computing MAX */
169 i__1 = 0, i__2 = ilaenv_(&c__3, "SGELQF", " ", m, n, &c_n1, &c_n1, (
170 ftnlen)6, (ftnlen)1);
171 nx = max(i__1,i__2);
172 if (nx < k) {
173
174 /* Determine if workspace is large enough for blocked code. */
175
176 ldwork = *m;
177 iws = ldwork * nb;
178 if (*lwork < iws) {
179
180 /* Not enough workspace to use optimal NB: reduce NB and */
181 /* determine the minimum value of NB. */
182
183 nb = *lwork / ldwork;
184 /* Computing MAX */
185 i__1 = 2, i__2 = ilaenv_(&c__2, "SGELQF", " ", m, n, &c_n1, &
186 c_n1, (ftnlen)6, (ftnlen)1);
187 nbmin = max(i__1,i__2);
188 }
189 }
190 }
191
192 if (nb >= nbmin && nb < k && nx < k) {
193
194 /* Use blocked code initially */
195
196 i__1 = k - nx;
197 i__2 = nb;
198 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
199 /* Computing MIN */
200 i__3 = k - i__ + 1;
201 ib = min(i__3,nb);
202
203 /* Compute the LQ factorization of the current block */
204 /* A(i:i+ib-1,i:n) */
205
206 i__3 = *n - i__ + 1;
207 sgelq2_(&ib, &i__3, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
208 1], &iinfo);
209 if (i__ + ib <= *m) {
210
211 /* Form the triangular factor of the block reflector */
212 /* H = H(i) H(i+1) . . . H(i+ib-1) */
213
214 i__3 = *n - i__ + 1;
215 slarft_("Forward", "Rowwise", &i__3, &ib, &a[i__ + i__ *
216 a_dim1], lda, &tau[i__], &work[1], &ldwork, (ftnlen)7,
217 (ftnlen)7);
218
219 /* Apply H to A(i+ib:m,i:n) from the right */
220
221 i__3 = *m - i__ - ib + 1;
222 i__4 = *n - i__ + 1;
223 slarfb_("Right", "No transpose", "Forward", "Rowwise", &i__3,
224 &i__4, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
225 ldwork, &a[i__ + ib + i__ * a_dim1], lda, &work[ib +
226 1], &ldwork, (ftnlen)5, (ftnlen)12, (ftnlen)7, (
227 ftnlen)7);
228 }
229 /* L10: */
230 }
231 } else {
232 i__ = 1;
233 }
234
235 /* Use unblocked code to factor the last or only block. */
236
237 if (i__ <= k) {
238 i__2 = *m - i__ + 1;
239 i__1 = *n - i__ + 1;
240 sgelq2_(&i__2, &i__1, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[1]
241 , &iinfo);
242 }
243
244 work[1] = (real) iws;
245 return 0;
246
247 /* End of SGELQF */
248
249 } /* sgelqf_ */
250
251