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 dormqr_(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 dorm2r_(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 *);
40 integer ldwork, lwkopt;
44 /* -- LAPACK routine (version 3.1) -- */
45 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
48 /* .. Scalar Arguments .. */
50 /* .. Array Arguments .. */
56 /* DORMQR overwrites the general real M-by-N matrix C with */
58 /* SIDE = 'L' SIDE = 'R' */
59 /* TRANS = 'N': Q * C C * Q */
60 /* TRANS = 'T': Q**T * C C * Q**T */
62 /* where Q is a real orthogonal matrix defined as the product of k */
63 /* elementary reflectors */
65 /* Q = H(1) H(2) . . . H(k) */
67 /* as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N */
73 /* SIDE (input) CHARACTER*1 */
74 /* = 'L': apply Q or Q**T from the Left; */
75 /* = 'R': apply Q or Q**T from the Right. */
77 /* TRANS (input) CHARACTER*1 */
78 /* = 'N': No transpose, apply Q; */
79 /* = 'T': Transpose, apply Q**T. */
81 /* M (input) INTEGER */
82 /* The number of rows of the matrix C. M >= 0. */
84 /* N (input) INTEGER */
85 /* The number of columns of the matrix C. N >= 0. */
87 /* K (input) INTEGER */
88 /* The number of elementary reflectors whose product defines */
90 /* If SIDE = 'L', M >= K >= 0; */
91 /* if SIDE = 'R', N >= K >= 0. */
93 /* A (input) DOUBLE PRECISION array, dimension (LDA,K) */
94 /* The i-th column must contain the vector which defines the */
95 /* elementary reflector H(i), for i = 1,2,...,k, as returned by */
96 /* DGEQRF in the first k columns of its array argument A. */
97 /* A is modified by the routine but restored on exit. */
99 /* LDA (input) INTEGER */
100 /* The leading dimension of the array A. */
101 /* If SIDE = 'L', LDA >= max(1,M); */
102 /* if SIDE = 'R', LDA >= max(1,N). */
104 /* TAU (input) DOUBLE PRECISION array, dimension (K) */
105 /* TAU(i) must contain the scalar factor of the elementary */
106 /* reflector H(i), as returned by DGEQRF. */
108 /* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
109 /* On entry, the M-by-N matrix C. */
110 /* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. */
112 /* LDC (input) INTEGER */
113 /* The leading dimension of the array C. LDC >= max(1,M). */
115 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
116 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
118 /* LWORK (input) INTEGER */
119 /* The dimension of the array WORK. */
120 /* If SIDE = 'L', LWORK >= max(1,N); */
121 /* if SIDE = 'R', LWORK >= max(1,M). */
122 /* For optimum performance LWORK >= N*NB if SIDE = 'L', and */
123 /* LWORK >= M*NB if SIDE = 'R', where NB is the optimal */
126 /* If LWORK = -1, then a workspace query is assumed; the routine */
127 /* only calculates the optimal size of the WORK array, returns */
128 /* this value as the first entry of the WORK array, and no error */
129 /* message related to LWORK is issued by XERBLA. */
131 /* INFO (output) INTEGER */
132 /* = 0: successful exit */
133 /* < 0: if INFO = -i, the i-th argument had an illegal value */
135 /* ===================================================================== */
137 /* .. Parameters .. */
139 /* .. Local Scalars .. */
141 /* .. Local Arrays .. */
143 /* .. External Functions .. */
145 /* .. External Subroutines .. */
147 /* .. Intrinsic Functions .. */
149 /* .. Executable Statements .. */
151 /* Test the input arguments */
153 /* Parameter adjustments */
155 a_offset = 1 + a_dim1;
159 c_offset = 1 + c_dim1;
165 left = lsame_(side, "L");
166 notran = lsame_(trans, "N");
167 lquery = *lwork == -1;
169 /* NQ is the order of Q and NW is the minimum dimension of WORK */
178 if (! left && ! lsame_(side, "R")) {
180 } else if (! notran && ! lsame_(trans, "T")) {
186 } else if (*k < 0 || *k > nq) {
188 } else if (*lda < max(1,nq)) {
190 } else if (*ldc < max(1,*m)) {
192 } else if (*lwork < max(1,nw) && ! lquery) {
198 /* Determine the block size. NB may be at most NBMAX, where NBMAX */
199 /* is used to define the local array T. */
202 /* Writing concatenation */
203 i__3[0] = 1, a__1[0] = side;
204 i__3[1] = 1, a__1[1] = trans;
205 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
206 i__1 = 64, i__2 = ilaenv_(&c__1, "DORMQR", ch__1, m, n, k, &c_n1);
208 lwkopt = max(1,nw) * nb;
209 work[1] = (doublereal) lwkopt;
214 xerbla_("DORMQR", &i__1);
220 /* Quick return if possible */
222 if (*m == 0 || *n == 0 || *k == 0) {
229 if (nb > 1 && nb < *k) {
232 nb = *lwork / ldwork;
234 /* Writing concatenation */
235 i__3[0] = 1, a__1[0] = side;
236 i__3[1] = 1, a__1[1] = trans;
237 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
238 i__1 = 2, i__2 = ilaenv_(&c__2, "DORMQR", ch__1, m, n, k, &c_n1);
239 nbmin = max(i__1,i__2);
245 if (nb < nbmin || nb >= *k) {
247 /* Use unblocked code */
249 dorm2r_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
250 c_offset], ldc, &work[1], &iinfo);
253 /* Use blocked code */
255 if (left && ! notran || ! left && notran) {
260 i1 = (*k - 1) / nb * nb + 1;
275 for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
277 i__4 = nb, i__5 = *k - i__ + 1;
280 /* Form the triangular factor of the block reflector */
281 /* H = H(i) H(i+1) . . . H(i+ib-1) */
284 dlarft_("Forward", "Columnwise", &i__4, &ib, &a[i__ + i__ *
285 a_dim1], lda, &tau[i__], t, &c__65)
289 /* H or H' is applied to C(i:m,1:n) */
295 /* H or H' is applied to C(1:m,i:n) */
303 dlarfb_(side, trans, "Forward", "Columnwise", &mi, &ni, &ib, &a[
304 i__ + i__ * a_dim1], lda, t, &c__65, &c__[ic + jc *
305 c_dim1], ldc, &work[1], &ldwork);
309 work[1] = (doublereal) lwkopt;