3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c_n1 = -1;
7 static doublereal c_b33 = 0.;
8 static integer c__0 = 0;
10 /* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer *
11 nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
12 doublereal *work, 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;
19 doublereal anrm, bnrm;
23 extern logical lsame_(char *, char *);
26 extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
27 extern doublereal dlamch_(char *), dlange_(char *, integer *,
28 integer *, doublereal *, integer *, doublereal *);
29 extern /* Subroutine */ int dgelqf_(integer *, integer *, doublereal *,
30 integer *, doublereal *, doublereal *, integer *, integer *),
31 dlascl_(char *, integer *, integer *, doublereal *, doublereal *,
32 integer *, integer *, doublereal *, integer *, integer *),
33 dgeqrf_(integer *, integer *, doublereal *, integer *,
34 doublereal *, doublereal *, integer *, integer *), dlaset_(char *,
35 integer *, integer *, doublereal *, doublereal *, doublereal *,
36 integer *), xerbla_(char *, integer *);
37 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
38 integer *, integer *);
41 extern /* Subroutine */ int dormlq_(char *, char *, integer *, integer *,
42 integer *, doublereal *, integer *, doublereal *, doublereal *,
43 integer *, doublereal *, integer *, integer *),
44 dormqr_(char *, char *, integer *, integer *, integer *,
45 doublereal *, integer *, doublereal *, doublereal *, integer *,
46 doublereal *, integer *, integer *);
49 extern /* Subroutine */ int dtrtrs_(char *, char *, char *, integer *,
50 integer *, doublereal *, integer *, doublereal *, integer *,
54 /* -- LAPACK driver routine (version 3.1) -- */
55 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
58 /* .. Scalar Arguments .. */
60 /* .. Array Arguments .. */
66 /* DGELS solves overdetermined or underdetermined real linear systems */
67 /* involving an M-by-N matrix A, or its transpose, using a QR or LQ */
68 /* factorization of A. It is assumed that A has full rank. */
70 /* The following options are provided: */
72 /* 1. If TRANS = 'N' and m >= n: find the least squares solution of */
73 /* an overdetermined system, i.e., solve the least squares problem */
74 /* minimize || B - A*X ||. */
76 /* 2. If TRANS = 'N' and m < n: find the minimum norm solution of */
77 /* an underdetermined system A * X = B. */
79 /* 3. If TRANS = 'T' and m >= n: find the minimum norm solution of */
80 /* an undetermined system A**T * X = B. */
82 /* 4. If TRANS = 'T' and m < n: find the least squares solution of */
83 /* an overdetermined system, i.e., solve the least squares problem */
84 /* minimize || B - A**T * X ||. */
86 /* Several right hand side vectors b and solution vectors x can be */
87 /* handled in a single call; they are stored as the columns of the */
88 /* M-by-NRHS right hand side matrix B and the N-by-NRHS solution */
94 /* TRANS (input) CHARACTER*1 */
95 /* = 'N': the linear system involves A; */
96 /* = 'T': the linear system involves A**T. */
98 /* M (input) INTEGER */
99 /* The number of rows of the matrix A. M >= 0. */
101 /* N (input) INTEGER */
102 /* The number of columns of the matrix A. N >= 0. */
104 /* NRHS (input) INTEGER */
105 /* The number of right hand sides, i.e., the number of */
106 /* columns of the matrices B and X. NRHS >=0. */
108 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
109 /* On entry, the M-by-N matrix A. */
111 /* if M >= N, A is overwritten by details of its QR */
112 /* factorization as returned by DGEQRF; */
113 /* if M < N, A is overwritten by details of its LQ */
114 /* factorization as returned by DGELQF. */
116 /* LDA (input) INTEGER */
117 /* The leading dimension of the array A. LDA >= max(1,M). */
119 /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
120 /* On entry, the matrix B of right hand side vectors, stored */
121 /* columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS */
122 /* if TRANS = 'T'. */
123 /* On exit, if INFO = 0, B is overwritten by the solution */
124 /* vectors, stored columnwise: */
125 /* if TRANS = 'N' and m >= n, rows 1 to n of B contain the least */
126 /* squares solution vectors; the residual sum of squares for the */
127 /* solution in each column is given by the sum of squares of */
128 /* elements N+1 to M in that column; */
129 /* if TRANS = 'N' and m < n, rows 1 to N 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 /* minimum norm solution vectors; */
133 /* if TRANS = 'T' and m < n, rows 1 to M of B contain the */
134 /* least squares solution vectors; the residual sum of squares */
135 /* for the solution in each column is given by the sum of */
136 /* squares of elements M+1 to N in that column. */
138 /* LDB (input) INTEGER */
139 /* The leading dimension of the array B. LDB >= MAX(1,M,N). */
141 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
142 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
144 /* LWORK (input) INTEGER */
145 /* The dimension of the array WORK. */
146 /* LWORK >= max( 1, MN + max( MN, NRHS ) ). */
147 /* For optimal performance, */
148 /* LWORK >= max( 1, MN + max( MN, NRHS )*NB ). */
149 /* where MN = min(M,N) and NB is the optimum block size. */
151 /* If LWORK = -1, then a workspace query is assumed; the routine */
152 /* only calculates the optimal size of the WORK array, returns */
153 /* this value as the first entry of the WORK array, and no error */
154 /* message related to LWORK is issued by XERBLA. */
156 /* INFO (output) INTEGER */
157 /* = 0: successful exit */
158 /* < 0: if INFO = -i, the i-th argument had an illegal value */
159 /* > 0: if INFO = i, the i-th diagonal element of the */
160 /* triangular factor of A is zero, so that A does not have */
161 /* full rank; the least squares solution could not be */
164 /* ===================================================================== */
166 /* .. Parameters .. */
168 /* .. Local Scalars .. */
170 /* .. Local Arrays .. */
172 /* .. External Functions .. */
174 /* .. External Subroutines .. */
176 /* .. Intrinsic Functions .. */
178 /* .. Executable Statements .. */
180 /* Test the input arguments. */
182 /* Parameter adjustments */
184 a_offset = 1 + a_dim1;
187 b_offset = 1 + b_dim1;
194 lquery = *lwork == -1;
195 if (! (lsame_(trans, "N") || lsame_(trans, "T"))) {
201 } else if (*nrhs < 0) {
203 } else if (*lda < max(1,*m)) {
205 } else /* if(complicated condition) */ {
208 if (*ldb < max(i__1,*n)) {
210 } else /* if(complicated condition) */ {
212 i__1 = 1, i__2 = mn + max(mn,*nrhs);
213 if (*lwork < max(i__1,i__2) && ! lquery) {
219 /* Figure out optimal block size */
221 if (*info == 0 || *info == -10) {
224 if (lsame_(trans, "N")) {
229 nb = ilaenv_(&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1);
232 i__1 = nb, i__2 = ilaenv_(&c__1, "DORMQR", "LN", m, nrhs, n, &
237 i__1 = nb, i__2 = ilaenv_(&c__1, "DORMQR", "LT", m, nrhs, n, &
242 nb = ilaenv_(&c__1, "DGELQF", " ", m, n, &c_n1, &c_n1);
245 i__1 = nb, i__2 = ilaenv_(&c__1, "DORMLQ", "LT", n, nrhs, m, &
250 i__1 = nb, i__2 = ilaenv_(&c__1, "DORMLQ", "LN", n, nrhs, m, &
257 i__1 = 1, i__2 = mn + max(mn,*nrhs) * nb;
258 wsize = max(i__1,i__2);
259 work[1] = (doublereal) wsize;
265 xerbla_("DGELS ", &i__1);
271 /* Quick return if possible */
275 if (min(i__1,*nrhs) == 0) {
277 dlaset_("Full", &i__1, nrhs, &c_b33, &c_b33, &b[b_offset], ldb);
281 /* Get machine parameters */
283 smlnum = dlamch_("S") / dlamch_("P");
284 bignum = 1. / smlnum;
285 dlabad_(&smlnum, &bignum);
287 /* Scale A, B if max element outside range [SMLNUM,BIGNUM] */
289 anrm = dlange_("M", m, n, &a[a_offset], lda, rwork);
291 if (anrm > 0. && anrm < smlnum) {
293 /* Scale matrix norm up to SMLNUM */
295 dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda,
298 } else if (anrm > bignum) {
300 /* Scale matrix norm down to BIGNUM */
302 dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda,
305 } else if (anrm == 0.) {
307 /* Matrix all zero. Return zero solution. */
310 dlaset_("F", &i__1, nrhs, &c_b33, &c_b33, &b[b_offset], ldb);
318 bnrm = dlange_("M", &brow, nrhs, &b[b_offset], ldb, rwork);
320 if (bnrm > 0. && bnrm < smlnum) {
322 /* Scale matrix norm up to SMLNUM */
324 dlascl_("G", &c__0, &c__0, &bnrm, &smlnum, &brow, nrhs, &b[b_offset],
327 } else if (bnrm > bignum) {
329 /* Scale matrix norm down to BIGNUM */
331 dlascl_("G", &c__0, &c__0, &bnrm, &bignum, &brow, nrhs, &b[b_offset],
338 /* compute QR factorization of A */
341 dgeqrf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
344 /* workspace at least N, optimally N*NB */
348 /* Least-Squares Problem min || A * X - B || */
350 /* B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
353 dormqr_("Left", "Transpose", m, nrhs, n, &a[a_offset], lda, &work[
354 1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
356 /* workspace at least NRHS, optimally NRHS*NB */
358 /* B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS) */
360 dtrtrs_("Upper", "No transpose", "Non-unit", n, nrhs, &a[a_offset]
361 , lda, &b[b_offset], ldb, info);
371 /* Overdetermined system of equations A' * X = B */
373 /* B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS) */
375 dtrtrs_("Upper", "Transpose", "Non-unit", n, nrhs, &a[a_offset],
376 lda, &b[b_offset], ldb, info);
382 /* B(N+1:M,1:NRHS) = ZERO */
385 for (j = 1; j <= i__1; ++j) {
387 for (i__ = *n + 1; i__ <= i__2; ++i__) {
388 b[i__ + j * b_dim1] = 0.;
394 /* B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS) */
397 dormqr_("Left", "No transpose", m, nrhs, n, &a[a_offset], lda, &
398 work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
400 /* workspace at least NRHS, optimally NRHS*NB */
408 /* Compute LQ factorization of A */
411 dgelqf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
414 /* workspace at least M, optimally M*NB. */
418 /* underdetermined system of equations A * X = B */
420 /* B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS) */
422 dtrtrs_("Lower", "No transpose", "Non-unit", m, nrhs, &a[a_offset]
423 , lda, &b[b_offset], ldb, info);
429 /* B(M+1:N,1:NRHS) = 0 */
432 for (j = 1; j <= i__1; ++j) {
434 for (i__ = *m + 1; i__ <= i__2; ++i__) {
435 b[i__ + j * b_dim1] = 0.;
441 /* B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS) */
444 dormlq_("Left", "Transpose", n, nrhs, m, &a[a_offset], lda, &work[
445 1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
447 /* workspace at least NRHS, optimally NRHS*NB */
453 /* overdetermined system min || A' * X - B || */
455 /* B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) */
458 dormlq_("Left", "No transpose", n, nrhs, m, &a[a_offset], lda, &
459 work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
461 /* workspace at least NRHS, optimally NRHS*NB */
463 /* B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS) */
465 dtrtrs_("Lower", "Transpose", "Non-unit", m, nrhs, &a[a_offset],
466 lda, &b[b_offset], ldb, info);
481 dlascl_("G", &c__0, &c__0, &anrm, &smlnum, &scllen, nrhs, &b[b_offset]
483 } else if (iascl == 2) {
484 dlascl_("G", &c__0, &c__0, &anrm, &bignum, &scllen, nrhs, &b[b_offset]
488 dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, &scllen, nrhs, &b[b_offset]
490 } else if (ibscl == 2) {
491 dlascl_("G", &c__0, &c__0, &bignum, &bnrm, &scllen, nrhs, &b[b_offset]
496 work[1] = (doublereal) wsize;