#include "ESMF_LapackBlas.inc" !> \brief <b> DGELSY solves overdetermined or underdetermined systems for GE matrices</b> ! ! =========== DOCUMENTATION =========== ! ! Online html documentation available at ! http://www.netlib.org/lapack/explore-html/ ! !> \htmlonly !> Download DGELSY + dependencies !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsy.f"> !> [TGZ]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsy.f"> !> [ZIP]</a> !> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsy.f"> !> [TXT]</a> !> \endhtmlonly ! ! Definition: ! =========== ! ! SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, ! WORK, LWORK, INFO ) ! ! .. Scalar Arguments .. ! INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK ! DOUBLE PRECISION RCOND ! .. ! .. Array Arguments .. ! INTEGER JPVT( * ) ! DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) ! .. ! ! !> \par Purpose: ! ============= !> !> \verbatim !> !> DGELSY computes the minimum-norm solution to a real linear least !> squares problem: !> minimize || A * X - B || !> using a complete orthogonal factorization of A. A is an M-by-N !> matrix which may be rank-deficient. !> !> Several right hand side vectors b and solution vectors x can be !> handled in a single call; they are stored as the columns of the !> M-by-NRHS right hand side matrix B and the N-by-NRHS solution !> matrix X. !> !> The routine first computes a QR factorization with column pivoting: !> A * P = Q * [ R11 R12 ] !> [ 0 R22 ] !> with R11 defined as the largest leading submatrix whose estimated !> condition number is less than 1/RCOND. The order of R11, RANK, !> is the effective rank of A. !> !> Then, R22 is considered to be negligible, and R12 is annihilated !> by orthogonal transformations from the right, arriving at the !> complete orthogonal factorization: !> A * P = Q * [ T11 0 ] * Z !> [ 0 0 ] !> The minimum-norm solution is then !> X = P * Z**T [ inv(T11)*Q1**T*B ] !> [ 0 ] !> where Q1 consists of the first RANK columns of Q. !> !> This routine is basically identical to the original xGELSX except !> three differences: !> o The call to the subroutine xGEQPF has been substituted by the !> the call to the subroutine xGEQP3. This subroutine is a Blas-3 !> version of the QR factorization with column pivoting. !> o Matrix B (the right hand side) is updated with Blas-3. !> o The permutation of matrix B (the right hand side) is faster and !> more simple. !> \endverbatim ! ! Arguments: ! ========== ! !> \param[in] M !> \verbatim !> M is INTEGER !> The number of rows of the matrix A. M >= 0. !> \endverbatim !> !> \param[in] N !> \verbatim !> N is INTEGER !> The number of columns of the matrix A. N >= 0. !> \endverbatim !> !> \param[in] NRHS !> \verbatim !> NRHS is INTEGER !> The number of right hand sides, i.e., the number of !> columns of matrices B and X. NRHS >= 0. !> \endverbatim !> !> \param[in,out] A !> \verbatim !> A is DOUBLE PRECISION array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> On exit, A has been overwritten by details of its !> complete orthogonal factorization. !> \endverbatim !> !> \param[in] LDA !> \verbatim !> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !> \endverbatim !> !> \param[in,out] B !> \verbatim !> B is DOUBLE PRECISION array, dimension (LDB,NRHS) !> On entry, the M-by-NRHS right hand side matrix B. !> On exit, the N-by-NRHS solution matrix X. !> \endverbatim !> !> \param[in] LDB !> \verbatim !> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,M,N). !> \endverbatim !> !> \param[in,out] JPVT !> \verbatim !> JPVT is INTEGER array, dimension (N) !> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted !> to the front of AP, otherwise column i is a free column. !> On exit, if JPVT(i) = k, then the i-th column of AP !> was the k-th column of A. !> \endverbatim !> !> \param[in] RCOND !> \verbatim !> RCOND is DOUBLE PRECISION !> RCOND is used to determine the effective rank of A, which !> is defined as the order of the largest leading triangular !> submatrix R11 in the QR factorization with pivoting of A, !> whose estimated condition number < 1/RCOND. !> \endverbatim !> !> \param[out] RANK !> \verbatim !> RANK is INTEGER !> The effective rank of A, i.e., the order of the submatrix !> R11. This is the same as the order of the submatrix T11 !> in the complete orthogonal factorization of A. !> \endverbatim !> !> \param[out] WORK !> \verbatim !> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) !> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. !> \endverbatim !> !> \param[in] LWORK !> \verbatim !> LWORK is INTEGER !> The dimension of the array WORK. !> The unblocked strategy requires that: !> LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), !> where MN = min( M, N ). !> The block algorithm requires that: !> LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), !> where NB is an upper bound on the blocksize returned !> by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR, !> and DORMRZ. !> !> If LWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal size of the WORK array, returns !> this value as the first entry of the WORK array, and no error !> message related to LWORK is issued by XERBLA. !> \endverbatim !> !> \param[out] INFO !> \verbatim !> INFO is INTEGER !> = 0: successful exit !> < 0: If INFO = -i, the i-th argument had an illegal value. !> \endverbatim ! ! Authors: ! ======== ! !> \author Univ. of Tennessee !> \author Univ. of California Berkeley !> \author Univ. of Colorado Denver !> \author NAG Ltd. ! !> \date November 2011 ! !> \ingroup doubleGEsolve ! !> \par Contributors: ! ================== !> !> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n !> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n !> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n !> ! ===================================================================== SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, & & WORK, LWORK, INFO ) ! ! -- LAPACK driver routine (version 3.4.0) -- ! -- LAPACK is a software package provided by Univ. of Tennessee, -- ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- ! November 2011 ! ! .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK DOUBLE PRECISION RCOND ! .. ! .. Array Arguments .. INTEGER JPVT( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) ! .. ! ! ===================================================================== ! ! .. Parameters .. INTEGER IMAX, IMIN PARAMETER ( IMAX = 1, IMIN = 2 ) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) ! .. ! .. Local Scalars .. LOGICAL LQUERY INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN, & & LWKOPT, MN, NB, NB1, NB2, NB3, NB4 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX, & & SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE ! .. ! .. External Functions .. INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLANGE EXTERNAL ILAENV, DLAMCH, DLANGE ! .. ! .. External Subroutines .. EXTERNAL DCOPY, DGEQP3, DLABAD, DLAIC1, DLASCL, DLASET, & & DORMQR, DORMRZ, DTRSM, DTZRZF, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN ! .. ! .. Executable Statements .. ! MN = MIN( M, N ) ISMIN = MN + 1 ISMAX = 2*MN + 1 ! ! Test the input arguments. ! INFO = 0 LQUERY = ( LWORK.EQ.-1 ) IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( NRHS.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -5 ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN INFO = -7 END IF ! ! Figure out optimal block size ! IF( INFO.EQ.0 ) THEN IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN LWKMIN = 1 LWKOPT = 1 ELSE NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 ) NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, NRHS, -1 ) NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, NRHS, -1 ) NB = MAX( NB1, NB2, NB3, NB4 ) LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS ) LWKOPT = MAX( LWKMIN, & & MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS ) END IF WORK( 1 ) = LWKOPT ! IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN INFO = -12 END IF END IF ! IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGELSY', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF ! ! Quick return if possible ! IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN RANK = 0 RETURN END IF ! ! Get machine parameters ! SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) BIGNUM = ONE / SMLNUM CALL DLABAD( SMLNUM, BIGNUM ) ! ! Scale A, B if max entries outside range [SMLNUM,BIGNUM] ! ANRM = DLANGE( 'M', M, N, A, LDA, WORK ) IASCL = 0 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN ! ! Scale matrix norm up to SMLNUM ! CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) IASCL = 1 ELSE IF( ANRM.GT.BIGNUM ) THEN ! ! Scale matrix norm down to BIGNUM ! CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) IASCL = 2 ELSE IF( ANRM.EQ.ZERO ) THEN ! ! Matrix all zero. Return zero solution. ! CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) RANK = 0 GO TO 70 END IF ! BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK ) IBSCL = 0 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN ! ! Scale matrix norm up to SMLNUM ! CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) IBSCL = 1 ELSE IF( BNRM.GT.BIGNUM ) THEN ! ! Scale matrix norm down to BIGNUM ! CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) IBSCL = 2 END IF ! ! Compute QR factorization with column pivoting of A: ! A * P = Q * R ! CALL DGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), & & LWORK-MN, INFO ) WSIZE = MN + WORK( MN+1 ) ! ! workspace: MN+2*N+NB*(N+1). ! Details of Householder rotations stored in WORK(1:MN). ! ! Determine RANK using incremental condition estimation ! WORK( ISMIN ) = ONE WORK( ISMAX ) = ONE SMAX = ABS( A( 1, 1 ) ) SMIN = SMAX IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN RANK = 0 CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB ) GO TO 70 ELSE RANK = 1 END IF ! 10 CONTINUE IF( RANK.LT.MN ) THEN I = RANK + 1 CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), & & A( I, I ), SMINPR, S1, C1 ) CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), & & A( I, I ), SMAXPR, S2, C2 ) ! IF( SMAXPR*RCOND.LE.SMINPR ) THEN DO 20 I = 1, RANK WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 20 CONTINUE WORK( ISMIN+RANK ) = C1 WORK( ISMAX+RANK ) = C2 SMIN = SMINPR SMAX = SMAXPR RANK = RANK + 1 GO TO 10 END IF END IF ! ! workspace: 3*MN. ! ! Logically partition R = [ R11 R12 ] ! [ 0 R22 ] ! where R11 = R(1:RANK,1:RANK) ! ! [R11,R12] = [ T11, 0 ] * Y ! IF( RANK.LT.N ) & & CALL DTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ), & & LWORK-2*MN, INFO ) ! ! workspace: 2*MN. ! Details of Householder rotations stored in WORK(MN+1:2*MN) ! ! B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS) ! CALL DORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ), & & B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) ) ! ! workspace: 2*MN+NB*NRHS. ! ! B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) ! CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, & & NRHS, ONE, A, LDA, B, LDB ) ! DO 40 J = 1, NRHS DO 30 I = RANK + 1, N B( I, J ) = ZERO 30 CONTINUE 40 CONTINUE ! ! B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS) ! IF( RANK.LT.N ) THEN CALL DORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A, & & LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ), & & LWORK-2*MN, INFO ) END IF ! ! workspace: 2*MN+NRHS. ! ! B(1:N,1:NRHS) := P * B(1:N,1:NRHS) ! DO 60 J = 1, NRHS DO 50 I = 1, N WORK( JPVT( I ) ) = B( I, J ) 50 CONTINUE CALL DCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 ) 60 CONTINUE ! ! workspace: N. ! ! Undo scaling ! IF( IASCL.EQ.1 ) THEN CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, & & INFO ) ELSE IF( IASCL.EQ.2 ) THEN CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, & & INFO ) END IF IF( IBSCL.EQ.1 ) THEN CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) ELSE IF( IBSCL.EQ.2 ) THEN CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) END IF ! 70 CONTINUE WORK( 1 ) = LWKOPT ! RETURN ! ! End of DGELSY ! END