3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static real c_b33 = 0.f;
8 static integer c__0 = 0;
10 /* Subroutine */ int sgels_(char *trans, integer *m, integer *n, integer *
11 nrhs, real *a, integer *lda, real *b, integer *ldb, real *work,
12 integer *lwork, integer *info)
14 /* System generated locals */
15 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
18 integer i__, j, nb, mn;
23 extern logical lsame_(char *, char *);
26 extern /* Subroutine */ int slabad_(real *, real *);
27 extern doublereal slamch_(char *), slange_(char *, integer *,
28 integer *, real *, integer *, real *);
29 extern /* Subroutine */ int xerbla_(char *, integer *);
30 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
31 integer *, integer *);
34 extern /* Subroutine */ int sgelqf_(integer *, integer *, real *, integer
35 *, real *, real *, integer *, integer *), slascl_(char *, integer
36 *, integer *, real *, real *, integer *, integer *, real *,
37 integer *, integer *), sgeqrf_(integer *, integer *, real
38 *, integer *, real *, real *, integer *, integer *), slaset_(char
39 *, integer *, integer *, real *, real *, real *, integer *);
41 extern /* Subroutine */ int sormlq_(char *, char *, integer *, integer *,
42 integer *, real *, integer *, real *, real *, integer *, real *,
43 integer *, integer *);
45 extern /* Subroutine */ int sormqr_(char *, char *, integer *, integer *,
46 integer *, real *, integer *, real *, real *, integer *, real *,
47 integer *, integer *), strtrs_(char *, char *,
48 char *, integer *, integer *, real *, integer *, real *, integer *
52 /* -- LAPACK driver routine (version 3.1) -- */
53 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
56 /* .. Scalar Arguments .. */
58 /* .. Array Arguments .. */
64 /* SGELS solves overdetermined or underdetermined real linear systems */
65 /* involving an M-by-N matrix A, or its transpose, using a QR or LQ */
66 /* factorization of A. It is assumed that A has full rank. */
68 /* The following options are provided: */
70 /* 1. If TRANS = 'N' and m >= n: find the least squares solution of */
71 /* an overdetermined system, i.e., solve the least squares problem */
72 /* minimize || B - A*X ||. */
74 /* 2. If TRANS = 'N' and m < n: find the minimum norm solution of */
75 /* an underdetermined system A * X = B. */
77 /* 3. If TRANS = 'T' and m >= n: find the minimum norm solution of */
78 /* an undetermined system A**T * X = B. */
80 /* 4. If TRANS = 'T' and m < n: find the least squares solution of */
81 /* an overdetermined system, i.e., solve the least squares problem */
82 /* minimize || B - A**T * X ||. */
84 /* Several right hand side vectors b and solution vectors x can be */
85 /* handled in a single call; they are stored as the columns of the */
86 /* M-by-NRHS right hand side matrix B and the N-by-NRHS solution */
92 /* TRANS (input) CHARACTER*1 */
93 /* = 'N': the linear system involves A; */
94 /* = 'T': the linear system involves A**T. */
96 /* M (input) INTEGER */
97 /* The number of rows of the matrix A. M >= 0. */
99 /* N (input) INTEGER */
100 /* The number of columns of the matrix A. N >= 0. */
102 /* NRHS (input) INTEGER */
103 /* The number of right hand sides, i.e., the number of */
104 /* columns of the matrices B and X. NRHS >=0. */
106 /* A (input/output) REAL array, dimension (LDA,N) */
107 /* On entry, the M-by-N matrix A. */
109 /* if M >= N, A is overwritten by details of its QR */
110 /* factorization as returned by SGEQRF; */
111 /* if M < N, A is overwritten by details of its LQ */
112 /* factorization as returned by SGELQF. */
114 /* LDA (input) INTEGER */
115 /* The leading dimension of the array A. LDA >= max(1,M). */
117 /* B (input/output) REAL array, dimension (LDB,NRHS) */
118 /* On entry, the matrix B of right hand side vectors, stored */
119 /* columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS */
120 /* if TRANS = 'T'. */
121 /* On exit, if INFO = 0, B is overwritten by the solution */
122 /* vectors, stored columnwise: */
123 /* if TRANS = 'N' and m >= n, rows 1 to n of B contain the least */
124 /* squares solution vectors; the residual sum of squares for the */
125 /* solution in each column is given by the sum of squares of */
126 /* elements N+1 to M in that column; */
127 /* if TRANS = 'N' and m < n, rows 1 to N of B contain the */
128 /* minimum norm solution vectors; */
129 /* if TRANS = 'T' and m >= n, rows 1 to M of B contain the */
130 /* minimum norm solution vectors; */
131 /* if TRANS = 'T' and m < n, rows 1 to M of B contain the */
132 /* least squares solution vectors; the residual sum of squares */
133 /* for the solution in each column is given by the sum of */
134 /* squares of elements M+1 to N in that column. */
136 /* LDB (input) INTEGER */
137 /* The leading dimension of the array B. LDB >= MAX(1,M,N). */
139 /* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
140 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
142 /* LWORK (input) INTEGER */
143 /* The dimension of the array WORK. */
144 /* LWORK >= max( 1, MN + max( MN, NRHS ) ). */
145 /* For optimal performance, */
146 /* LWORK >= max( 1, MN + max( MN, NRHS )*NB ). */
147 /* where MN = min(M,N) and NB is the optimum block size. */
149 /* If LWORK = -1, then a workspace query is assumed; the routine */
150 /* only calculates the optimal size of the WORK array, returns */
151 /* this value as the first entry of the WORK array, and no error */
152 /* message related to LWORK is issued by XERBLA. */
154 /* INFO (output) INTEGER */
155 /* = 0: successful exit */
156 /* < 0: if INFO = -i, the i-th argument had an illegal value */
157 /* > 0: if INFO = i, the i-th diagonal element of the */
158 /* triangular factor of A is zero, so that A does not have */
159 /* full rank; the least squares solution could not be */
162 /* ===================================================================== */
164 /* .. Parameters .. */
166 /* .. Local Scalars .. */
168 /* .. Local Arrays .. */
170 /* .. External Functions .. */
172 /* .. External Subroutines .. */
174 /* .. Intrinsic Functions .. */
176 /* .. Executable Statements .. */
178 /* Test the input arguments. */
180 /* Parameter adjustments */
182 a_offset = 1 + a_dim1;
185 b_offset = 1 + b_dim1;
192 lquery = *lwork == -1;
193 if (! (lsame_(trans, "N") || lsame_(trans, "T"))) {
199 } else if (*nrhs < 0) {
201 } else if (*lda < max(1,*m)) {
203 } else /* if(complicated condition) */ {
206 if (*ldb < max(i__1,*n)) {
208 } else /* if(complicated condition) */ {
210 i__1 = 1, i__2 = mn + max(mn,*nrhs);
211 if (*lwork < max(i__1,i__2) && ! lquery) {
217 /* Figure out optimal block size */
219 if (*info == 0 || *info == -10) {
222 if (lsame_(trans, "N")) {
227 nb = ilaenv_(&c__1, "SGEQRF", " ", m, n, &c_n1, &c_n1);
230 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMQR", "LN", m, nrhs, n, &
235 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMQR", "LT", m, nrhs, n, &
240 nb = ilaenv_(&c__1, "SGELQF", " ", m, n, &c_n1, &c_n1);
243 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMLQ", "LT", n, nrhs, m, &
248 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMLQ", "LN", n, nrhs, m, &
255 i__1 = 1, i__2 = mn + max(mn,*nrhs) * nb;
256 wsize = max(i__1,i__2);
257 work[1] = (real) wsize;
263 xerbla_("SGELS ", &i__1);
269 /* Quick return if possible */
273 if (min(i__1,*nrhs) == 0) {
275 slaset_("Full", &i__1, nrhs, &c_b33, &c_b33, &b[b_offset], ldb);
279 /* Get machine parameters */
281 smlnum = slamch_("S") / slamch_("P");
282 bignum = 1.f / smlnum;
283 slabad_(&smlnum, &bignum);
285 /* Scale A, B if max element outside range [SMLNUM,BIGNUM] */
287 anrm = slange_("M", m, n, &a[a_offset], lda, rwork);
289 if (anrm > 0.f && anrm < smlnum) {
291 /* Scale matrix norm up to SMLNUM */
293 slascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda,
296 } else if (anrm > bignum) {
298 /* Scale matrix norm down to BIGNUM */
300 slascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda,
303 } else if (anrm == 0.f) {
305 /* Matrix all zero. Return zero solution. */
308 slaset_("F", &i__1, nrhs, &c_b33, &c_b33, &b[b_offset], ldb);
316 bnrm = slange_("M", &brow, nrhs, &b[b_offset], ldb, rwork);
318 if (bnrm > 0.f && bnrm < smlnum) {
320 /* Scale matrix norm up to SMLNUM */
322 slascl_("G", &c__0, &c__0, &bnrm, &smlnum, &brow, nrhs, &b[b_offset],
325 } else if (bnrm > bignum) {
327 /* Scale matrix norm down to BIGNUM */
329 slascl_("G", &c__0, &c__0, &bnrm, &bignum, &brow, nrhs, &b[b_offset],
336 /* compute QR factorization of A */
339 sgeqrf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
342 /* workspace at least N, optimally N*NB */
346 /* Least-Squares Problem min || A * X - B || */
348 /* B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
351 sormqr_("Left", "Transpose", m, nrhs, n, &a[a_offset], lda, &work[
352 1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
354 /* workspace at least NRHS, optimally NRHS*NB */
356 /* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS) */
358 strtrs_("Upper", "No transpose", "Non-unit", n, nrhs, &a[a_offset]
359 , lda, &b[b_offset], ldb, info);
369 /* Overdetermined system of equations A' * X = B */
371 /* B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS) */
373 strtrs_("Upper", "Transpose", "Non-unit", n, nrhs, &a[a_offset],
374 lda, &b[b_offset], ldb, info);
380 /* B(N+1:M,1:NRHS) = ZERO */
383 for (j = 1; j <= i__1; ++j) {
385 for (i__ = *n + 1; i__ <= i__2; ++i__) {
386 b[i__ + j * b_dim1] = 0.f;
392 /* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS) */
395 sormqr_("Left", "No transpose", m, nrhs, n, &a[a_offset], lda, &
396 work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
398 /* workspace at least NRHS, optimally NRHS*NB */
406 /* Compute LQ factorization of A */
409 sgelqf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
412 /* workspace at least M, optimally M*NB. */
416 /* underdetermined system of equations A * X = B */
418 /* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS) */
420 strtrs_("Lower", "No transpose", "Non-unit", m, nrhs, &a[a_offset]
421 , lda, &b[b_offset], ldb, info);
427 /* B(M+1:N,1:NRHS) = 0 */
430 for (j = 1; j <= i__1; ++j) {
432 for (i__ = *m + 1; i__ <= i__2; ++i__) {
433 b[i__ + j * b_dim1] = 0.f;
439 /* B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS) */
442 sormlq_("Left", "Transpose", n, nrhs, m, &a[a_offset], lda, &work[
443 1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
445 /* workspace at least NRHS, optimally NRHS*NB */
451 /* overdetermined system min || A' * X - B || */
453 /* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) */
456 sormlq_("Left", "No transpose", n, nrhs, m, &a[a_offset], lda, &
457 work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
459 /* workspace at least NRHS, optimally NRHS*NB */
461 /* B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS) */
463 strtrs_("Lower", "Transpose", "Non-unit", m, nrhs, &a[a_offset],
464 lda, &b[b_offset], ldb, info);
479 slascl_("G", &c__0, &c__0, &anrm, &smlnum, &scllen, nrhs, &b[b_offset]
481 } else if (iascl == 2) {
482 slascl_("G", &c__0, &c__0, &anrm, &bignum, &scllen, nrhs, &b[b_offset]
486 slascl_("G", &c__0, &c__0, &smlnum, &bnrm, &scllen, nrhs, &b[b_offset]
488 } else if (ibscl == 2) {
489 slascl_("G", &c__0, &c__0, &bignum, &bnrm, &scllen, nrhs, &b[b_offset]
494 work[1] = (real) wsize;