349 lines
10 KiB
Fortran
349 lines
10 KiB
Fortran
*> \brief \b DGETSQRHRT
|
|
*
|
|
* =========== DOCUMENTATION ===========
|
|
*
|
|
* Online html documentation available at
|
|
* http://www.netlib.org/lapack/explore-html/
|
|
*
|
|
*> \htmlonly
|
|
*> Download DGETSQRHRT + dependencies
|
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgetsqrhrt.f">
|
|
*> [TGZ]</a>
|
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgetsqrhrt.f">
|
|
*> [ZIP]</a>
|
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgetsqrhrt.f">
|
|
*> [TXT]</a>
|
|
*> \endhtmlonly
|
|
*
|
|
* Definition:
|
|
* ===========
|
|
*
|
|
* SUBROUTINE DGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
|
|
* $ LWORK, INFO )
|
|
* IMPLICIT NONE
|
|
*
|
|
* .. Scalar Arguments ..
|
|
* INTEGER INFO, LDA, LDT, LWORK, M, N, NB1, NB2, MB1
|
|
* ..
|
|
* .. Array Arguments ..
|
|
* DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
|
|
* ..
|
|
*
|
|
*
|
|
*> \par Purpose:
|
|
* =============
|
|
*>
|
|
*> \verbatim
|
|
*>
|
|
*> DGETSQRHRT computes a NB2-sized column blocked QR-factorization
|
|
*> of a real M-by-N matrix A with M >= N,
|
|
*>
|
|
*> A = Q * R.
|
|
*>
|
|
*> The routine uses internally a NB1-sized column blocked and MB1-sized
|
|
*> row blocked TSQR-factorization and perfors the reconstruction
|
|
*> of the Householder vectors from the TSQR output. The routine also
|
|
*> converts the R_tsqr factor from the TSQR-factorization output into
|
|
*> the R factor that corresponds to the Householder QR-factorization,
|
|
*>
|
|
*> A = Q_tsqr * R_tsqr = Q * R.
|
|
*>
|
|
*> The output Q and R factors are stored in the same format as in DGEQRT
|
|
*> (Q is in blocked compact WY-representation). See the documentation
|
|
*> of DGEQRT for more details on the format.
|
|
*> \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. M >= N >= 0.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] MB1
|
|
*> \verbatim
|
|
*> MB1 is INTEGER
|
|
*> The row block size to be used in the blocked TSQR.
|
|
*> MB1 > N.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NB1
|
|
*> \verbatim
|
|
*> NB1 is INTEGER
|
|
*> The column block size to be used in the blocked TSQR.
|
|
*> N >= NB1 >= 1.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NB2
|
|
*> \verbatim
|
|
*> NB2 is INTEGER
|
|
*> The block size to be used in the blocked QR that is
|
|
*> output. NB2 >= 1.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in,out] A
|
|
*> \verbatim
|
|
*> A is DOUBLE PRECISION array, dimension (LDA,N)
|
|
*>
|
|
*> On entry: an M-by-N matrix A.
|
|
*>
|
|
*> On exit:
|
|
*> a) the elements on and above the diagonal
|
|
*> of the array contain the N-by-N upper-triangular
|
|
*> matrix R corresponding to the Householder QR;
|
|
*> b) the elements below the diagonal represent Q by
|
|
*> the columns of blocked V (compact WY-representation).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDA
|
|
*> \verbatim
|
|
*> LDA is INTEGER
|
|
*> The leading dimension of the array A. LDA >= max(1,M).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] T
|
|
*> \verbatim
|
|
*> T is DOUBLE PRECISION array, dimension (LDT,N))
|
|
*> The upper triangular block reflectors stored in compact form
|
|
*> as a sequence of upper triangular blocks.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDT
|
|
*> \verbatim
|
|
*> LDT is INTEGER
|
|
*> The leading dimension of the array T. LDT >= NB2.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WORK
|
|
*> \verbatim
|
|
*> (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
|
|
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LWORK
|
|
*> \verbatim
|
|
*> The dimension of the array WORK.
|
|
*> LWORK >= MAX( LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) ),
|
|
*> where
|
|
*> NUM_ALL_ROW_BLOCKS = CEIL((M-N)/(MB1-N)),
|
|
*> NB1LOCAL = MIN(NB1,N).
|
|
*> LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL,
|
|
*> LW1 = NB1LOCAL * N,
|
|
*> LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) ),
|
|
*> 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.
|
|
*
|
|
*> \ingroup doubleOTHERcomputational
|
|
*
|
|
*> \par Contributors:
|
|
* ==================
|
|
*>
|
|
*> \verbatim
|
|
*>
|
|
*> November 2020, Igor Kozachenko,
|
|
*> Computer Science Division,
|
|
*> University of California, Berkeley
|
|
*>
|
|
*> \endverbatim
|
|
*>
|
|
* =====================================================================
|
|
SUBROUTINE DGETSQRHRT( M, N, MB1, NB1, NB2, A, LDA, T, LDT, WORK,
|
|
$ LWORK, INFO )
|
|
IMPLICIT NONE
|
|
*
|
|
* -- LAPACK computational routine --
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, LDT, LWORK, M, N, NB1, NB2, MB1
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * ), T( LDT, * ), WORK( * )
|
|
* ..
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
PARAMETER ( ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL LQUERY
|
|
INTEGER I, IINFO, J, LW1, LW2, LWT, LDWT, LWORKOPT,
|
|
$ NB1LOCAL, NB2LOCAL, NUM_ALL_ROW_BLOCKS
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DCOPY, DLATSQR, DORGTSQR_ROW, DORHR_COL,
|
|
$ XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC CEILING, DBLE, MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input arguments
|
|
*
|
|
INFO = 0
|
|
LQUERY = LWORK.EQ.-1
|
|
IF( M.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 .OR. M.LT.N ) THEN
|
|
INFO = -2
|
|
ELSE IF( MB1.LE.N ) THEN
|
|
INFO = -3
|
|
ELSE IF( NB1.LT.1 ) THEN
|
|
INFO = -4
|
|
ELSE IF( NB2.LT.1 ) THEN
|
|
INFO = -5
|
|
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
|
|
INFO = -7
|
|
ELSE IF( LDT.LT.MAX( 1, MIN( NB2, N ) ) ) THEN
|
|
INFO = -9
|
|
ELSE
|
|
*
|
|
* Test the input LWORK for the dimension of the array WORK.
|
|
* This workspace is used to store array:
|
|
* a) Matrix T and WORK for DLATSQR;
|
|
* b) N-by-N upper-triangular factor R_tsqr;
|
|
* c) Matrix T and array WORK for DORGTSQR_ROW;
|
|
* d) Diagonal D for DORHR_COL.
|
|
*
|
|
IF( LWORK.LT.N*N+1 .AND. .NOT.LQUERY ) THEN
|
|
INFO = -11
|
|
ELSE
|
|
*
|
|
* Set block size for column blocks
|
|
*
|
|
NB1LOCAL = MIN( NB1, N )
|
|
*
|
|
NUM_ALL_ROW_BLOCKS = MAX( 1,
|
|
$ CEILING( DBLE( M - N ) / DBLE( MB1 - N ) ) )
|
|
*
|
|
* Length and leading dimension of WORK array to place
|
|
* T array in TSQR.
|
|
*
|
|
LWT = NUM_ALL_ROW_BLOCKS * N * NB1LOCAL
|
|
|
|
LDWT = NB1LOCAL
|
|
*
|
|
* Length of TSQR work array
|
|
*
|
|
LW1 = NB1LOCAL * N
|
|
*
|
|
* Length of DORGTSQR_ROW work array.
|
|
*
|
|
LW2 = NB1LOCAL * MAX( NB1LOCAL, ( N - NB1LOCAL ) )
|
|
*
|
|
LWORKOPT = MAX( LWT + LW1, MAX( LWT+N*N+LW2, LWT+N*N+N ) )
|
|
*
|
|
IF( ( LWORK.LT.MAX( 1, LWORKOPT ) ).AND.(.NOT.LQUERY) ) THEN
|
|
INFO = -11
|
|
END IF
|
|
*
|
|
END IF
|
|
END IF
|
|
*
|
|
* Handle error in the input parameters and return workspace query.
|
|
*
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DGETSQRHRT', -INFO )
|
|
RETURN
|
|
ELSE IF ( LQUERY ) THEN
|
|
WORK( 1 ) = DBLE( LWORKOPT )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( MIN( M, N ).EQ.0 ) THEN
|
|
WORK( 1 ) = DBLE( LWORKOPT )
|
|
RETURN
|
|
END IF
|
|
*
|
|
NB2LOCAL = MIN( NB2, N )
|
|
*
|
|
*
|
|
* (1) Perform TSQR-factorization of the M-by-N matrix A.
|
|
*
|
|
CALL DLATSQR( M, N, MB1, NB1LOCAL, A, LDA, WORK, LDWT,
|
|
$ WORK(LWT+1), LW1, IINFO )
|
|
*
|
|
* (2) Copy the factor R_tsqr stored in the upper-triangular part
|
|
* of A into the square matrix in the work array
|
|
* WORK(LWT+1:LWT+N*N) column-by-column.
|
|
*
|
|
DO J = 1, N
|
|
CALL DCOPY( J, A( 1, J ), 1, WORK( LWT + N*(J-1)+1 ), 1 )
|
|
END DO
|
|
*
|
|
* (3) Generate a M-by-N matrix Q with orthonormal columns from
|
|
* the result stored below the diagonal in the array A in place.
|
|
*
|
|
|
|
CALL DORGTSQR_ROW( M, N, MB1, NB1LOCAL, A, LDA, WORK, LDWT,
|
|
$ WORK( LWT+N*N+1 ), LW2, IINFO )
|
|
*
|
|
* (4) Perform the reconstruction of Householder vectors from
|
|
* the matrix Q (stored in A) in place.
|
|
*
|
|
CALL DORHR_COL( M, N, NB2LOCAL, A, LDA, T, LDT,
|
|
$ WORK( LWT+N*N+1 ), IINFO )
|
|
*
|
|
* (5) Copy the factor R_tsqr stored in the square matrix in the
|
|
* work array WORK(LWT+1:LWT+N*N) into the upper-triangular
|
|
* part of A.
|
|
*
|
|
* (6) Compute from R_tsqr the factor R_hr corresponding to
|
|
* the reconstructed Householder vectors, i.e. R_hr = S * R_tsqr.
|
|
* This multiplication by the sign matrix S on the left means
|
|
* changing the sign of I-th row of the matrix R_tsqr according
|
|
* to sign of the I-th diagonal element DIAG(I) of the matrix S.
|
|
* DIAG is stored in WORK( LWT+N*N+1 ) from the DORHR_COL output.
|
|
*
|
|
* (5) and (6) can be combined in a single loop, so the rows in A
|
|
* are accessed only once.
|
|
*
|
|
DO I = 1, N
|
|
IF( WORK( LWT+N*N+I ).EQ.-ONE ) THEN
|
|
DO J = I, N
|
|
A( I, J ) = -ONE * WORK( LWT+N*(J-1)+I )
|
|
END DO
|
|
ELSE
|
|
CALL DCOPY( N-I+1, WORK(LWT+N*(I-1)+I), N, A( I, I ), LDA )
|
|
END IF
|
|
END DO
|
|
*
|
|
WORK( 1 ) = DBLE( LWORKOPT )
|
|
RETURN
|
|
*
|
|
* End of DGETSQRHRT
|
|
*
|
|
END |