dlr_core.F File Reference

Go to the source code of this file.

Modules
module	dmumps_lr_core

Functions/Subroutines
subroutine	dmumps_lr_core::init_lrb (lrb_out, k, m, n, islr)
subroutine	dmumps_lr_core::is_front_blr_candidate (inode, niv, nfront, nass, blron, k489, k490, k491, k492, k20, k60, idad, k38, lrstatus, n, lrgroups)
subroutine	dmumps_lr_core::alloc_lrb (lrb_out, k, m, n, islr, iflag, ierror, keep8)
subroutine	dmumps_lr_core::alloc_lrb_from_acc (acc_lrb, lrb_out, k, m, n, loru, iflag, ierror, keep8)
subroutine	dmumps_lr_core::regrouping2 (cut, npartsass, nass, npartscb, ncb, ibcksz, onlycb, k472)
subroutine	dmumps_lr_core::dmumps_lrtrsm (a, la, poselt_local, nfront, lda, lrb, niv, sym, loru, iw, offset_iw)
subroutine	dmumps_lr_core::dmumps_lrgemm_scaling (lrb, scaled, a, la, diag, ld_diag, iw2, poseltt, nfront, block, maxi_cluster)
subroutine	dmumps_lr_core::dmumps_lrgemm4 (alpha, lrb1, lrb2, beta, a, la, poseltt, nfront, sym, iflag, ierror, midblk_compress, toleps, tol_opt, kpercent, rank, buildq, lua_activated, loru, lrb3, maxi_rank, maxi_cluster, diag, ld_diag, iw2, block)
subroutine	dmumps_lr_core::dmumps_decompress_acc (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, loru, count_flops)
subroutine	dmumps_lr_core::dmumps_compress_fr_updates (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, toleps, tol_opt, kpercent, buildq, loru, cb_compress)
subroutine	dmumps_lr_core::dmumps_recompress_acc (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, midblk_compress, toleps, tol_opt, kpercent_rmb, kpercent_lua, new_acc_rank)
recursive subroutine	dmumps_lr_core::dmumps_recompress_acc_narytree (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, keep8, nfront, niv, midblk_compress, toleps, tol_opt, kpercent_rmb, kpercent_lua, k478, rank_list, pos_list, nb_nodes, level, acc_tmp)
subroutine	dmumps_lr_core::dmumps_recompress_acc_v2 (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, midblk_compress, toleps, tol_opt, kpercent_rmb, kpercent_lua, new_acc_rank)
subroutine	dmumps_lr_core::max_cluster (cut, cut_size, maxi_cluster)
subroutine	dmumps_lr_core::dmumps_get_lua_order (nb_blocks, order, rank, iwhandler, sym, fs_or_cb, i, j, frfr_updates, lbandslave_in, k474, blr_u_col)
subroutine	dmumps_lr_core::dmumps_blr_asm_niv1 (a, la, posel1, nfront, nass1, iwhandler, son_iw, liw, lstk, nelim, k1, k2, sym, keep, keep8, opassw)
subroutine	dmumps_truncated_rrqr (m, n, a, lda, jpvt, tau, work, ldw, rwork, toleps, tol_opt, rank, maxrank, info, islr)

Function/Subroutine Documentation

dmumps_truncated_rrqr()

subroutine dmumps_truncated_rrqr	(	integer	m,
		integer	n,
		double precision, dimension(lda,*)	a,
		integer	lda,
		integer, dimension(*)	jpvt,
		double precision, dimension(*)	tau,
		double precision, dimension(ldw,*)	work,
		integer	ldw,
		double precision, dimension(*)	rwork,
		double precision	toleps,
		integer	tol_opt,
		integer	rank,
		integer	maxrank,
		integer	info,
		logical	islr )

Definition at line 1608 of file dlr_core.F.

