3 /* Subroutine */ int dorml2_(char *side, char *trans, integer *m, integer *n,
4 integer *k, doublereal *a, integer *lda, doublereal *tau, doublereal *
5 c__, integer *ldc, doublereal *work, integer *info)
7 /* System generated locals */
8 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2;
11 integer i__, i1, i2, i3, ic, jc, mi, ni, nq;
14 extern /* Subroutine */ int dlarf_(char *, integer *, integer *,
15 doublereal *, integer *, doublereal *, doublereal *, integer *,
17 extern logical lsame_(char *, char *);
18 extern /* Subroutine */ int xerbla_(char *, integer *);
22 /* -- LAPACK routine (version 3.1) -- */
23 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
26 /* .. Scalar Arguments .. */
28 /* .. Array Arguments .. */
34 /* DORML2 overwrites the general real m by n matrix C with */
36 /* Q * C if SIDE = 'L' and TRANS = 'N', or */
38 /* Q'* C if SIDE = 'L' and TRANS = 'T', or */
40 /* C * Q if SIDE = 'R' and TRANS = 'N', or */
42 /* C * Q' if SIDE = 'R' and TRANS = 'T', */
44 /* where Q is a real orthogonal matrix defined as the product of k */
45 /* elementary reflectors */
47 /* Q = H(k) . . . H(2) H(1) */
49 /* as returned by DGELQF. Q is of order m if SIDE = 'L' and of order n */
55 /* SIDE (input) CHARACTER*1 */
56 /* = 'L': apply Q or Q' from the Left */
57 /* = 'R': apply Q or Q' from the Right */
59 /* TRANS (input) CHARACTER*1 */
60 /* = 'N': apply Q (No transpose) */
61 /* = 'T': apply Q' (Transpose) */
63 /* M (input) INTEGER */
64 /* The number of rows of the matrix C. M >= 0. */
66 /* N (input) INTEGER */
67 /* The number of columns of the matrix C. N >= 0. */
69 /* K (input) INTEGER */
70 /* The number of elementary reflectors whose product defines */
72 /* If SIDE = 'L', M >= K >= 0; */
73 /* if SIDE = 'R', N >= K >= 0. */
75 /* A (input) DOUBLE PRECISION array, dimension */
76 /* (LDA,M) if SIDE = 'L', */
77 /* (LDA,N) if SIDE = 'R' */
78 /* The i-th row must contain the vector which defines the */
79 /* elementary reflector H(i), for i = 1,2,...,k, as returned by */
80 /* DGELQF in the first k rows of its array argument A. */
81 /* A is modified by the routine but restored on exit. */
83 /* LDA (input) INTEGER */
84 /* The leading dimension of the array A. LDA >= max(1,K). */
86 /* TAU (input) DOUBLE PRECISION array, dimension (K) */
87 /* TAU(i) must contain the scalar factor of the elementary */
88 /* reflector H(i), as returned by DGELQF. */
90 /* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) */
91 /* On entry, the m by n matrix C. */
92 /* On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. */
94 /* LDC (input) INTEGER */
95 /* The leading dimension of the array C. LDC >= max(1,M). */
97 /* WORK (workspace) DOUBLE PRECISION array, dimension */
98 /* (N) if SIDE = 'L', */
99 /* (M) if SIDE = 'R' */
101 /* INFO (output) INTEGER */
102 /* = 0: successful exit */
103 /* < 0: if INFO = -i, the i-th argument had an illegal value */
105 /* ===================================================================== */
107 /* .. Parameters .. */
109 /* .. Local Scalars .. */
111 /* .. External Functions .. */
113 /* .. External Subroutines .. */
115 /* .. Intrinsic Functions .. */
117 /* .. Executable Statements .. */
119 /* Test the input arguments */
121 /* Parameter adjustments */
123 a_offset = 1 + a_dim1;
127 c_offset = 1 + c_dim1;
133 left = lsame_(side, "L");
134 notran = lsame_(trans, "N");
136 /* NQ is the order of Q */
143 if (! left && ! lsame_(side, "R")) {
145 } else if (! notran && ! lsame_(trans, "T")) {
151 } else if (*k < 0 || *k > nq) {
153 } else if (*lda < max(1,*k)) {
155 } else if (*ldc < max(1,*m)) {
160 xerbla_("DORML2", &i__1);
164 /* Quick return if possible */
166 if (*m == 0 || *n == 0 || *k == 0) {
170 if (left && notran || ! left && ! notran) {
190 for (i__ = i1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
193 /* H(i) is applied to C(i:m,1:n) */
199 /* H(i) is applied to C(1:m,i:n) */
207 aii = a[i__ + i__ * a_dim1];
208 a[i__ + i__ * a_dim1] = 1.;
209 dlarf_(side, &mi, &ni, &a[i__ + i__ * a_dim1], lda, &tau[i__], &c__[
210 ic + jc * c_dim1], ldc, &work[1]);
211 a[i__ + i__ * a_dim1] = aii;