1 SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
2 $ INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER TRANS
11 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGELS solves overdetermined or underdetermined real linear systems
21 * involving an M-by-N matrix A, or its transpose, using a QR or LQ
22 * factorization of A. It is assumed that A has full rank.
23 *
24 * The following options are provided:
25 *
26 * 1. If TRANS = 'N' and m >= n: find the least squares solution of
27 * an overdetermined system, i.e., solve the least squares problem
28 * minimize || B - A*X ||.
29 *
30 * 2. If TRANS = 'N' and m < n: find the minimum norm solution of
31 * an underdetermined system A * X = B.
32 *
33 * 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
34 * an undetermined system A**T * X = B.
35 *
36 * 4. If TRANS = 'T' and m < n: find the least squares solution of
37 * an overdetermined system, i.e., solve the least squares problem
38 * minimize || B - A**T * X ||.
39 *
40 * Several right hand side vectors b and solution vectors x can be
41 * handled in a single call; they are stored as the columns of the
42 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
43 * matrix X.
44 *
45 * Arguments
46 * =========
47 *
48 * TRANS (input) CHARACTER*1
49 * = 'N': the linear system involves A;
50 * = 'T': the linear system involves A**T.
51 *
52 * M (input) INTEGER
53 * The number of rows of the matrix A. M >= 0.
54 *
55 * N (input) INTEGER
56 * The number of columns of the matrix A. N >= 0.
57 *
58 * NRHS (input) INTEGER
59 * The number of right hand sides, i.e., the number of
60 * columns of the matrices B and X. NRHS >=0.
61 *
62 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
63 * On entry, the M-by-N matrix A.
64 * On exit,
65 * if M >= N, A is overwritten by details of its QR
66 * factorization as returned by DGEQRF;
67 * if M < N, A is overwritten by details of its LQ
68 * factorization as returned by DGELQF.
69 *
70 * LDA (input) INTEGER
71 * The leading dimension of the array A. LDA >= max(1,M).
72 *
73 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
74 * On entry, the matrix B of right hand side vectors, stored
75 * columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
76 * if TRANS = 'T'.
77 * On exit, if INFO = 0, B is overwritten by the solution
78 * vectors, stored columnwise:
79 * if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
80 * squares solution vectors; the residual sum of squares for the
81 * solution in each column is given by the sum of squares of
82 * elements N+1 to M in that column;
83 * if TRANS = 'N' and m < n, rows 1 to N of B contain the
84 * minimum norm solution vectors;
85 * if TRANS = 'T' and m >= n, rows 1 to M of B contain the
86 * minimum norm solution vectors;
87 * if TRANS = 'T' and m < n, rows 1 to M of B contain the
88 * least squares solution vectors; the residual sum of squares
89 * for the solution in each column is given by the sum of
90 * squares of elements M+1 to N in that column.
91 *
92 * LDB (input) INTEGER
93 * The leading dimension of the array B. LDB >= MAX(1,M,N).
94 *
95 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
96 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
97 *
98 * LWORK (input) INTEGER
99 * The dimension of the array WORK.
100 * LWORK >= max( 1, MN + max( MN, NRHS ) ).
101 * For optimal performance,
102 * LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
103 * where MN = min(M,N) and NB is the optimum block size.
104 *
105 * If LWORK = -1, then a workspace query is assumed; the routine
106 * only calculates the optimal size of the WORK array, returns
107 * this value as the first entry of the WORK array, and no error
108 * message related to LWORK is issued by XERBLA.
109 *
110 * INFO (output) INTEGER
111 * = 0: successful exit
112 * < 0: if INFO = -i, the i-th argument had an illegal value
113 * > 0: if INFO = i, the i-th diagonal element of the
114 * triangular factor of A is zero, so that A does not have
115 * full rank; the least squares solution could not be
116 * computed.
117 *
118 * =====================================================================
119 *
120 * .. Parameters ..
121 DOUBLE PRECISION ZERO, ONE
122 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
123 * ..
124 * .. Local Scalars ..
125 LOGICAL LQUERY, TPSD
126 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
127 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
128 * ..
129 * .. Local Arrays ..
130 DOUBLE PRECISION RWORK( 1 )
131 * ..
132 * .. External Functions ..
133 LOGICAL LSAME
134 INTEGER ILAENV
135 DOUBLE PRECISION DLAMCH, DLANGE
136 EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
137 * ..
138 * .. External Subroutines ..
139 EXTERNAL DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR,
140 $ DTRTRS, XERBLA
141 * ..
142 * .. Intrinsic Functions ..
143 INTRINSIC DBLE, MAX, MIN
144 * ..
145 * .. Executable Statements ..
146 *
147 * Test the input arguments.
148 *
149 INFO = 0
150 MN = MIN( M, N )
151 LQUERY = ( LWORK.EQ.-1 )
152 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN
153 INFO = -1
154 ELSE IF( M.LT.0 ) THEN
155 INFO = -2
156 ELSE IF( N.LT.0 ) THEN
157 INFO = -3
158 ELSE IF( NRHS.LT.0 ) THEN
159 INFO = -4
160 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
161 INFO = -6
162 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
163 INFO = -8
164 ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
165 $ THEN
166 INFO = -10
167 END IF
168 *
169 * Figure out optimal block size
170 *
171 IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
172 *
173 TPSD = .TRUE.
174 IF( LSAME( TRANS, 'N' ) )
175 $ TPSD = .FALSE.
176 *
177 IF( M.GE.N ) THEN
178 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
179 IF( TPSD ) THEN
180 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N,
181 $ -1 ) )
182 ELSE
183 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N,
184 $ -1 ) )
185 END IF
186 ELSE
187 NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
188 IF( TPSD ) THEN
189 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M,
190 $ -1 ) )
191 ELSE
192 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M,
193 $ -1 ) )
194 END IF
195 END IF
196 *
197 WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
198 WORK( 1 ) = DBLE( WSIZE )
199 *
200 END IF
201 *
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'DGELS ', -INFO )
204 RETURN
205 ELSE IF( LQUERY ) THEN
206 RETURN
207 END IF
208 *
209 * Quick return if possible
210 *
211 IF( MIN( M, N, NRHS ).EQ.0 ) THEN
212 CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
213 RETURN
214 END IF
215 *
216 * Get machine parameters
217 *
218 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
219 BIGNUM = ONE / SMLNUM
220 CALL DLABAD( SMLNUM, BIGNUM )
221 *
222 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
223 *
224 ANRM = DLANGE( 'M', M, N, A, LDA, RWORK )
225 IASCL = 0
226 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
227 *
228 * Scale matrix norm up to SMLNUM
229 *
230 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
231 IASCL = 1
232 ELSE IF( ANRM.GT.BIGNUM ) THEN
233 *
234 * Scale matrix norm down to BIGNUM
235 *
236 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
237 IASCL = 2
238 ELSE IF( ANRM.EQ.ZERO ) THEN
239 *
240 * Matrix all zero. Return zero solution.
241 *
242 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
243 GO TO 50
244 END IF
245 *
246 BROW = M
247 IF( TPSD )
248 $ BROW = N
249 BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
250 IBSCL = 0
251 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
252 *
253 * Scale matrix norm up to SMLNUM
254 *
255 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
256 $ INFO )
257 IBSCL = 1
258 ELSE IF( BNRM.GT.BIGNUM ) THEN
259 *
260 * Scale matrix norm down to BIGNUM
261 *
262 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
263 $ INFO )
264 IBSCL = 2
265 END IF
266 *
267 IF( M.GE.N ) THEN
268 *
269 * compute QR factorization of A
270 *
271 CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
272 $ INFO )
273 *
274 * workspace at least N, optimally N*NB
275 *
276 IF( .NOT.TPSD ) THEN
277 *
278 * Least-Squares Problem min || A * X - B ||
279 *
280 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
281 *
282 CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA,
283 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
284 $ INFO )
285 *
286 * workspace at least NRHS, optimally NRHS*NB
287 *
288 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
289 *
290 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
291 $ A, LDA, B, LDB, INFO )
292 *
293 IF( INFO.GT.0 ) THEN
294 RETURN
295 END IF
296 *
297 SCLLEN = N
298 *
299 ELSE
300 *
301 * Overdetermined system of equations A**T * X = B
302 *
303 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
304 *
305 CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS,
306 $ A, LDA, B, LDB, INFO )
307 *
308 IF( INFO.GT.0 ) THEN
309 RETURN
310 END IF
311 *
312 * B(N+1:M,1:NRHS) = ZERO
313 *
314 DO 20 J = 1, NRHS
315 DO 10 I = N + 1, M
316 B( I, J ) = ZERO
317 10 CONTINUE
318 20 CONTINUE
319 *
320 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
321 *
322 CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
323 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
324 $ INFO )
325 *
326 * workspace at least NRHS, optimally NRHS*NB
327 *
328 SCLLEN = M
329 *
330 END IF
331 *
332 ELSE
333 *
334 * Compute LQ factorization of A
335 *
336 CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
337 $ INFO )
338 *
339 * workspace at least M, optimally M*NB.
340 *
341 IF( .NOT.TPSD ) THEN
342 *
343 * underdetermined system of equations A * X = B
344 *
345 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
346 *
347 CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
348 $ A, LDA, B, LDB, INFO )
349 *
350 IF( INFO.GT.0 ) THEN
351 RETURN
352 END IF
353 *
354 * B(M+1:N,1:NRHS) = 0
355 *
356 DO 40 J = 1, NRHS
357 DO 30 I = M + 1, N
358 B( I, J ) = ZERO
359 30 CONTINUE
360 40 CONTINUE
361 *
362 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
363 *
364 CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA,
365 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
366 $ INFO )
367 *
368 * workspace at least NRHS, optimally NRHS*NB
369 *
370 SCLLEN = N
371 *
372 ELSE
373 *
374 * overdetermined system min || A**T * X - B ||
375 *
376 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
377 *
378 CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
379 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
380 $ INFO )
381 *
382 * workspace at least NRHS, optimally NRHS*NB
383 *
384 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
385 *
386 CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
387 $ A, LDA, B, LDB, INFO )
388 *
389 IF( INFO.GT.0 ) THEN
390 RETURN
391 END IF
392 *
393 SCLLEN = M
394 *
395 END IF
396 *
397 END IF
398 *
399 * Undo scaling
400 *
401 IF( IASCL.EQ.1 ) THEN
402 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
403 $ INFO )
404 ELSE IF( IASCL.EQ.2 ) THEN
405 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
406 $ INFO )
407 END IF
408 IF( IBSCL.EQ.1 ) THEN
409 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
410 $ INFO )
411 ELSE IF( IBSCL.EQ.2 ) THEN
412 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
413 $ INFO )
414 END IF
415 *
416 50 CONTINUE
417 WORK( 1 ) = DBLE( WSIZE )
418 *
419 RETURN
420 *
421 * End of DGELS
422 *
423 END
2 $ INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER TRANS
11 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGELS solves overdetermined or underdetermined real linear systems
21 * involving an M-by-N matrix A, or its transpose, using a QR or LQ
22 * factorization of A. It is assumed that A has full rank.
23 *
24 * The following options are provided:
25 *
26 * 1. If TRANS = 'N' and m >= n: find the least squares solution of
27 * an overdetermined system, i.e., solve the least squares problem
28 * minimize || B - A*X ||.
29 *
30 * 2. If TRANS = 'N' and m < n: find the minimum norm solution of
31 * an underdetermined system A * X = B.
32 *
33 * 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
34 * an undetermined system A**T * X = B.
35 *
36 * 4. If TRANS = 'T' and m < n: find the least squares solution of
37 * an overdetermined system, i.e., solve the least squares problem
38 * minimize || B - A**T * X ||.
39 *
40 * Several right hand side vectors b and solution vectors x can be
41 * handled in a single call; they are stored as the columns of the
42 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
43 * matrix X.
44 *
45 * Arguments
46 * =========
47 *
48 * TRANS (input) CHARACTER*1
49 * = 'N': the linear system involves A;
50 * = 'T': the linear system involves A**T.
51 *
52 * M (input) INTEGER
53 * The number of rows of the matrix A. M >= 0.
54 *
55 * N (input) INTEGER
56 * The number of columns of the matrix A. N >= 0.
57 *
58 * NRHS (input) INTEGER
59 * The number of right hand sides, i.e., the number of
60 * columns of the matrices B and X. NRHS >=0.
61 *
62 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
63 * On entry, the M-by-N matrix A.
64 * On exit,
65 * if M >= N, A is overwritten by details of its QR
66 * factorization as returned by DGEQRF;
67 * if M < N, A is overwritten by details of its LQ
68 * factorization as returned by DGELQF.
69 *
70 * LDA (input) INTEGER
71 * The leading dimension of the array A. LDA >= max(1,M).
72 *
73 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
74 * On entry, the matrix B of right hand side vectors, stored
75 * columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
76 * if TRANS = 'T'.
77 * On exit, if INFO = 0, B is overwritten by the solution
78 * vectors, stored columnwise:
79 * if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
80 * squares solution vectors; the residual sum of squares for the
81 * solution in each column is given by the sum of squares of
82 * elements N+1 to M in that column;
83 * if TRANS = 'N' and m < n, rows 1 to N of B contain the
84 * minimum norm solution vectors;
85 * if TRANS = 'T' and m >= n, rows 1 to M of B contain the
86 * minimum norm solution vectors;
87 * if TRANS = 'T' and m < n, rows 1 to M of B contain the
88 * least squares solution vectors; the residual sum of squares
89 * for the solution in each column is given by the sum of
90 * squares of elements M+1 to N in that column.
91 *
92 * LDB (input) INTEGER
93 * The leading dimension of the array B. LDB >= MAX(1,M,N).
94 *
95 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
96 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
97 *
98 * LWORK (input) INTEGER
99 * The dimension of the array WORK.
100 * LWORK >= max( 1, MN + max( MN, NRHS ) ).
101 * For optimal performance,
102 * LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
103 * where MN = min(M,N) and NB is the optimum block size.
104 *
105 * If LWORK = -1, then a workspace query is assumed; the routine
106 * only calculates the optimal size of the WORK array, returns
107 * this value as the first entry of the WORK array, and no error
108 * message related to LWORK is issued by XERBLA.
109 *
110 * INFO (output) INTEGER
111 * = 0: successful exit
112 * < 0: if INFO = -i, the i-th argument had an illegal value
113 * > 0: if INFO = i, the i-th diagonal element of the
114 * triangular factor of A is zero, so that A does not have
115 * full rank; the least squares solution could not be
116 * computed.
117 *
118 * =====================================================================
119 *
120 * .. Parameters ..
121 DOUBLE PRECISION ZERO, ONE
122 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
123 * ..
124 * .. Local Scalars ..
125 LOGICAL LQUERY, TPSD
126 INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
127 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
128 * ..
129 * .. Local Arrays ..
130 DOUBLE PRECISION RWORK( 1 )
131 * ..
132 * .. External Functions ..
133 LOGICAL LSAME
134 INTEGER ILAENV
135 DOUBLE PRECISION DLAMCH, DLANGE
136 EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
137 * ..
138 * .. External Subroutines ..
139 EXTERNAL DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR,
140 $ DTRTRS, XERBLA
141 * ..
142 * .. Intrinsic Functions ..
143 INTRINSIC DBLE, MAX, MIN
144 * ..
145 * .. Executable Statements ..
146 *
147 * Test the input arguments.
148 *
149 INFO = 0
150 MN = MIN( M, N )
151 LQUERY = ( LWORK.EQ.-1 )
152 IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN
153 INFO = -1
154 ELSE IF( M.LT.0 ) THEN
155 INFO = -2
156 ELSE IF( N.LT.0 ) THEN
157 INFO = -3
158 ELSE IF( NRHS.LT.0 ) THEN
159 INFO = -4
160 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
161 INFO = -6
162 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
163 INFO = -8
164 ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
165 $ THEN
166 INFO = -10
167 END IF
168 *
169 * Figure out optimal block size
170 *
171 IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
172 *
173 TPSD = .TRUE.
174 IF( LSAME( TRANS, 'N' ) )
175 $ TPSD = .FALSE.
176 *
177 IF( M.GE.N ) THEN
178 NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
179 IF( TPSD ) THEN
180 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N,
181 $ -1 ) )
182 ELSE
183 NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N,
184 $ -1 ) )
185 END IF
186 ELSE
187 NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
188 IF( TPSD ) THEN
189 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M,
190 $ -1 ) )
191 ELSE
192 NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M,
193 $ -1 ) )
194 END IF
195 END IF
196 *
197 WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
198 WORK( 1 ) = DBLE( WSIZE )
199 *
200 END IF
201 *
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'DGELS ', -INFO )
204 RETURN
205 ELSE IF( LQUERY ) THEN
206 RETURN
207 END IF
208 *
209 * Quick return if possible
210 *
211 IF( MIN( M, N, NRHS ).EQ.0 ) THEN
212 CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
213 RETURN
214 END IF
215 *
216 * Get machine parameters
217 *
218 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
219 BIGNUM = ONE / SMLNUM
220 CALL DLABAD( SMLNUM, BIGNUM )
221 *
222 * Scale A, B if max element outside range [SMLNUM,BIGNUM]
223 *
224 ANRM = DLANGE( 'M', M, N, A, LDA, RWORK )
225 IASCL = 0
226 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
227 *
228 * Scale matrix norm up to SMLNUM
229 *
230 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
231 IASCL = 1
232 ELSE IF( ANRM.GT.BIGNUM ) THEN
233 *
234 * Scale matrix norm down to BIGNUM
235 *
236 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
237 IASCL = 2
238 ELSE IF( ANRM.EQ.ZERO ) THEN
239 *
240 * Matrix all zero. Return zero solution.
241 *
242 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
243 GO TO 50
244 END IF
245 *
246 BROW = M
247 IF( TPSD )
248 $ BROW = N
249 BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
250 IBSCL = 0
251 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
252 *
253 * Scale matrix norm up to SMLNUM
254 *
255 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
256 $ INFO )
257 IBSCL = 1
258 ELSE IF( BNRM.GT.BIGNUM ) THEN
259 *
260 * Scale matrix norm down to BIGNUM
261 *
262 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
263 $ INFO )
264 IBSCL = 2
265 END IF
266 *
267 IF( M.GE.N ) THEN
268 *
269 * compute QR factorization of A
270 *
271 CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
272 $ INFO )
273 *
274 * workspace at least N, optimally N*NB
275 *
276 IF( .NOT.TPSD ) THEN
277 *
278 * Least-Squares Problem min || A * X - B ||
279 *
280 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
281 *
282 CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA,
283 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
284 $ INFO )
285 *
286 * workspace at least NRHS, optimally NRHS*NB
287 *
288 * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
289 *
290 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
291 $ A, LDA, B, LDB, INFO )
292 *
293 IF( INFO.GT.0 ) THEN
294 RETURN
295 END IF
296 *
297 SCLLEN = N
298 *
299 ELSE
300 *
301 * Overdetermined system of equations A**T * X = B
302 *
303 * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
304 *
305 CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS,
306 $ A, LDA, B, LDB, INFO )
307 *
308 IF( INFO.GT.0 ) THEN
309 RETURN
310 END IF
311 *
312 * B(N+1:M,1:NRHS) = ZERO
313 *
314 DO 20 J = 1, NRHS
315 DO 10 I = N + 1, M
316 B( I, J ) = ZERO
317 10 CONTINUE
318 20 CONTINUE
319 *
320 * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
321 *
322 CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
323 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
324 $ INFO )
325 *
326 * workspace at least NRHS, optimally NRHS*NB
327 *
328 SCLLEN = M
329 *
330 END IF
331 *
332 ELSE
333 *
334 * Compute LQ factorization of A
335 *
336 CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
337 $ INFO )
338 *
339 * workspace at least M, optimally M*NB.
340 *
341 IF( .NOT.TPSD ) THEN
342 *
343 * underdetermined system of equations A * X = B
344 *
345 * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
346 *
347 CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
348 $ A, LDA, B, LDB, INFO )
349 *
350 IF( INFO.GT.0 ) THEN
351 RETURN
352 END IF
353 *
354 * B(M+1:N,1:NRHS) = 0
355 *
356 DO 40 J = 1, NRHS
357 DO 30 I = M + 1, N
358 B( I, J ) = ZERO
359 30 CONTINUE
360 40 CONTINUE
361 *
362 * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
363 *
364 CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA,
365 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
366 $ INFO )
367 *
368 * workspace at least NRHS, optimally NRHS*NB
369 *
370 SCLLEN = N
371 *
372 ELSE
373 *
374 * overdetermined system min || A**T * X - B ||
375 *
376 * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
377 *
378 CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
379 $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
380 $ INFO )
381 *
382 * workspace at least NRHS, optimally NRHS*NB
383 *
384 * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
385 *
386 CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
387 $ A, LDA, B, LDB, INFO )
388 *
389 IF( INFO.GT.0 ) THEN
390 RETURN
391 END IF
392 *
393 SCLLEN = M
394 *
395 END IF
396 *
397 END IF
398 *
399 * Undo scaling
400 *
401 IF( IASCL.EQ.1 ) THEN
402 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
403 $ INFO )
404 ELSE IF( IASCL.EQ.2 ) THEN
405 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
406 $ INFO )
407 END IF
408 IF( IBSCL.EQ.1 ) THEN
409 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
410 $ INFO )
411 ELSE IF( IBSCL.EQ.2 ) THEN
412 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
413 $ INFO )
414 END IF
415 *
416 50 CONTINUE
417 WORK( 1 ) = DBLE( WSIZE )
418 *
419 RETURN
420 *
421 * End of DGELS
422 *
423 END