C     This routine computes a Rank-Revealing QR factorization of a dense
C     matrix A. The factorization is truncated when the absolute value of
C     a diagonal coefficient of the R factor becomes smaller than a
C     prescribed threshold TOLEPS. The resulting partial Q and R factors
C     provide a rank-k approximation of the input matrix A with accuracy
C     TOLEPS.
C     
C     This routine is obtained by merging the LAPACK
C     (http://www.netlib.org/lapack/) CGEQP3 and CLAQPS routines and by
C     applying a minor modification to the outer factorization loop in
C     order to stop computations as soon as possible when the required
C     accuracy is reached.
C
C     Copyright (c) 1992-2017 The University of Tennessee and The 
C     University of Tennessee Research Foundation.  All rights reserved.
C     Copyright (c) 2000-2017 The University of California Berkeley. 
C     All rights reserved.
C     Copyright (c) 2006-2017 The University of Colorado Denver.  
C     All rights reserved.
C
C     Redistribution and use in source and binary forms, with or without
C     modification, are permitted provided that the following conditions
C     are met:
C
C      - Redistributions of source code must retain the above copyright
C        notice, this list of conditions and the following disclaimer.
C
C      - Redistributions in binary form must reproduce the above 
C        copyright notice, this list of conditions and the following 
C        disclaimer listed in this license in the documentation and/or 
C        other materials provided with the distribution.
C
C      - Neither the name of the copyright holders nor the names of its
C        contributors may be used to endorse or promote products derived from
C        this software without specific prior written permission.
C
C      The copyright holders provide no reassurances that the source code
C      provided does not infringe any patent, copyright, or any other
C      intellectual property rights of third parties.  The copyright holders
C      disclaim any liability to any recipient for claims brought against
C      recipient by any third party for infringement of that parties
C      intellectual property rights.
C
C      THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C      "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C      LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C      A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C      OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C      SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C      LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C      DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C      THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C      (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C      OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C
      IMPLICIT NONE
C
      INTEGER            ::  INFO, LDA, LDW, M, N, RANK, MAXRANK
C     TOL_OPT controls the tolerance option used      
C       >0 => use 2-norm (||.||_X = ||.||_2)
C       <0 => use Frobenius-norm (||.||_X = ||.||_F)
C     Furthermore, depending on abs(TOL_OPT):      
C       1 => absolute: ||B_{I(k+1:end),J(k+1:end)}||_X <= TOLEPS     
C       2 => relative to 2-norm of the compressed block: 
C        ||B_{I(k+1:end),J(k+1:end)}||_X <= TOLEPS*||B_{I,J}||_2
C       3 => relative to the max of the 2-norms of the row and column diagonal blocks 
C        ||B_{I(k+1:end),J{k+1:end}}||_X <= TOLEPS*max(||B_{I,I}||_2,||B_{J,J}||_2)
C       4 => relative to the sqrt of product of the 2-norms of the row and column diagonal blocks 
C        ||B_{I(k+1:end),J{k+1:end}}||_X <= TOLEPS*sqrt(||B_{I,I}||_2*||B_{J,J}||_2)
      INTEGER            ::  TOL_OPT
      DOUBLE PRECISION               ::  TOLEPS
      INTEGER            ::  JPVT(*)
      DOUBLE PRECISION               ::  RWORK(*)
      DOUBLE PRECISION            ::  A(LDA,*), TAU(*)
      DOUBLE PRECISION            ::  WORK(LDW,*)
      LOGICAL            ::  ISLR
      DOUBLE PRECISION               ::  TOLEPS_EFF, TRUNC_ERR
      INTEGER, PARAMETER ::  INB=1, inbmin=2
      INTEGER            :: J, JB, MINMN, NB
      INTEGER            :: OFFSET, ITEMP
      INTEGER            :: LSTICC, PVT, K, RK
      DOUBLE PRECISION               :: TEMP, TEMP2, TOL3Z
      DOUBLE PRECISION            :: AKK
      LOGICAL INADMISSIBLE
      DOUBLE PRECISION, PARAMETER    :: RZERO=0.0d+0, rone=1.0d+0
      DOUBLE PRECISION :: ZERO
      DOUBLE PRECISION :: ONE
      parameter( one = 1.0d+0 )
      parameter( zero = 0.0d+0 ) 
      DOUBLE PRECISION               :: dlamch
      INTEGER            :: ilaenv, idamax
      EXTERNAL           :: idamax, dlamch
      EXTERNAL           dgeqrf, dormqr, xerbla
      EXTERNAL           ilaenv
      EXTERNAL           dgemm, dgemv, dlarfg, dswap
      DOUBLE PRECISION, EXTERNAL :: dnrm2
      info = 0
      islr = .false.
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 ) THEN
         info = -2
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -4
      END IF
      IF( info.EQ.0 ) THEN
         IF( ldw.LT.n ) THEN
            info = -8
         END IF
      END IF
      IF( info.NE.0 ) THEN
         WRITE(*,999) -info
         RETURN
      END IF
      minmn = min(m,n)
      IF( minmn.EQ.0 ) THEN
         rank = 0
         RETURN
      END IF
      nb = ilaenv( inb, 'CGEQRF', ' ', m, n, -1, -1 )
      SELECT CASE(abs(tol_opt))
      CASE(1)
        toleps_eff = toleps
      CASE(2)
C      TOLEPS_EFF will be computed at step K=1 below        
      CASE DEFAULT
        write(*,*) 'Internal error in DMUMPS_TRUNCATED_RRQR: TOL_OPT =',
     &        tol_opt
        CALL mumps_abort()
      END SELECT
      toleps_eff = toleps
