1 *> \brief <b> DGELSY solves overdetermined or underdetermined systems for GE matrices</b>
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download DGELSY + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsy.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsy.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsy.f">
21 * SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
24 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26 * DOUBLE PRECISION RCOND
28 * .. Array Arguments ..
30 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
39 *> DGELSY computes the minimum-norm solution to a real linear least
41 *> minimize || A * X - B ||
42 *> using a complete orthogonal factorization of A. A is an M-by-N
43 *> matrix which may be rank-deficient.
45 *> Several right hand side vectors b and solution vectors x can be
46 *> handled in a single call; they are stored as the columns of the
47 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
50 *> The routine first computes a QR factorization with column pivoting:
51 *> A * P = Q * [ R11 R12 ]
53 *> with R11 defined as the largest leading submatrix whose estimated
54 *> condition number is less than 1/RCOND. The order of R11, RANK,
55 *> is the effective rank of A.
57 *> Then, R22 is considered to be negligible, and R12 is annihilated
58 *> by orthogonal transformations from the right, arriving at the
59 *> complete orthogonal factorization:
60 *> A * P = Q * [ T11 0 ] * Z
62 *> The minimum-norm solution is then
63 *> X = P * Z**T [ inv(T11)*Q1**T*B ]
65 *> where Q1 consists of the first RANK columns of Q.
67 *> This routine is basically identical to the original xGELSX except
69 *> o The call to the subroutine xGEQPF has been substituted by the
70 *> the call to the subroutine xGEQP3. This subroutine is a Blas-3
71 *> version of the QR factorization with column pivoting.
72 *> o Matrix B (the right hand side) is updated with Blas-3.
73 *> o The permutation of matrix B (the right hand side) is faster and
83 *> The number of rows of the matrix A. M >= 0.
89 *> The number of columns of the matrix A. N >= 0.
95 *> The number of right hand sides, i.e., the number of
96 *> columns of matrices B and X. NRHS >= 0.
101 *> A is DOUBLE PRECISION array, dimension (LDA,N)
102 *> On entry, the M-by-N matrix A.
103 *> On exit, A has been overwritten by details of its
104 *> complete orthogonal factorization.
110 *> The leading dimension of the array A. LDA >= max(1,M).
115 *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
116 *> On entry, the M-by-NRHS right hand side matrix B.
117 *> On exit, the N-by-NRHS solution matrix X.
123 *> The leading dimension of the array B. LDB >= max(1,M,N).
126 *> \param[in,out] JPVT
128 *> JPVT is INTEGER array, dimension (N)
129 *> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
130 *> to the front of AP, otherwise column i is a free column.
131 *> On exit, if JPVT(i) = k, then the i-th column of AP
132 *> was the k-th column of A.
137 *> RCOND is DOUBLE PRECISION
138 *> RCOND is used to determine the effective rank of A, which
139 *> is defined as the order of the largest leading triangular
140 *> submatrix R11 in the QR factorization with pivoting of A,
141 *> whose estimated condition number < 1/RCOND.
147 *> The effective rank of A, i.e., the order of the submatrix
148 *> R11. This is the same as the order of the submatrix T11
149 *> in the complete orthogonal factorization of A.
154 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
155 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
161 *> The dimension of the array WORK.
162 *> The unblocked strategy requires that:
163 *> LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
164 *> where MN = min( M, N ).
165 *> The block algorithm requires that:
166 *> LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
167 *> where NB is an upper bound on the blocksize returned
168 *> by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
171 *> If LWORK = -1, then a workspace query is assumed; the routine
172 *> only calculates the optimal size of the WORK array, returns
173 *> this value as the first entry of the WORK array, and no error
174 *> message related to LWORK is issued by XERBLA.
180 *> = 0: successful exit
181 *> < 0: If INFO = -i, the i-th argument had an illegal value.
187 *> \author Univ. of Tennessee
188 *> \author Univ. of California Berkeley
189 *> \author Univ. of Colorado Denver
192 *> \date November 2011
194 *> \ingroup doubleGEsolve
196 *> \par Contributors:
199 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
200 *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
201 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
203 * =====================================================================
204 SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205 $ WORK, LWORK, INFO )
207 * -- LAPACK driver routine (version 3.4.0) --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212 * .. Scalar Arguments ..
213 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
214 DOUBLE PRECISION RCOND
216 * .. Array Arguments ..
218 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
221 * =====================================================================
225 PARAMETER ( IMAX = 1, IMIN = 2 )
226 DOUBLE PRECISION ZERO, ONE
227 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
229 * .. Local Scalars ..
231 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
232 $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4
233 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
234 $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
236 * .. External Functions ..
238 DOUBLE PRECISION DLAMCH, DLANGE
239 EXTERNAL ILAENV, DLAMCH, DLANGE
241 * .. External Subroutines ..
242 EXTERNAL DCOPY, DGEQP3, DLABAD, DLAIC1, DLASCL, DLASET,
243 $ DORMQR, DORMRZ, DTRSM, DTZRZF, XERBLA
245 * .. Intrinsic Functions ..
246 INTRINSIC ABS, MAX, MIN
248 * .. Executable Statements ..
254 * Test the input arguments.
257 LQUERY = ( LWORK.EQ.-1 )
260 ELSE IF( N.LT.0 ) THEN
262 ELSE IF( NRHS.LT.0 ) THEN
264 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
266 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
270 * Figure out optimal block size
273 IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
277 NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
278 NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
279 NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, NRHS, -1 )
280 NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, NRHS, -1 )
281 NB = MAX( NB1, NB2, NB3, NB4 )
282 LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS )
283 LWKOPT = MAX( LWKMIN,
284 $ MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS )
288 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
294 CALL XERBLA( 'DGELSY', -INFO )
296 ELSE IF( LQUERY ) THEN
300 * Quick return if possible
302 IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
307 * Get machine parameters
309 SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
310 BIGNUM = ONE / SMLNUM
311 CALL DLABAD( SMLNUM, BIGNUM )
313 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
315 ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
317 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
319 * Scale matrix norm up to SMLNUM
321 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
323 ELSE IF( ANRM.GT.BIGNUM ) THEN
325 * Scale matrix norm down to BIGNUM
327 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
329 ELSE IF( ANRM.EQ.ZERO ) THEN
331 * Matrix all zero. Return zero solution.
333 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
338 BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
340 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
342 * Scale matrix norm up to SMLNUM
344 CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
346 ELSE IF( BNRM.GT.BIGNUM ) THEN
348 * Scale matrix norm down to BIGNUM
350 CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
354 * Compute QR factorization with column pivoting of A:
357 CALL DGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
359 WSIZE = MN + WORK( MN+1 )
361 * workspace: MN+2*N+NB*(N+1).
362 * Details of Householder rotations stored in WORK(1:MN).
364 * Determine RANK using incremental condition estimation
368 SMAX = ABS( A( 1, 1 ) )
370 IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
372 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
379 IF( RANK.LT.MN ) THEN
381 CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
382 $ A( I, I ), SMINPR, S1, C1 )
383 CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
384 $ A( I, I ), SMAXPR, S2, C2 )
386 IF( SMAXPR*RCOND.LE.SMINPR ) THEN
388 WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
389 WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
391 WORK( ISMIN+RANK ) = C1
392 WORK( ISMAX+RANK ) = C2
402 * Logically partition R = [ R11 R12 ]
404 * where R11 = R(1:RANK,1:RANK)
406 * [R11,R12] = [ T11, 0 ] * Y
409 $ CALL DTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
413 * Details of Householder rotations stored in WORK(MN+1:2*MN)
415 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
417 CALL DORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
418 $ B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
419 WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) )
421 * workspace: 2*MN+NB*NRHS.
423 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
425 CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
426 $ NRHS, ONE, A, LDA, B, LDB )
429 DO 30 I = RANK + 1, N
434 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
437 CALL DORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A,
438 $ LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ),
442 * workspace: 2*MN+NRHS.
444 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
448 WORK( JPVT( I ) ) = B( I, J )
450 CALL DCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
457 IF( IASCL.EQ.1 ) THEN
458 CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
459 CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
461 ELSE IF( IASCL.EQ.2 ) THEN
462 CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
463 CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
466 IF( IBSCL.EQ.1 ) THEN
467 CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
468 ELSE IF( IBSCL.EQ.2 ) THEN
469 CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )