clr_core.F File Reference

Go to the source code of this file.

Modules
module	cmumps_lr_core

Functions/Subroutines
subroutine	cmumps_lr_core::init_lrb (lrb_out, k, m, n, islr)
subroutine	cmumps_lr_core::is_front_blr_candidate (inode, niv, nfront, nass, blron, k489, k490, k491, k492, k20, k60, idad, k38, lrstatus, n, lrgroups)
subroutine	cmumps_lr_core::alloc_lrb (lrb_out, k, m, n, islr, iflag, ierror, keep8)
subroutine	cmumps_lr_core::alloc_lrb_from_acc (acc_lrb, lrb_out, k, m, n, loru, iflag, ierror, keep8)
subroutine	cmumps_lr_core::regrouping2 (cut, npartsass, nass, npartscb, ncb, ibcksz, onlycb, k472)
subroutine	cmumps_lr_core::cmumps_lrtrsm (a, la, poselt_local, nfront, lda, lrb, niv, sym, loru, iw, offset_iw)
subroutine	cmumps_lr_core::cmumps_lrgemm_scaling (lrb, scaled, a, la, diag, ld_diag, iw2, poseltt, nfront, block, maxi_cluster)
subroutine	cmumps_lr_core::cmumps_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	cmumps_lr_core::cmumps_decompress_acc (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, loru, count_flops)
subroutine	cmumps_lr_core::cmumps_compress_fr_updates (acc_lrb, maxi_cluster, maxi_rank, a, la, poseltt, nfront, niv, toleps, tol_opt, kpercent, buildq, loru, cb_compress)
subroutine	cmumps_lr_core::cmumps_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	cmumps_lr_core::cmumps_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	cmumps_lr_core::cmumps_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	cmumps_lr_core::max_cluster (cut, cut_size, maxi_cluster)
subroutine	cmumps_lr_core::cmumps_get_lua_order (nb_blocks, order, rank, iwhandler, sym, fs_or_cb, i, j, frfr_updates, lbandslave_in, k474, blr_u_col)
subroutine	cmumps_lr_core::cmumps_blr_asm_niv1 (a, la, posel1, nfront, nass1, iwhandler, son_iw, liw, lstk, nelim, k1, k2, sym, keep, keep8, opassw)
subroutine	cmumps_truncated_rrqr (m, n, a, lda, jpvt, tau, work, ldw, rwork, toleps, tol_opt, rank, maxrank, info, islr)

Function/Subroutine Documentation

cmumps_truncated_rrqr()

subroutine cmumps_truncated_rrqr	(	integer	m,
		integer	n,
		complex, dimension(lda,*)	a,
		integer	lda,
		integer, dimension(*)	jpvt,
		complex, dimension(*)	tau,
		complex, dimension(ldw,*)	work,
		integer	ldw,
		real, dimension(*)	rwork,
		real	toleps,
		integer	tol_opt,
		integer	rank,
		integer	maxrank,
		integer	info,
		logical	islr )

