dgetsqrhrt()

subroutine dgetsqrhrt	(	integer	m,
		integer	n,
		integer	mb1,
		integer	nb1,
		integer	nb2,
		double precision, dimension( lda, * )	a,
		integer	lda,
		double precision, dimension( ldt, * )	t,
		integer	ldt,
		double precision, dimension( * )	work,
		integer	lwork,
		integer	info
	)

DGETSQRHRT

Download DGETSQRHRT + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 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.

Parameters

[in]	M	M is INTEGER The number of rows of the matrix A. M >= 0.
[in]	N	N is INTEGER The number of columns of the matrix A. M >= N >= 0.
[in]	MB1	MB1 is INTEGER The row block size to be used in the blocked TSQR. MB1 > N.
[in]	NB1	NB1 is INTEGER The column block size to be used in the blocked TSQR. N >= NB1 >= 1.
[in]	NB2	NB2 is INTEGER The block size to be used in the blocked QR that is output. NB2 >= 1.
[in,out]	A	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).
[in]	LDA	LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).
[out]	T	T is DOUBLE PRECISION array, dimension (LDT,N)) The upper triangular block reflectors stored in compact form as a sequence of upper triangular blocks.
[in]	LDT	LDT is INTEGER The leading dimension of the array T. LDT >= NB2.
[out]	WORK	(workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]	LWORK	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.
[out]	INFO	INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

 November 2020, Igor Kozachenko,
                Computer Science Division,
                University of California, Berkeley

Definition at line 177 of file dgetsqrhrt.f.

      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
*

Here is the call graph for this function:

Here is the caller graph for this function: