1 /* ./src_f77/dorglq.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
dorglq_(integer * m,integer * n,integer * k,doublereal * a,integer * lda,doublereal * tau,doublereal * work,integer * lwork,integer * info)15 /* Subroutine */ int dorglq_(integer *m, integer *n, integer *k, doublereal *
16 a, integer *lda, doublereal *tau, doublereal *work, integer *lwork,
17 integer *info)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, i__1, i__2, i__3;
21
22 /* Local variables */
23 static integer i__, j, l, ib, nb, ki, kk, nx, iws, nbmin, iinfo;
24 extern /* Subroutine */ int dorgl2_(integer *, integer *, integer *,
25 doublereal *, integer *, doublereal *, doublereal *, integer *),
26 dlarfb_(char *, char *, char *, char *, integer *, integer *,
27 integer *, doublereal *, integer *, doublereal *, integer *,
28 doublereal *, integer *, doublereal *, integer *, ftnlen, ftnlen,
29 ftnlen, ftnlen), dlarft_(char *, char *, integer *, integer *,
30 doublereal *, integer *, doublereal *, doublereal *, integer *,
31 ftnlen, ftnlen), xerbla_(char *, integer *, ftnlen);
32 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
33 integer *, integer *, ftnlen, ftnlen);
34 static integer ldwork, lwkopt;
35 static logical lquery;
36
37
38 /* -- LAPACK routine (version 3.0) -- */
39 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
40 /* Courant Institute, Argonne National Lab, and Rice University */
41 /* June 30, 1999 */
42
43 /* .. Scalar Arguments .. */
44 /* .. */
45 /* .. Array Arguments .. */
46 /* .. */
47
48 /* Purpose */
49 /* ======= */
50
51 /* DORGLQ generates an M-by-N real matrix Q with orthonormal rows, */
52 /* which is defined as the first M rows of a product of K elementary */
53 /* reflectors of order N */
54
55 /* Q = H(k) . . . H(2) H(1) */
56
57 /* as returned by DGELQF. */
58
59 /* Arguments */
60 /* ========= */
61
62 /* M (input) INTEGER */
63 /* The number of rows of the matrix Q. M >= 0. */
64
65 /* N (input) INTEGER */
66 /* The number of columns of the matrix Q. N >= M. */
67
68 /* K (input) INTEGER */
69 /* The number of elementary reflectors whose product defines the */
70 /* matrix Q. M >= K >= 0. */
71
72 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
73 /* On entry, the i-th row must contain the vector which defines */
74 /* the elementary reflector H(i), for i = 1,2,...,k, as returned */
75 /* by DGELQF in the first k rows of its array argument A. */
76 /* On exit, the M-by-N matrix Q. */
77
78 /* LDA (input) INTEGER */
79 /* The first dimension of the array A. LDA >= max(1,M). */
80
81 /* TAU (input) DOUBLE PRECISION array, dimension (K) */
82 /* TAU(i) must contain the scalar factor of the elementary */
83 /* reflector H(i), as returned by DGELQF. */
84
85 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
86 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
87
88 /* LWORK (input) INTEGER */
89 /* The dimension of the array WORK. LWORK >= max(1,M). */
90 /* For optimum performance LWORK >= M*NB, where NB is */
91 /* the optimal blocksize. */
92
93 /* If LWORK = -1, then a workspace query is assumed; the routine */
94 /* only calculates the optimal size of the WORK array, returns */
95 /* this value as the first entry of the WORK array, and no error */
96 /* message related to LWORK is issued by XERBLA. */
97
98 /* INFO (output) INTEGER */
99 /* = 0: successful exit */
100 /* < 0: if INFO = -i, the i-th argument has an illegal value */
101
102 /* ===================================================================== */
103
104 /* .. Parameters .. */
105 /* .. */
106 /* .. Local Scalars .. */
107 /* .. */
108 /* .. External Subroutines .. */
109 /* .. */
110 /* .. Intrinsic Functions .. */
111 /* .. */
112 /* .. External Functions .. */
113 /* .. */
114 /* .. Executable Statements .. */
115
116 /* Test the input arguments */
117
118 /* Parameter adjustments */
119 a_dim1 = *lda;
120 a_offset = 1 + a_dim1;
121 a -= a_offset;
122 --tau;
123 --work;
124
125 /* Function Body */
126 *info = 0;
127 nb = ilaenv_(&c__1, "DORGLQ", " ", m, n, k, &c_n1, (ftnlen)6, (ftnlen)1);
128 lwkopt = max(1,*m) * nb;
129 work[1] = (doublereal) lwkopt;
130 lquery = *lwork == -1;
131 if (*m < 0) {
132 *info = -1;
133 } else if (*n < *m) {
134 *info = -2;
135 } else if (*k < 0 || *k > *m) {
136 *info = -3;
137 } else if (*lda < max(1,*m)) {
138 *info = -5;
139 } else if (*lwork < max(1,*m) && ! lquery) {
140 *info = -8;
141 }
142 if (*info != 0) {
143 i__1 = -(*info);
144 xerbla_("DORGLQ", &i__1, (ftnlen)6);
145 return 0;
146 } else if (lquery) {
147 return 0;
148 }
149
150 /* Quick return if possible */
151
152 if (*m <= 0) {
153 work[1] = 1.;
154 return 0;
155 }
156
157 nbmin = 2;
158 nx = 0;
159 iws = *m;
160 if (nb > 1 && nb < *k) {
161
162 /* Determine when to cross over from blocked to unblocked code. */
163
164 /* Computing MAX */
165 i__1 = 0, i__2 = ilaenv_(&c__3, "DORGLQ", " ", m, n, k, &c_n1, (
166 ftnlen)6, (ftnlen)1);
167 nx = max(i__1,i__2);
168 if (nx < *k) {
169
170 /* Determine if workspace is large enough for blocked code. */
171
172 ldwork = *m;
173 iws = ldwork * nb;
174 if (*lwork < iws) {
175
176 /* Not enough workspace to use optimal NB: reduce NB and */
177 /* determine the minimum value of NB. */
178
179 nb = *lwork / ldwork;
180 /* Computing MAX */
181 i__1 = 2, i__2 = ilaenv_(&c__2, "DORGLQ", " ", m, n, k, &c_n1,
182 (ftnlen)6, (ftnlen)1);
183 nbmin = max(i__1,i__2);
184 }
185 }
186 }
187
188 if (nb >= nbmin && nb < *k && nx < *k) {
189
190 /* Use blocked code after the last block. */
191 /* The first kk rows are handled by the block method. */
192
193 ki = (*k - nx - 1) / nb * nb;
194 /* Computing MIN */
195 i__1 = *k, i__2 = ki + nb;
196 kk = min(i__1,i__2);
197
198 /* Set A(kk+1:m,1:kk) to zero. */
199
200 i__1 = kk;
201 for (j = 1; j <= i__1; ++j) {
202 i__2 = *m;
203 for (i__ = kk + 1; i__ <= i__2; ++i__) {
204 a[i__ + j * a_dim1] = 0.;
205 /* L10: */
206 }
207 /* L20: */
208 }
209 } else {
210 kk = 0;
211 }
212
213 /* Use unblocked code for the last or only block. */
214
215 if (kk < *m) {
216 i__1 = *m - kk;
217 i__2 = *n - kk;
218 i__3 = *k - kk;
219 dorgl2_(&i__1, &i__2, &i__3, &a[kk + 1 + (kk + 1) * a_dim1], lda, &
220 tau[kk + 1], &work[1], &iinfo);
221 }
222
223 if (kk > 0) {
224
225 /* Use blocked code */
226
227 i__1 = -nb;
228 for (i__ = ki + 1; i__1 < 0 ? i__ >= 1 : i__ <= 1; i__ += i__1) {
229 /* Computing MIN */
230 i__2 = nb, i__3 = *k - i__ + 1;
231 ib = min(i__2,i__3);
232 if (i__ + ib <= *m) {
233
234 /* Form the triangular factor of the block reflector */
235 /* H = H(i) H(i+1) . . . H(i+ib-1) */
236
237 i__2 = *n - i__ + 1;
238 dlarft_("Forward", "Rowwise", &i__2, &ib, &a[i__ + i__ *
239 a_dim1], lda, &tau[i__], &work[1], &ldwork, (ftnlen)7,
240 (ftnlen)7);
241
242 /* Apply H' to A(i+ib:m,i:n) from the right */
243
244 i__2 = *m - i__ - ib + 1;
245 i__3 = *n - i__ + 1;
246 dlarfb_("Right", "Transpose", "Forward", "Rowwise", &i__2, &
247 i__3, &ib, &a[i__ + i__ * a_dim1], lda, &work[1], &
248 ldwork, &a[i__ + ib + i__ * a_dim1], lda, &work[ib +
249 1], &ldwork, (ftnlen)5, (ftnlen)9, (ftnlen)7, (ftnlen)
250 7);
251 }
252
253 /* Apply H' to columns i:n of current block */
254
255 i__2 = *n - i__ + 1;
256 dorgl2_(&ib, &i__2, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &
257 work[1], &iinfo);
258
259 /* Set columns 1:i-1 of current block to zero */
260
261 i__2 = i__ - 1;
262 for (j = 1; j <= i__2; ++j) {
263 i__3 = i__ + ib - 1;
264 for (l = i__; l <= i__3; ++l) {
265 a[l + j * a_dim1] = 0.;
266 /* L30: */
267 }
268 /* L40: */
269 }
270 /* L50: */
271 }
272 }
273
274 work[1] = (doublereal) iws;
275 return 0;
276
277 /* End of DORGLQ */
278
279 } /* dorglq_ */
280
281