Definition at line 1608 of file clr_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
      REAL               ::  TOLEPS
      INTEGER            ::  JPVT(*)
      REAL               ::  RWORK(*)
      COMPLEX            ::  A(LDA,*), TAU(*)
      COMPLEX            ::  WORK(LDW,*)
      LOGICAL            ::  ISLR
      REAL               ::  TOLEPS_EFF, TRUNC_ERR
      INTEGER, PARAMETER ::  INB=1, inbmin=2
      INTEGER            :: J, JB, MINMN, NB
      INTEGER            :: OFFSET, ITEMP
      INTEGER            :: LSTICC, PVT, K, RK
      REAL               :: TEMP, TEMP2, TOL3Z
      COMPLEX            :: AKK
      LOGICAL INADMISSIBLE
      REAL, PARAMETER    :: RZERO=0.0e+0, rone=1.0e+0
      COMPLEX :: ZERO
      COMPLEX :: ONE
      parameter( one = ( 1.0e+0, 0.0e+0 ) )
      parameter( zero = ( 0.0e+0, 0.0e+0 ) )
      REAL               :: slamch
      INTEGER            :: ilaenv, isamax
      EXTERNAL           :: isamax, slamch
      EXTERNAL           cgeqrf, cunmqr, xerbla
      EXTERNAL           ilaenv
      EXTERNAL           cgemm, cgemv, clarfg, cswap
      REAL, EXTERNAL :: scnrm2
      REAL, EXTERNAL :: snrm2
      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 CMUMPS_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 ) = scnrm2( M, A( 1, J ), 1 )
         rwork( j ) = scnrm2( 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 = snrm2( N, VN1( 1 ), 1 )
        trunc_err = snrm2( n, rwork( 1 ), 1 )
      ENDIF
      offset = 0
      tol3z  = sqrt(slamch('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 ) + ISAMAX( N-RK+1, VN1( RK ), 1 )
            pvt = ( rk-1 ) + isamax( 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 cswap( m, a( 1, pvt ), 1, a( 1, rk ), 1 )
c              CALL cswap( K-1, F( PVT-OFFSET, 1 ), LDF,
c    &              F( K, 1 ), LDF )
               CALL cswap( 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
               DO j = 1, k-1
C                 F( K, J ) = CONJG( F( K, J ) )
                  work( k, j+1 ) = conjg( work( k, j+1 ) )
               END DO
               CALL cgemv( '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 )
               DO j = 1, k - 1
C                 F( K, J ) = CONJG( F( K, J ) )
                  work( k, j + 1 ) = conjg( work( k, j + 1 ) )
               END DO
            END IF
*     Generate elementary reflector H(k).
            IF( rk.LT.m ) THEN
               CALL clarfg( m-rk+1, a(rk,rk), a(rk+1,rk), 1, tau(rk) )
            ELSE
               CALL clarfg( 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 cgemv( 'Conjugate 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 cgemv( 'Conjugate 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 cgemv( '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.
.LT.            IF( RKN ) THEN
               CALL cgemm( 'no transpose', 'conjugate transpose',
     &              1, N-RK,
C    &              K, -ONE, A( RK, OFFSET+1 ), LDA, F( K+1, 1 ), LDF,
     &              K, -ONE, A( RK, OFFSET+1 ), LDA, WORK( K+1,2 ), LDW,
     &              ONE, A( RK, RK+1 ), LDA )
            END IF
*     Update partial column norms.
*     
.LT.            IF( RKMINMN ) THEN
               DO J = RK + 1, N
C                 IF( VN1( J ).NE.RZERO ) THEN
.NE.                  IF( RWORK( J )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
.LE.                     IF( TEMP2  TOL3Z ) THEN
C                       VN2( J ) = REAL( LSTICC )
                        RWORK( N+J ) = REAL( 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
.NE.            IF (LSTICC0) EXIT
.LT.            IF (TOL_OPT0) THEN
C             Compute TRUNC_ERR for next step              
C             TRUNC_ERR = snrm2( N-RK, VN1( RK+1 ), 1 )
              TRUNC_ERR = snrm2( 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.
.LT.         IF( RKMIN(N,M) ) THEN
            CALL cgemm( 'no transpose', 'conjugate 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.
.GT.         DO WHILE( LSTICC0 ) 
C           ITEMP = NINT( VN2( LSTICC ) )
            ITEMP = NINT( RWORK( N + LSTICC ) )
C           VN1( LSTICC ) = scnrm2( M-RK, A( RK+1, LSTICC ), 1 )
            RWORK( LSTICC ) = scnrm2( 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
.GE.         IF(RKMINMN) EXIT
         OFFSET = RK
.LT.         IF (TOL_OPT0) THEN
C          Compute TRUNC_ERR for next step              
C          TRUNC_ERR = snrm2( N-RK, VN1( RK+1 ), 1 )
           TRUNC_ERR = snrm2( N-RK, RWORK( RK+1 ), 1 )
         ENDIF
      END DO
      RANK = RK
.NOT..GT.        ISLR = (RKMAXRANK)
      RETURN
 999  FORMAT ('on entry to cmumps_truncated_rrqr, parameter number',
     &            I2,' had an illegal value')