3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
8 /* Subroutine */ int sorgbr_(char *vect, integer *m, integer *n, integer *k,
9 real *a, integer *lda, real *tau, real *work, integer *lwork, integer
12 /* System generated locals */
13 integer a_dim1, a_offset, i__1, i__2, i__3;
16 integer i__, j, nb, mn;
17 extern logical lsame_(char *, char *);
20 extern /* Subroutine */ int xerbla_(char *, integer *);
21 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
22 integer *, integer *);
23 extern /* Subroutine */ int sorglq_(integer *, integer *, integer *, real
24 *, integer *, real *, real *, integer *, integer *), sorgqr_(
25 integer *, integer *, integer *, real *, integer *, real *, real *
26 , integer *, integer *);
31 /* -- LAPACK routine (version 3.1) -- */
32 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
35 /* .. Scalar Arguments .. */
37 /* .. Array Arguments .. */
43 /* SORGBR generates one of the real orthogonal matrices Q or P**T */
44 /* determined by SGEBRD when reducing a real matrix A to bidiagonal */
45 /* form: A = Q * B * P**T. Q and P**T are defined as products of */
46 /* elementary reflectors H(i) or G(i) respectively. */
48 /* If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q */
50 /* if m >= k, Q = H(1) H(2) . . . H(k) and SORGBR returns the first n */
51 /* columns of Q, where m >= n >= k; */
52 /* if m < k, Q = H(1) H(2) . . . H(m-1) and SORGBR returns Q as an */
55 /* If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T */
57 /* if k < n, P**T = G(k) . . . G(2) G(1) and SORGBR returns the first m */
58 /* rows of P**T, where n >= m >= k; */
59 /* if k >= n, P**T = G(n-1) . . . G(2) G(1) and SORGBR returns P**T as */
60 /* an N-by-N matrix. */
65 /* VECT (input) CHARACTER*1 */
66 /* Specifies whether the matrix Q or the matrix P**T is */
67 /* required, as defined in the transformation applied by SGEBRD: */
68 /* = 'Q': generate Q; */
69 /* = 'P': generate P**T. */
71 /* M (input) INTEGER */
72 /* The number of rows of the matrix Q or P**T to be returned. */
75 /* N (input) INTEGER */
76 /* The number of columns of the matrix Q or P**T to be returned. */
78 /* If VECT = 'Q', M >= N >= min(M,K); */
79 /* if VECT = 'P', N >= M >= min(N,K). */
81 /* K (input) INTEGER */
82 /* If VECT = 'Q', the number of columns in the original M-by-K */
83 /* matrix reduced by SGEBRD. */
84 /* If VECT = 'P', the number of rows in the original K-by-N */
85 /* matrix reduced by SGEBRD. */
88 /* A (input/output) REAL array, dimension (LDA,N) */
89 /* On entry, the vectors which define the elementary reflectors, */
90 /* as returned by SGEBRD. */
91 /* On exit, the M-by-N matrix Q or P**T. */
93 /* LDA (input) INTEGER */
94 /* The leading dimension of the array A. LDA >= max(1,M). */
96 /* TAU (input) REAL array, dimension */
97 /* (min(M,K)) if VECT = 'Q' */
98 /* (min(N,K)) if VECT = 'P' */
99 /* TAU(i) must contain the scalar factor of the elementary */
100 /* reflector H(i) or G(i), which determines Q or P**T, as */
101 /* returned by SGEBRD in its array argument TAUQ or TAUP. */
103 /* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
104 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
106 /* LWORK (input) INTEGER */
107 /* The dimension of the array WORK. LWORK >= max(1,min(M,N)). */
108 /* For optimum performance LWORK >= min(M,N)*NB, where NB */
109 /* is the optimal blocksize. */
111 /* If LWORK = -1, then a workspace query is assumed; the routine */
112 /* only calculates the optimal size of the WORK array, returns */
113 /* this value as the first entry of the WORK array, and no error */
114 /* message related to LWORK is issued by XERBLA. */
116 /* INFO (output) INTEGER */
117 /* = 0: successful exit */
118 /* < 0: if INFO = -i, the i-th argument had an illegal value */
120 /* ===================================================================== */
122 /* .. Parameters .. */
124 /* .. Local Scalars .. */
126 /* .. External Functions .. */
128 /* .. External Subroutines .. */
130 /* .. Intrinsic Functions .. */
132 /* .. Executable Statements .. */
134 /* Test the input arguments */
136 /* Parameter adjustments */
138 a_offset = 1 + a_dim1;
145 wantq = lsame_(vect, "Q");
147 lquery = *lwork == -1;
148 if (! wantq && ! lsame_(vect, "P")) {
152 } else if (*n < 0 || wantq && (*n > *m || *n < min(*m,*k)) || ! wantq && (
153 *m > *n || *m < min(*n,*k))) {
157 } else if (*lda < max(1,*m)) {
159 } else if (*lwork < max(1,mn) && ! lquery) {
165 nb = ilaenv_(&c__1, "SORGQR", " ", m, n, k, &c_n1);
167 nb = ilaenv_(&c__1, "SORGLQ", " ", m, n, k, &c_n1);
169 lwkopt = max(1,mn) * nb;
170 work[1] = (real) lwkopt;
175 xerbla_("SORGBR", &i__1);
181 /* Quick return if possible */
183 if (*m == 0 || *n == 0) {
190 /* Form Q, determined by a call to SGEBRD to reduce an m-by-k */
195 /* If m >= k, assume m >= n >= k */
197 sorgqr_(m, n, k, &a[a_offset], lda, &tau[1], &work[1], lwork, &
202 /* If m < k, assume m = n */
204 /* Shift the vectors which define the elementary reflectors one */
205 /* column to the right, and set the first row and column of Q */
206 /* to those of the unit matrix */
208 for (j = *m; j >= 2; --j) {
209 a[j * a_dim1 + 1] = 0.f;
211 for (i__ = j + 1; i__ <= i__1; ++i__) {
212 a[i__ + j * a_dim1] = a[i__ + (j - 1) * a_dim1];
219 for (i__ = 2; i__ <= i__1; ++i__) {
220 a[i__ + a_dim1] = 0.f;
225 /* Form Q(2:m,2:m) */
230 sorgqr_(&i__1, &i__2, &i__3, &a[(a_dim1 << 1) + 2], lda, &tau[
231 1], &work[1], lwork, &iinfo);
236 /* Form P', determined by a call to SGEBRD to reduce a k-by-n */
241 /* If k < n, assume k <= m <= n */
243 sorglq_(m, n, k, &a[a_offset], lda, &tau[1], &work[1], lwork, &
248 /* If k >= n, assume m = n */
250 /* Shift the vectors which define the elementary reflectors one */
251 /* row downward, and set the first row and column of P' to */
252 /* those of the unit matrix */
256 for (i__ = 2; i__ <= i__1; ++i__) {
257 a[i__ + a_dim1] = 0.f;
261 for (j = 2; j <= i__1; ++j) {
262 for (i__ = j - 1; i__ >= 2; --i__) {
263 a[i__ + j * a_dim1] = a[i__ - 1 + j * a_dim1];
266 a[j * a_dim1 + 1] = 0.f;
271 /* Form P'(2:n,2:n) */
276 sorglq_(&i__1, &i__2, &i__3, &a[(a_dim1 << 1) + 2], lda, &tau[
277 1], &work[1], lwork, &iinfo);
281 work[1] = (real) lwkopt;