LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgelst.f
Go to the documentation of this file.
1*> \brief <b> DGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGELST + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelst.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelst.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelst.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER TRANS
26* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DGELST solves overdetermined or underdetermined real linear systems
39*> involving an M-by-N matrix A, or its transpose, using a QR or LQ
40*> factorization of A with compact WY representation of Q.
41*> It is assumed that A has full rank.
42*>
43*> The following options are provided:
44*>
45*> 1. If TRANS = 'N' and m >= n: find the least squares solution of
46*> an overdetermined system, i.e., solve the least squares problem
47*> minimize || B - A*X ||.
48*>
49*> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
50*> an underdetermined system A * X = B.
51*>
52*> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
53*> an underdetermined system A**T * X = B.
54*>
55*> 4. If TRANS = 'T' and m < n: find the least squares solution of
56*> an overdetermined system, i.e., solve the least squares problem
57*> minimize || B - A**T * X ||.
58*>
59*> Several right hand side vectors b and solution vectors x can be
60*> handled in a single call; they are stored as the columns of the
61*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
62*> matrix X.
63*> \endverbatim
64*
65* Arguments:
66* ==========
67*
68*> \param[in] TRANS
69*> \verbatim
70*> TRANS is CHARACTER*1
71*> = 'N': the linear system involves A;
72*> = 'T': the linear system involves A**T.
73*> \endverbatim
74*>
75*> \param[in] M
76*> \verbatim
77*> M is INTEGER
78*> The number of rows of the matrix A. M >= 0.
79*> \endverbatim
80*>
81*> \param[in] N
82*> \verbatim
83*> N is INTEGER
84*> The number of columns of the matrix A. N >= 0.
85*> \endverbatim
86*>
87*> \param[in] NRHS
88*> \verbatim
89*> NRHS is INTEGER
90*> The number of right hand sides, i.e., the number of
91*> columns of the matrices B and X. NRHS >=0.
92*> \endverbatim
93*>
94*> \param[in,out] A
95*> \verbatim
96*> A is DOUBLE PRECISION array, dimension (LDA,N)
97*> On entry, the M-by-N matrix A.
98*> On exit,
99*> if M >= N, A is overwritten by details of its QR
100*> factorization as returned by DGEQRT;
101*> if M < N, A is overwritten by details of its LQ
102*> factorization as returned by DGELQT.
103*> \endverbatim
104*>
105*> \param[in] LDA
106*> \verbatim
107*> LDA is INTEGER
108*> The leading dimension of the array A. LDA >= max(1,M).
109*> \endverbatim
110*>
111*> \param[in,out] B
112*> \verbatim
113*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
114*> On entry, the matrix B of right hand side vectors, stored
115*> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
116*> if TRANS = 'T'.
117*> On exit, if INFO = 0, B is overwritten by the solution
118*> vectors, stored columnwise:
119*> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
120*> squares solution vectors; the residual sum of squares for the
121*> solution in each column is given by the sum of squares of
122*> elements N+1 to M in that column;
123*> if TRANS = 'N' and m < n, rows 1 to N of B contain the
124*> minimum norm solution vectors;
125*> if TRANS = 'T' and m >= n, rows 1 to M of B contain the
126*> minimum norm solution vectors;
127*> if TRANS = 'T' and m < n, rows 1 to M of B contain the
128*> least squares solution vectors; the residual sum of squares
129*> for the solution in each column is given by the sum of
130*> squares of elements M+1 to N in that column.
131*> \endverbatim
132*>
133*> \param[in] LDB
134*> \verbatim
135*> LDB is INTEGER
136*> The leading dimension of the array B. LDB >= MAX(1,M,N).
137*> \endverbatim
138*>
139*> \param[out] WORK
140*> \verbatim
141*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
142*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143*> \endverbatim
144*>
145*> \param[in] LWORK
146*> \verbatim
147*> LWORK is INTEGER
148*> The dimension of the array WORK.
149*> LWORK >= max( 1, MN + max( MN, NRHS ) ).
150*> For optimal performance,
151*> LWORK >= max( 1, (MN + max( MN, NRHS ))*NB ).
152*> where MN = min(M,N) and NB is the optimum block size.
153*>
154*> If LWORK = -1, then a workspace query is assumed; the routine
155*> only calculates the optimal size of the WORK array, returns
156*> this value as the first entry of the WORK array, and no error
157*> message related to LWORK is issued by XERBLA.
158*> \endverbatim
159*>
160*> \param[out] INFO
161*> \verbatim
162*> INFO is INTEGER
163*> = 0: successful exit
164*> < 0: if INFO = -i, the i-th argument had an illegal value
165*> > 0: if INFO = i, the i-th diagonal element of the
166*> triangular factor of A is zero, so that A does not have
167*> full rank; the least squares solution could not be
168*> computed.
169*> \endverbatim
170*
171* Authors:
172* ========
173*
174*> \author Univ. of Tennessee
175*> \author Univ. of California Berkeley
176*> \author Univ. of Colorado Denver
177*> \author NAG Ltd.
178*
179*> \ingroup gelst
180*
181*> \par Contributors:
182* ==================
183*>
184*> \verbatim
185*>
186*> November 2022, Igor Kozachenko,
187*> Computer Science Division,
188*> University of California, Berkeley
189*> \endverbatim
190*
191* =====================================================================
192 SUBROUTINE dgelst( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
193 $ INFO )
194*
195* -- LAPACK driver routine --
196* -- LAPACK is a software package provided by Univ. of Tennessee, --
197* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198*
199* .. Scalar Arguments ..
200 CHARACTER TRANS
201 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
202* ..
203* .. Array Arguments ..
204 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 DOUBLE PRECISION ZERO, ONE
211 parameter( zero = 0.0d+0, one = 1.0d+0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL LQUERY, TPSD
215 INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
216 $ nb, nbmin, scllen
217 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
218* ..
219* .. Local Arrays ..
220 DOUBLE PRECISION RWORK( 1 )
221* ..
222* .. External Functions ..
223 LOGICAL LSAME
224 INTEGER ILAENV
225 DOUBLE PRECISION DLAMCH, DLANGE
226 EXTERNAL lsame, ilaenv, dlamch, dlange
227* ..
228* .. External Subroutines ..
229 EXTERNAL dgelqt, dgeqrt, dgemlqt, dgemqrt, dlascl,
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC dble, max, min
234* ..
235* .. Executable Statements ..
236*
237* Test the input arguments.
238*
239 info = 0
240 mn = min( m, n )
241 lquery = ( lwork.EQ.-1 )
242 IF( .NOT.( lsame( trans, 'N' ) .OR. lsame( trans, 'T' ) ) ) THEN
243 info = -1
244 ELSE IF( m.LT.0 ) THEN
245 info = -2
246 ELSE IF( n.LT.0 ) THEN
247 info = -3
248 ELSE IF( nrhs.LT.0 ) THEN
249 info = -4
250 ELSE IF( lda.LT.max( 1, m ) ) THEN
251 info = -6
252 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
253 info = -8
254 ELSE IF( lwork.LT.max( 1, mn+max( mn, nrhs ) ) .AND. .NOT.lquery )
255 $ THEN
256 info = -10
257 END IF
258*
259* Figure out optimal block size and optimal workspace size
260*
261 IF( info.EQ.0 .OR. info.EQ.-10 ) THEN
262*
263 tpsd = .true.
264 IF( lsame( trans, 'N' ) )
265 $ tpsd = .false.
266*
267 nb = ilaenv( 1, 'DGELST', ' ', m, n, -1, -1 )
268*
269 mnnrhs = max( mn, nrhs )
270 lwopt = max( 1, (mn+mnnrhs)*nb )
271 work( 1 ) = dble( lwopt )
272*
273 END IF
274*
275 IF( info.NE.0 ) THEN
276 CALL xerbla( 'DGELST ', -info )
277 RETURN
278 ELSE IF( lquery ) THEN
279 RETURN
280 END IF
281*
282* Quick return if possible
283*
284 IF( min( m, n, nrhs ).EQ.0 ) THEN
285 CALL dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
286 work( 1 ) = dble( lwopt )
287 RETURN
288 END IF
289*
290* *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
291*
292 IF( nb.GT.mn ) nb = mn
293*
294* Determine the block size from the supplied LWORK
295* ( at this stage we know that LWORK >= (minimum required workspace,
296* but it may be less than optimal)
297*
298 nb = min( nb, lwork/( mn + mnnrhs ) )
299*
300* The minimum value of NB, when blocked code is used
301*
302 nbmin = max( 2, ilaenv( 2, 'DGELST', ' ', m, n, -1, -1 ) )
303*
304 IF( nb.LT.nbmin ) THEN
305 nb = 1
306 END IF
307*
308* Get machine parameters
309*
310 smlnum = dlamch( 'S' ) / dlamch( 'P' )
311 bignum = one / smlnum
312*
313* Scale A, B if max element outside range [SMLNUM,BIGNUM]
314*
315 anrm = dlange( 'M', m, n, a, lda, rwork )
316 iascl = 0
317 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
318*
319* Scale matrix norm up to SMLNUM
320*
321 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
322 iascl = 1
323 ELSE IF( anrm.GT.bignum ) THEN
324*
325* Scale matrix norm down to BIGNUM
326*
327 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
328 iascl = 2
329 ELSE IF( anrm.EQ.zero ) THEN
330*
331* Matrix all zero. Return zero solution.
332*
333 CALL dlaset( 'Full', max( m, n ), nrhs, zero, zero, b, ldb )
334 work( 1 ) = dble( lwopt )
335 RETURN
336 END IF
337*
338 brow = m
339 IF( tpsd )
340 $ brow = n
341 bnrm = dlange( 'M', brow, nrhs, b, ldb, rwork )
342 ibscl = 0
343 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
344*
345* Scale matrix norm up to SMLNUM
346*
347 CALL dlascl( 'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
348 $ info )
349 ibscl = 1
350 ELSE IF( bnrm.GT.bignum ) THEN
351*
352* Scale matrix norm down to BIGNUM
353*
354 CALL dlascl( 'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
355 $ info )
356 ibscl = 2
357 END IF
358*
359 IF( m.GE.n ) THEN
360*
361* M > N:
362* Compute the blocked QR factorization of A,
363* using the compact WY representation of Q,
364* workspace at least N, optimally N*NB.
365*
366 CALL dgeqrt( m, n, nb, a, lda, work( 1 ), nb,
367 $ work( mn*nb+1 ), info )
368*
369 IF( .NOT.tpsd ) THEN
370*
371* M > N, A is not transposed:
372* Overdetermined system of equations,
373* least-squares problem, min || A * X - B ||.
374*
375* Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
376* using the compact WY representation of Q,
377* workspace at least NRHS, optimally NRHS*NB.
378*
379 CALL dgemqrt( 'Left', 'Transpose', m, nrhs, n, nb, a, lda,
380 $ work( 1 ), nb, b, ldb, work( mn*nb+1 ),
381 $ info )
382*
383* Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
384*
385 CALL dtrtrs( 'Upper', 'No transpose', 'Non-unit', n, nrhs,
386 $ a, lda, b, ldb, info )
387*
388 IF( info.GT.0 ) THEN
389 RETURN
390 END IF
391*
392 scllen = n
393*
394 ELSE
395*
396* M > N, A is transposed:
397* Underdetermined system of equations,
398* minimum norm solution of A**T * X = B.
399*
400* Compute B := inv(R**T) * B in two row blocks of B.
401*
402* Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
403*
404 CALL dtrtrs( 'Upper', 'Transpose', 'Non-unit', n, nrhs,
405 $ a, lda, b, ldb, info )
406*
407 IF( info.GT.0 ) THEN
408 RETURN
409 END IF
410*
411* Block 2: Zero out all rows below the N-th row in B:
412* B(N+1:M,1:NRHS) = ZERO
413*
414 DO j = 1, nrhs
415 DO i = n + 1, m
416 b( i, j ) = zero
417 END DO
418 END DO
419*
420* Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
421* using the compact WY representation of Q,
422* workspace at least NRHS, optimally NRHS*NB.
423*
424 CALL dgemqrt( 'Left', 'No transpose', m, nrhs, n, nb,
425 $ a, lda, work( 1 ), nb, b, ldb,
426 $ work( mn*nb+1 ), info )
427*
428 scllen = m
429*
430 END IF
431*
432 ELSE
433*
434* M < N:
435* Compute the blocked LQ factorization of A,
436* using the compact WY representation of Q,
437* workspace at least M, optimally M*NB.
438*
439 CALL dgelqt( m, n, nb, a, lda, work( 1 ), nb,
440 $ work( mn*nb+1 ), info )
441*
442 IF( .NOT.tpsd ) THEN
443*
444* M < N, A is not transposed:
445* Underdetermined system of equations,
446* minimum norm solution of A * X = B.
447*
448* Compute B := inv(L) * B in two row blocks of B.
449*
450* Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
451*
452 CALL dtrtrs( 'Lower', 'No transpose', 'Non-unit', m, nrhs,
453 $ a, lda, b, ldb, info )
454*
455 IF( info.GT.0 ) THEN
456 RETURN
457 END IF
458*
459* Block 2: Zero out all rows below the M-th row in B:
460* B(M+1:N,1:NRHS) = ZERO
461*
462 DO j = 1, nrhs
463 DO i = m + 1, n
464 b( i, j ) = zero
465 END DO
466 END DO
467*
468* Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
469* using the compact WY representation of Q,
470* workspace at least NRHS, optimally NRHS*NB.
471*
472 CALL dgemlqt( 'Left', 'Transpose', n, nrhs, m, nb, a, lda,
473 $ work( 1 ), nb, b, ldb,
474 $ work( mn*nb+1 ), info )
475*
476 scllen = n
477*
478 ELSE
479*
480* M < N, A is transposed:
481* Overdetermined system of equations,
482* least-squares problem, min || A**T * X - B ||.
483*
484* Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
485* using the compact WY representation of Q,
486* workspace at least NRHS, optimally NRHS*NB.
487*
488 CALL dgemlqt( 'Left', 'No transpose', n, nrhs, m, nb,
489 $ a, lda, work( 1 ), nb, b, ldb,
490 $ work( mn*nb+1), info )
491*
492* Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
493*
494 CALL dtrtrs( 'Lower', 'Transpose', 'Non-unit', m, nrhs,
495 $ a, lda, b, ldb, info )
496*
497 IF( info.GT.0 ) THEN
498 RETURN
499 END IF
500*
501 scllen = m
502*
503 END IF
504*
505 END IF
506*
507* Undo scaling
508*
509 IF( iascl.EQ.1 ) THEN
510 CALL dlascl( 'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
511 $ info )
512 ELSE IF( iascl.EQ.2 ) THEN
513 CALL dlascl( 'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
514 $ info )
515 END IF
516 IF( ibscl.EQ.1 ) THEN
517 CALL dlascl( 'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
518 $ info )
519 ELSE IF( ibscl.EQ.2 ) THEN
520 CALL dlascl( 'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
521 $ info )
522 END IF
523*
524 work( 1 ) = dble( lwopt )
525*
526 RETURN
527*
528* End of DGELST
529*
530 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgelqt(m, n, mb, a, lda, t, ldt, work, info)
DGELQT
Definition dgelqt.f:139
subroutine dgelst(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
DGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization ...
Definition dgelst.f:194
subroutine dgemlqt(side, trans, m, n, k, mb, v, ldv, t, ldt, c, ldc, work, info)
DGEMLQT
Definition dgemlqt.f:168
subroutine dgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
DGEMQRT
Definition dgemqrt.f:168
subroutine dgeqrt(m, n, nb, a, lda, t, ldt, work, info)
DGEQRT
Definition dgeqrt.f:141
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dtrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
DTRTRS
Definition dtrtrs.f:140