3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static integer c__2 = 2;
9 /* Subroutine */ int sormbr_(char *vect, char *side, char *trans, integer *m,
10 integer *n, integer *k, real *a, integer *lda, real *tau, real *c__,
11 integer *ldc, real *work, integer *lwork, integer *info)
13 /* System generated locals */
15 integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3[2];
18 /* Builtin functions */
19 /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
22 integer i1, i2, nb, mi, ni, nq, nw;
24 extern logical lsame_(char *, char *);
26 extern /* Subroutine */ int xerbla_(char *, integer *);
27 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
28 integer *, integer *);
29 logical notran, applyq;
31 extern /* Subroutine */ int sormlq_(char *, char *, integer *, integer *,
32 integer *, real *, integer *, real *, real *, integer *, real *,
33 integer *, integer *);
36 extern /* Subroutine */ int sormqr_(char *, char *, integer *, integer *,
37 integer *, real *, integer *, real *, real *, integer *, real *,
38 integer *, integer *);
41 /* -- LAPACK routine (version 3.1) -- */
42 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
45 /* .. Scalar Arguments .. */
47 /* .. Array Arguments .. */
53 /* If VECT = 'Q', SORMBR overwrites the general real M-by-N matrix C */
55 /* SIDE = 'L' SIDE = 'R' */
56 /* TRANS = 'N': Q * C C * Q */
57 /* TRANS = 'T': Q**T * C C * Q**T */
59 /* If VECT = 'P', SORMBR overwrites the general real M-by-N matrix C */
61 /* SIDE = 'L' SIDE = 'R' */
62 /* TRANS = 'N': P * C C * P */
63 /* TRANS = 'T': P**T * C C * P**T */
65 /* Here Q and P**T are the orthogonal matrices determined by SGEBRD when */
66 /* reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and */
67 /* P**T are defined as products of elementary reflectors H(i) and G(i) */
70 /* Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the */
71 /* order of the orthogonal matrix Q or P**T that is applied. */
73 /* If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: */
74 /* if nq >= k, Q = H(1) H(2) . . . H(k); */
75 /* if nq < k, Q = H(1) H(2) . . . H(nq-1). */
77 /* If VECT = 'P', A is assumed to have been a K-by-NQ matrix: */
78 /* if k < nq, P = G(1) G(2) . . . G(k); */
79 /* if k >= nq, P = G(1) G(2) . . . G(nq-1). */
84 /* VECT (input) CHARACTER*1 */
85 /* = 'Q': apply Q or Q**T; */
86 /* = 'P': apply P or P**T. */
88 /* SIDE (input) CHARACTER*1 */
89 /* = 'L': apply Q, Q**T, P or P**T from the Left; */
90 /* = 'R': apply Q, Q**T, P or P**T from the Right. */
92 /* TRANS (input) CHARACTER*1 */
93 /* = 'N': No transpose, apply Q or P; */
94 /* = 'T': Transpose, apply Q**T or P**T. */
96 /* M (input) INTEGER */
97 /* The number of rows of the matrix C. M >= 0. */
99 /* N (input) INTEGER */
100 /* The number of columns of the matrix C. N >= 0. */
102 /* K (input) INTEGER */
103 /* If VECT = 'Q', the number of columns in the original */
104 /* matrix reduced by SGEBRD. */
105 /* If VECT = 'P', the number of rows in the original */
106 /* matrix reduced by SGEBRD. */
109 /* A (input) REAL array, dimension */
110 /* (LDA,min(nq,K)) if VECT = 'Q' */
111 /* (LDA,nq) if VECT = 'P' */
112 /* The vectors which define the elementary reflectors H(i) and */
113 /* G(i), whose products determine the matrices Q and P, as */
114 /* returned by SGEBRD. */
116 /* LDA (input) INTEGER */
117 /* The leading dimension of the array A. */
118 /* If VECT = 'Q', LDA >= max(1,nq); */
119 /* if VECT = 'P', LDA >= max(1,min(nq,K)). */
121 /* TAU (input) REAL array, dimension (min(nq,K)) */
122 /* TAU(i) must contain the scalar factor of the elementary */
123 /* reflector H(i) or G(i) which determines Q or P, as returned */
124 /* by SGEBRD in the array argument TAUQ or TAUP. */
126 /* C (input/output) REAL array, dimension (LDC,N) */
127 /* On entry, the M-by-N matrix C. */
128 /* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q */
129 /* or P*C or P**T*C or C*P or C*P**T. */
131 /* LDC (input) INTEGER */
132 /* The leading dimension of the array C. LDC >= max(1,M). */
134 /* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
135 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
137 /* LWORK (input) INTEGER */
138 /* The dimension of the array WORK. */
139 /* If SIDE = 'L', LWORK >= max(1,N); */
140 /* if SIDE = 'R', LWORK >= max(1,M). */
141 /* For optimum performance LWORK >= N*NB if SIDE = 'L', and */
142 /* LWORK >= M*NB if SIDE = 'R', where NB is the optimal */
145 /* If LWORK = -1, then a workspace query is assumed; the routine */
146 /* only calculates the optimal size of the WORK array, returns */
147 /* this value as the first entry of the WORK array, and no error */
148 /* message related to LWORK is issued by XERBLA. */
150 /* INFO (output) INTEGER */
151 /* = 0: successful exit */
152 /* < 0: if INFO = -i, the i-th argument had an illegal value */
154 /* ===================================================================== */
156 /* .. Local Scalars .. */
158 /* .. External Functions .. */
160 /* .. External Subroutines .. */
162 /* .. Intrinsic Functions .. */
164 /* .. Executable Statements .. */
166 /* Test the input arguments */
168 /* Parameter adjustments */
170 a_offset = 1 + a_dim1;
174 c_offset = 1 + c_dim1;
180 applyq = lsame_(vect, "Q");
181 left = lsame_(side, "L");
182 notran = lsame_(trans, "N");
183 lquery = *lwork == -1;
185 /* NQ is the order of Q or P and NW is the minimum dimension of WORK */
194 if (! applyq && ! lsame_(vect, "P")) {
196 } else if (! left && ! lsame_(side, "R")) {
198 } else if (! notran && ! lsame_(trans, "T")) {
206 } else /* if(complicated condition) */ {
208 i__1 = 1, i__2 = min(nq,*k);
209 if (applyq && *lda < max(1,nq) || ! applyq && *lda < max(i__1,i__2)) {
211 } else if (*ldc < max(1,*m)) {
213 } else if (*lwork < max(1,nw) && ! lquery) {
221 /* Writing concatenation */
222 i__3[0] = 1, a__1[0] = side;
223 i__3[1] = 1, a__1[1] = trans;
224 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
227 nb = ilaenv_(&c__1, "SORMQR", ch__1, &i__1, n, &i__2, &c_n1);
229 /* Writing concatenation */
230 i__3[0] = 1, a__1[0] = side;
231 i__3[1] = 1, a__1[1] = trans;
232 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
235 nb = ilaenv_(&c__1, "SORMQR", ch__1, m, &i__1, &i__2, &c_n1);
239 /* Writing concatenation */
240 i__3[0] = 1, a__1[0] = side;
241 i__3[1] = 1, a__1[1] = trans;
242 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
245 nb = ilaenv_(&c__1, "SORMLQ", ch__1, &i__1, n, &i__2, &c_n1);
247 /* Writing concatenation */
248 i__3[0] = 1, a__1[0] = side;
249 i__3[1] = 1, a__1[1] = trans;
250 s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
253 nb = ilaenv_(&c__1, "SORMLQ", ch__1, m, &i__1, &i__2, &c_n1);
256 lwkopt = max(1,nw) * nb;
257 work[1] = (real) lwkopt;
262 xerbla_("SORMBR", &i__1);
268 /* Quick return if possible */
271 if (*m == 0 || *n == 0) {
281 /* Q was determined by a call to SGEBRD with nq >= k */
283 sormqr_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
284 c_offset], ldc, &work[1], lwork, &iinfo);
287 /* Q was determined by a call to SGEBRD with nq < k */
301 sormqr_(side, trans, &mi, &ni, &i__1, &a[a_dim1 + 2], lda, &tau[1]
302 , &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
309 *(unsigned char *)transt = 'T';
311 *(unsigned char *)transt = 'N';
315 /* P was determined by a call to SGEBRD with nq > k */
317 sormlq_(side, transt, m, n, k, &a[a_offset], lda, &tau[1], &c__[
318 c_offset], ldc, &work[1], lwork, &iinfo);
321 /* P was determined by a call to SGEBRD with nq <= k */
335 sormlq_(side, transt, &mi, &ni, &i__1, &a[(a_dim1 << 1) + 1], lda,
336 &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &
340 work[1] = (real) lwkopt;