3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static integer c__2 = 2;
8 static integer c__65 = 65;
10 /* Subroutine */ int dormlq_(char *side, char *trans, integer *m, integer *n,
11 integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *
12 c__, integer *ldc, doublereal *work, integer *lwork, integer *info)
14 /* System generated locals */
16 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3[2], i__4,
20 /* Builtin functions */
21 /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
25 doublereal t[4160] /* was [65][64] */;
26 integer i1, i2, i3, ib, ic, jc, nb, mi, ni, nq, nw, iws;
28 extern logical lsame_(char *, char *);
30 extern /* Subroutine */ int dorml2_(char *, char *, integer *, integer *,
31 integer *, doublereal *, integer *, doublereal *, doublereal *,
32 integer *, doublereal *, integer *), dlarfb_(char
33 *, char *, char *, char *, integer *, integer *, integer *,
34 doublereal *, integer *, doublereal *, integer *, doublereal *,
35 integer *, doublereal *, integer *), dlarft_(char *, char *, integer *, integer *, doublereal
36 *, integer *, doublereal *, doublereal *, integer *), xerbla_(char *, integer *);
37 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
38 integer *, integer *);
46 /* -- LAPACK routine (version 3.1) -- */
47 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
50 /* .. Scalar Arguments .. */
52 /* .. Array Arguments .. */
58 /* DORMLQ overwrites the general real M-by-N matrix C with */
60 /* SIDE = 'L' SIDE = 'R' */
61 /* TRANS = 'N': Q * C C * Q */
62 /* TRANS = 'T': Q**T * C C * Q**T */
64 /* where Q is a real orthogonal matrix defined as the product of k */
65 /* elementary reflectors */
67 /* Q = H(k) . . . H(2) H(1) */
69 /* as returned by DGELQF. Q is of order M if SIDE = 'L' and of order N */
75 /* SIDE (input) CHARACTER*1 */
76 /* = 'L': apply Q or Q**T from the Left; */
77 /* = 'R': apply Q or Q**T from the Right. */
79 /* TRANS (input) CHARACTER*1 */
80 /* = 'N': No transpose, apply Q; */
81 /* = 'T': Transpose, apply Q**T. */
83 /* M (input) INTEGER */
84 /* The number of rows of the matrix C. M >= 0. */
86 /* N (input) INTEGER */
87 /* The number of columns of the matrix C. N >= 0. */
89 /* K (input) INTEGER */
90 /* The number of elementary reflectors whose product defines */
92 /* If SIDE = 'L', M >= K >= 0; */
93 /* if SIDE = 'R', N >= K >= 0. */
95 /* A (input) DOUBLE PRECISION array, dimension */
96 /* (LDA,M) if SIDE = 'L', */
97 /* (LDA,N) if SIDE = 'R' */
98 /* The i-th row must contain the vector which defines the */
99 /* elementary reflector H(i), for i = 1,2,...,k, as returned by */
100 /* DGELQF in the first k rows of its array argument A. */
101 /* A is modified by the routine but restored on exit. */
103 /* LDA (input) INTEGER */
104 /* The leading dimension of the array A. LDA >= max(1,K). */
106 /* TAU (input) DOUBLE PRECISION array, dimension (K) */
107 /* TAU(i) must contain the scalar factor of the elementary */
108 /* reflector H(i), as returned by DGELQF. */
110 /* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
111 /* On entry, the M-by-N matrix C. */
112 /* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. */
114 /* LDC (input) INTEGER */
115 /* The leading dimension of the array C. LDC >= max(1,M). */
117 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
118 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
120 /* LWORK (input) INTEGER */
121 /* The dimension of the array WORK. */
122 /* If SIDE = 'L', LWORK >= max(1,N); */
123 /* if SIDE = 'R', LWORK >= max(1,M). */
124 /* For optimum performance LWORK >= N*NB if SIDE = 'L', and */
125 /* LWORK >= M*NB if SIDE = 'R', where NB is the optimal */
128 /* If LWORK = -1, then a workspace query is assumed; the routine */
129 /* only calculates the optimal size of the WORK array, returns */
130 /* this value as the first entry of the WORK array, and no error */
131 /* message related to LWORK is issued by XERBLA. */
133 /* INFO (output) INTEGER */
134 /* = 0: successful exit */
135 /* < 0: if INFO = -i, the i-th argument had an illegal value */
137 /* ===================================================================== */
139 /* .. Parameters .. */
141 /* .. Local Scalars .. */
143 /* .. Local Arrays .. */
145 /* .. External Functions .. */
147 /* .. External Subroutines .. */
149 /* .. Intrinsic Functions .. */
151 /* .. Executable Statements .. */
153 /* Test the input arguments */
155 /* Parameter adjustments */
157 a_offset = 1 + a_dim1;
161 c_offset = 1 + c_dim1;
167 left = lsame_(side, "L");
168 notran = lsame_(trans, "N");
169 lquery = *lwork == -1;
171 /* NQ is the order of Q and NW is the minimum dimension of WORK */
180 if (! left && ! lsame_(side, "R")) {
182 } else if (! notran && ! lsame_(trans, "T")) {
188 } else if (*k < 0 || *k > nq) {
190 } else if (*lda < max(1,*k)) {
192 } else if (*ldc < max(1,*m)) {
194 } else if (*lwork < max(1,nw) && ! lquery) {
200 /* Determine the block size. NB may be at most NBMAX, where NBMAX */
201 /* is used to define the local array T. */
204 /* Writing concatenation */
205 i__3[0] = 1, a__1[0] = side;
206 i__3[1] = 1, a__1[1] = trans;
207 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
208 i__1 = 64, i__2 = ilaenv_(&c__1, "DORMLQ", ch__1, m, n, k, &c_n1);
210 lwkopt = max(1,nw) * nb;
211 work[1] = (doublereal) lwkopt;
216 xerbla_("DORMLQ", &i__1);
222 /* Quick return if possible */
224 if (*m == 0 || *n == 0 || *k == 0) {
231 if (nb > 1 && nb < *k) {
234 nb = *lwork / ldwork;
236 /* Writing concatenation */
237 i__3[0] = 1, a__1[0] = side;
238 i__3[1] = 1, a__1[1] = trans;
239 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
240 i__1 = 2, i__2 = ilaenv_(&c__2, "DORMLQ", ch__1, m, n, k, &c_n1);
241 nbmin = max(i__1,i__2);
247 if (nb < nbmin || nb >= *k) {
249 /* Use unblocked code */
251 dorml2_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
252 c_offset], ldc, &work[1], &iinfo);
255 /* Use blocked code */
257 if (left && notran || ! left && ! notran) {
262 i1 = (*k - 1) / nb * nb + 1;
276 *(unsigned char *)transt = 'T';
278 *(unsigned char *)transt = 'N';
283 for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
285 i__4 = nb, i__5 = *k - i__ + 1;
288 /* Form the triangular factor of the block reflector */
289 /* H = H(i) H(i+1) . . . H(i+ib-1) */
292 dlarft_("Forward", "Rowwise", &i__4, &ib, &a[i__ + i__ * a_dim1],
293 lda, &tau[i__], t, &c__65);
296 /* H or H' is applied to C(i:m,1:n) */
302 /* H or H' is applied to C(1:m,i:n) */
310 dlarfb_(side, transt, "Forward", "Rowwise", &mi, &ni, &ib, &a[i__
311 + i__ * a_dim1], lda, t, &c__65, &c__[ic + jc * c_dim1],
312 ldc, &work[1], &ldwork);
316 work[1] = (doublereal) lwkopt;