C
C     Avoid pointers (and TARGET attribute on RWORK/WORK)
C     because of implicit interface. An implicit interface
C     is needed to avoid intermediate array copies
C     VN1  => RWORK(1:N)
C     VN2  => RWORK(N+1:2*N)
C     AUXV => WORK(1:LDW,1:1)
C     F    => WORK(1:LDW,2:NB+1)
C     LDF  =  LDW
*     Initialize partial column norms. The first N elements of work
*     store the exact column norms.
      DO j = 1, n
C        VN1( J ) = dnrm2( M, A( 1, J ), 1 )
         rwork( j ) = dnrm2( m, a( 1, j ), 1 )
C        VN2( J ) = VN1( J )
         rwork( n + j ) = rwork( j )
         jpvt(j) = j
      END DO
      IF (tol_opt.LT.0) THEN
C       Compute TRUNC_ERR for first step              
C       TRUNC_ERR = dnrm2( N, VN1( 1 ), 1 )
        trunc_err = dnrm2( n, rwork( 1 ), 1 )
      ENDIF
      offset = 0
      tol3z  = sqrt(dlamch('Epsilon'))
      DO 
         jb     = min(nb,minmn-offset)
         lsticc = 0
         k      = 0
         DO 
            IF(k.EQ.jb) EXIT
            k   = k+1
            rk  = offset+k
C           PVT = ( RK-1 ) + IDAMAX( N-RK+1, VN1( RK ), 1 )
            pvt = ( rk-1 ) + idamax( n-rk+1, rwork( rk ), 1 )
            IF (rk.EQ.1) THEN 
C             IF (abs(TOL_OPT).EQ.2) TOLEPS_EFF = VN1(PVT)*TOLEPS
              IF (abs(tol_opt).EQ.2) toleps_eff = rwork(pvt)*toleps
            ENDIF
            IF (tol_opt.GT.0) THEN
C             TRUNC_ERR = VN1(PVT)
              trunc_err = rwork(pvt)
C           ELSE
C             TRUNC_ERR has been already computed at previous step
            ENDIF
            IF(trunc_err.LT.toleps_eff) THEN
              rank = rk-1
              islr = .true.  
              RETURN
            ENDIF
              inadmissible = (rk.GT.maxrank)
            IF (inadmissible) THEN
               rank = rk
               info = rk
               islr = .false.
               RETURN
            END IF
            IF( pvt.NE.rk ) THEN
               CALL dswap( m, a( 1, pvt ), 1, a( 1, rk ), 1 )
c              CALL dswap( K-1, F( PVT-OFFSET, 1 ), LDF,
c    &              F( K, 1 ), LDF )
               CALL dswap( k-1, work( pvt-offset, 2 ), ldw,
     &              work( k, 2 ), ldw )
               itemp     = jpvt(pvt)
               jpvt(pvt) = jpvt(rk)
               jpvt(rk)  = itemp
C              VN1(PVT)  = VN1(RK)
C              VN2(PVT)  = VN2(RK)
               rwork(pvt)    = rwork(rk)
               rwork(n+pvt)  = rwork(n+rk)
            END IF
*     Apply previous Householder reflectors to column K:
*     A(RK:M,RK) := A(RK:M,RK) - A(RK:M,OFFSET+1:RK-1)*F(K,1:K-1)**H.
            IF( k.GT.1 ) THEN
               CALL dgemv( 'No transpose', m-rk+1, k-1, -one,
C    &              A(RK,OFFSET+1), LDA, F(K,1), LDF,
     &              a(rk,offset+1), lda, work(k,2), ldw,
     &              one, a(rk,rk), 1 )
            END IF
*     Generate elementary reflector H(k).
            IF( rk.LT.m ) THEN
               CALL dlarfg( m-rk+1, a(rk,rk), a(rk+1,rk), 1, tau(rk) )
            ELSE
               CALL dlarfg( 1, a(rk,rk), a(rk,rk), 1, tau(rk) )
            END IF
            akk      = a(rk,rk)
            a(rk,rk) = one
*     Compute Kth column of F:
*     F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)**H*A(RK:M,K).
            IF( rk.LT.n ) THEN
               CALL dgemv( 'Transpose', m-rk+1, n-rk, tau(rk),
     &              a(rk,rk+1), lda, a(rk,rk), 1, zero,
C    &              F( K+1, K ), 1 )
     &              work( k+1, k+1 ), 1 )
            END IF
*     Padding F(1:K,K) with zeros.
            DO j = 1, k
C              F( J, K ) = ZERO
               work( j, k+1 ) = zero
            END DO
