slr_core.F File Reference

Go to the source code of this file.

Modules
module	smumps_lr_core

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

Function/Subroutine Documentation

smumps_truncated_rrqr()

subroutine smumps_truncated_rrqr	(	integer	m,
		integer	n,
		real, dimension(lda,*)	a,
		integer	lda,
		integer, dimension(*)	jpvt,
		real, dimension(*)	tau,
		real, 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 slr_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(*)
      REAL            ::  A(LDA,*), TAU(*)
      REAL            ::  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
      REAL            :: AKK
      LOGICAL INADMISSIBLE
      REAL, PARAMETER    :: RZERO=0.0e+0, rone=1.0e+0
      REAL :: ZERO
      REAL :: ONE
      parameter( one = 1.0e+0 )
      parameter( zero = 0.0e+0 ) 
      REAL               :: slamch
      INTEGER            :: ilaenv, isamax
      EXTERNAL           :: isamax, slamch
      EXTERNAL           sgeqrf, sormqr, xerbla
      EXTERNAL           ilaenv
      EXTERNAL           sgemm, sgemv, slarfg, sswap
      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 smumps_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 ) = snrm2( M, A( 1, J ), 1 )
         RWORK( J ) = snrm2( M, A( 1, J ), 1 )
C        VN2( J ) = VN1( J )
         RWORK( N + J ) = RWORK( J )
         JPVT(J) = J
      END DO
.LT.      IF (TOL_OPT0) 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 
.EQ.            IF(KJB) 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 )
.EQ.            IF (RK1) THEN 
C             IF (abs(TOL_OPT).EQ.2) TOLEPS_EFF = VN1(PVT)*TOLEPS
.EQ.              IF (abs(TOL_OPT)2) TOLEPS_EFF = RWORK(PVT)*TOLEPS
            ENDIF
.GT.            IF (TOL_OPT0) THEN
C             TRUNC_ERR = VN1(PVT)
              TRUNC_ERR = RWORK(PVT)
C           ELSE
C             TRUNC_ERR has been already computed at previous step
            ENDIF
.LT.            IF(TRUNC_ERRTOLEPS_EFF) THEN
              RANK = RK-1
              ISLR = .TRUE.  
              RETURN
            ENDIF
.GT.              INADMISSIBLE = (RKMAXRANK)
            IF (INADMISSIBLE) THEN
               RANK = RK
               INFO = RK
               ISLR = .FALSE.
               RETURN
            END IF
.NE.            IF( PVTRK ) THEN
               CALL sswap( M, A( 1, PVT ), 1, A( 1, RK ), 1 )
c              CALL sswap( K-1, F( PVT-OFFSET, 1 ), LDF,
c    &              F( K, 1 ), LDF )
               CALL sswap( 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.
.GT.            IF( K1 ) THEN
               CALL sgemv( '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).
.LT.            IF( RKM ) THEN
               CALL slarfg( M-RK+1, A(RK,RK), A(RK+1,RK), 1, TAU(RK) )
            ELSE
               CALL slarfg( 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).
.LT.            IF( RKN ) THEN
               CALL sgemv( '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).
.GT.            IF( K1 ) THEN
               CALL sgemv( '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 sgemv( '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
C              CALL sgemv( 'No Transpose', N-RK, K, -ONE, F( K+1, 1 ), 
               CALL sgemv( '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.
*     
.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 sgemm( '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.
.GT.         DO WHILE( LSTICC0 ) 
C           ITEMP = NINT( VN2( LSTICC ) )
            ITEMP = NINT( RWORK( N + LSTICC ) )
C           VN1( LSTICC ) = snrm2( M-RK, A( RK+1, LSTICC ), 1 )
            RWORK( LSTICC ) = snrm2( 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 smumps_truncated_rrqr, parameter number',
     &            I2,' had an illegal value')