*     Incremental updating of F:
*     F(1:N,K) := F(1:N-OFFSET,K) - 
*             tau(RK)*F(1:N,1:K-1)*A(RK:M,OFFSET+1:RK-1)**H*A(RK:M,RK).
            IF( k.GT.1 ) THEN
               CALL dgemv( 'Transpose', m-rk+1, k-1, -tau(rk),
     &              a(rk,offset+1), lda, a(rk,rk), 1, zero,
     &              work(1,1), 1 )
C    &              AUXV(1,1), 1 )
               CALL dgemv( 'No transpose', n-offset, k-1, one,
     &              work(1,2), ldw, work(1,1), 1, one, work(1,k+1), 1 )
C    &              F(1,1), LDF, AUXV(1,1), 1, ONE, F(1,K), 1 )
            END IF
*     Update the current row of A:
*     A(RK,RK+1:N) := A(RK,RK+1:N) - A(RK,OFFSET+1:RK)*F(K+1:N,1:K)**H.
            IF( rk.LT.n ) THEN
C              CALL dgemv( 'No Transpose', N-RK, K, -ONE, F( K+1, 1 ), 
               CALL dgemv( 'No Transpose', n-rk, k, -one, work( k+1,2 ),
     &              ldw,
     &              a( rk, offset+1 ), lda, one, a( rk, rk+1 ), lda )
            END IF
*     Update partial column norms.
*     
            IF( rk.LT.minmn ) THEN
               DO j = rk + 1, n
C                 IF( VN1( J ).NE.RZERO ) THEN
                  IF( rwork( j ).NE.rzero ) THEN
*     
*     NOTE: The following 4 lines follow from the analysis in
*     Lapack Working Note 176.
*
C                    TEMP = ABS( A( RK, J ) ) / VN1( J )
                     temp = abs( a( rk, j ) ) / rwork( j )
                     temp = max( rzero, ( rone+temp )*( rone-temp ) )
C                    TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
                     temp2 = temp*( rwork( j ) / rwork( n+j ) )**2
                     IF( temp2 .LE. tol3z ) THEN
C                       VN2( J ) = dble( LSTICC )
                        rwork( n+j ) = dble( lsticc )
                        lsticc = j
                     ELSE
C                       VN1( J ) = VN1( J )*SQRT( TEMP )
                        rwork( j ) = rwork( j )*sqrt( temp )
                     END IF
                  END IF
               END DO
            END IF
            a( rk, rk ) = akk
            IF (lsticc.NE.0) EXIT
            IF (tol_opt.LT.0) THEN
C             Compute TRUNC_ERR for next step              
C             TRUNC_ERR = dnrm2( N-RK, VN1( RK+1 ), 1 )
              trunc_err = dnrm2( n-rk, rwork( rk+1 ), 1 )
            ENDIF
         END DO
*     Apply the block reflector to the rest of the matrix:
*     A(RK+1:M,RK+1:N) := A(RK+1:M,RK+1:N) -
*     A(RK+1:M,OFFSET+1:RK)*F(K+1:N-OFFSET,1:K)**H.
         IF( rk.LT.min(n,m) ) THEN
            CALL dgemm( 'No transpose', 'Transpose', m-rk,
     &           n-rk, k, -one, a(rk+1,offset+1), lda,
C    &           F(K+1,1), LDF, ONE, A(RK+1,RK+1), LDA )
     &           work(k+1,2), ldw, one, a(rk+1,rk+1), lda )
         END IF
*     Recomputation of difficult columns.
         DO WHILE( lsticc.GT.0 ) 
C           ITEMP = NINT( VN2( LSTICC ) )
            itemp = nint( rwork( n + lsticc ) )
C           VN1( LSTICC ) = dnrm2( M-RK, A( RK+1, LSTICC ), 1 )
            rwork( lsticc ) = dnrm2( m-rk, a( rk+1, lsticc ), 1 )
*     
*     NOTE: The computation of RWORK( LSTICC ) relies on the fact that 
*     SNRM2 does not fail on vectors with norm below the value of
*     SQRT(DLAMCH('S')) 
*     
C           VN2( LSTICC ) = VN1( LSTICC )
            rwork( n + lsticc ) = rwork( lsticc )
            lsticc = itemp
         END DO
         IF(rk.GE.minmn) EXIT
         offset = rk
         IF (tol_opt.LT.0) THEN
C          Compute TRUNC_ERR for next step              
C          TRUNC_ERR = dnrm2( N-RK, VN1( RK+1 ), 1 )
           trunc_err = dnrm2( n-rk, rwork( rk+1 ), 1 )
         ENDIF
      END DO
      rank = rk
        islr = .NOT.(rk.GT.maxrank)
      RETURN
 999  FORMAT ('On entry to DMUMPS_TRUNCATED_RRQR, parameter number',
     &            i2,' had an illegal value')