getPartialInv.F90

Go to the documentation of this file.

!======================================================================
!     Module to calculate the partial inverse of a real/complex sparse matrix via the PARDISO libraries.
!     This file is part of the EQuUs package. (for details see: http://eqt.elte.hu/equus/home)
!     Copyright (C) 2016 Peter Rakyta, Ph.D.
!
!     This program is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
! 
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program.  If not, see http://www.gnu.org/licenses/.
!     
!======================================================================
! MODULE getPartialInv
MODULE getpartialinv

#ifndef MAXNRHS
#define MAXNRHS 100
#endif

#ifndef MAXMEM
#define MAXMEM 16000000 !In kilobytes
#endif

!****************************************************************      
  CONTAINS
!****************************************************************

! DESCRIPTION:
#ifdef MIC
!dir$ attributes offload:mic :: dgetPartialInv
#endif
  subroutine dgetpartialinv(sizeA, nonzerosA, a, ia, ja, sizeInv, invA, error)
    use mkl_pardiso
    use mkl_service

    IMPLICIT NONE


!     Internal solver memory pointer 
    TYPE(mkl_pardiso_handle), pointer :: PT(:) => null()
    integer*4, pointer :: idum(:) => null(), pia(:) => null(), pja(:) => null()
    real(KIND=8), pointer :: ddum(:) => null(), pa(:) => null(), pb(:,:) => null(), px(:,:) => null()
    INTEGER*4 maxfct, mnum, mtype, phase, nrhs, error, msglvl, i, nnz, numIterations, idxIterations, offset
    INTEGER*4 iparm(64)


!     The matrix is given in the CSR format A = (a, ia, ja)
      INTEGER*4 sizeA, nonzerosA !matrix A is of size sizeA x sizeA
      INTEGER*4, TARGET ::  ia(sizeA+1), ja(nonzerosA) 
      real*8, TARGET :: a(nonzerosa)

      INTEGER*4, ALLOCATABLE :: perm(:)

      INTEGER*4 sizeInv !size of the inverse matrix
      real*8 :: inva(sizeinv,sizeinv) 
      
      INTEGER*4 idx, jdx
      integer maxmem

      DATA nrhs /1/, maxfct /1/, mnum /1/
      allocate( pt(64), idum(0), ddum(0) )

      do i = 1, 64
        pt(i)%DUMMY = 0
      end do
      pa => a
      pia => ia
      pja => ja
! initialize error flag
      error = 0

#ifdef DEBUG
      write(0,*)
      write(0,*) 'The input matrix: (getPartialInv)'
      idx = 0
      do jdx = 1, nonzerosa
        if (jdx >= ia(idx+1)) then
          idx=idx+1
        end if

        if (jdx < ia(idx+1)) then  
          !write(0,*) 'row', idx, 'col', ja(jdx), 'element of a = ', a(jdx)
        end if
      end do
      write(0,*) 'The size of the Partial inverse to be calculated:', sizeinv
#endif

!
!  .. Setup Pardiso control parameters und initialize the solvers     
!     internal adress pointers. This is only necessary for the FIRST   
!     call of the PARDISO solver.                                     
!       
!     choose the matrix type for pardiso factorization
      mtype     = 1  ! structurally symmetric real matrix

!  .. PARDISO license check and initialize solver
      call pardisoinit(pt, mtype, iparm)
      




!     Reordering and Factorization, This step also allocates
!     all memory that is necessary for the factorization
#ifdef DEBUG      
      msglvl    = 1      ! with statistical information
      iparm(27) = 1      ! perform matrix check
#else
      msglvl    = 0      ! without statistical information
      iparm(27) = 0      ! perform matrix check
#endif
      iparm(6) = 0  ! The solution is written into x  
      iparm(28) = 0 ! double real/complex

#ifdef MIC
      iparm(2) = 3 ! parallel OpenMP version of the nested dissection
      iparm(24) = 1 ! Intel MKL PARDISO uses a two-level factorization algorithm
      iparm(25) = 0 ! Intel MKL PARDISO uses a parallel algorithm for the solve step.
      iparm(34) = 0 ! 
      iparm(60) = 0 ! in-core (IC) mode
#else
      iparm(2) = 2 ! 
      iparm(24) = 0 ! Intel MKL PARDISO uses a single-level factorization algorithm
      iparm(25) = 0 ! Intel MKL PARDISO uses a parallel algorithm for the solve step.
      iparm(34) = 0 ! 
      iparm(60) = 0 ! in-core (IC) mode
#endif

      phase     = 11  ! analysis

      ! setting the appropriate number of right hand sides.
#ifdef MIC
      nrhs = maxnrhs      
      if (nrhs.GT.sizeinv) then
        nrhs = sizeinv
      end if
#else
      nrhs = 1
#endif

#ifdef DEBUG
      write(0,*) 'number of right sides: ', nrhs
#endif 
      
      CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, pa, pia, pja, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error)
      maxmem = (iparm(16)+iparm(17))

#ifdef DEBUG
      write(0,*) 'maxmem: ', maxmem, ' kbytes'
      write(0,*) 'maxmem total: ', maxmem+sizea*nrhs*16/1000, ' kbytes'
#endif

      IF ((maxmem+sizea*nrhs*16/1000).GT.maxmem) then
        nrhs = floor(real(maxmem-maxmem)/real(sizea*16)*1000)
        write(0,*) 'number of right sides: ', nrhs
        IF (nrhs.EQ.0) THEN
          error = -2
          RETURN
        ENDIF
      ENDIF

      phase     = 22  ! factorization 
      CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, pa, pia, pja, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error)

#ifdef DEBUG
      write(0,*) 'Reordering and factorization completed ... '
#endif
      IF (error .NE. 0) THEN
         write(0,*) 'The following ERROR was detected at the factorization step: ', error
        RETURN
      ENDIF 



!     Set parameters to calculate only the last sizeInv components of the solution
      iparm(31) = 3
      phase = 33 ! solver
      allocate( perm(sizea) )
      perm = 0
      do idx = 1, sizeinv
        perm(sizea-sizeinv+idx) = 1
      end do

!     Calculating the chunk of the solution to be calculated at once
      if (modulo(sizeinv, nrhs).EQ.0) then
        numiterations = sizeinv/nrhs
      else
        numiterations = floor(real(sizeinv)/real(nrhs)) + 1
      end if

#ifdef DEBUG
      write(0,*) 'Number of iterations to obtain the right solutions ', numiterations
#endif




!     Allocating the right hand side
      allocate( pb(sizea, nrhs) )
      allocate( px(sizea, nrhs) )

      do idxiterations = 1, numiterations

        offset = (idxiterations-1)*nrhs

        if ((sizeinv-offset).LT.nrhs) then
          nrhs = sizeinv-offset
          deallocate(pb, px) 
          allocate( pb(sizea, nrhs) )
          allocate( px(sizea, nrhs) )
        end if

        if ((sizeinv-offset).LT.nrhs) then
          nrhs = sizeinv-offset
        end if

!       Creating the right hand side
        pb = 0.d0
        px = 0.d0
        do idx = 1, nrhs
          pb(sizea-sizeinv+offset+idx, idx) = 1.d0
        end do
!       Solver        
        CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, pa, pia, pja, &
                   perm, nrhs, iparm, msglvl, pb, px, error)    

        IF (error .NE. 0) THEN
          exit
        ENDIF 


        do idx = 1, nrhs 
          do jdx = 1,sizeinv
            inva(offset+idx,jdx) = px(sizea-sizeinv+jdx, idx)
          end do
        end do  

      end do

      deallocate(pb, px, perm)


      IF (error .NE. 0) THEN
        write(0,*) 'The following ERROR was detected at the solution step: ', error
        RETURN
      ENDIF 

 
#ifdef DEBUG
      do idx = 1,sizeinv
        do jdx = 1,sizeinv
          write(0,*) 'row', idx, 'col', jdx, 'element of a = ', inva(jdx, idx)
        end do
      end do
#endif

      

!  
!   Termination and release of memory
      phase     = -1           ! release internal memory
      CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, ddum, &
                   idum, idum, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error) 

      deallocate(idum, ddum, pt)

      IF (error .NE. 0) THEN
        write(0,*) 'The following ERROR was detected at releasing the memory: ', error
        RETURN
      ENDIF 


  end



! DESCRIPTION:
#ifdef MIC
!dir$ attributes offload:mic :: zgetPartialInv
#endif
  subroutine zgetpartialinv(sizeA, nonzerosA, a, ia, ja, sizeInv, real_invA, imag_invA, error)

    use mkl_pardiso

    IMPLICIT NONE

!     Internal solver memory pointer 
    TYPE(mkl_pardiso_handle), pointer :: PT(:) => null()
    integer*4, pointer :: idum(:) => null(), pia(:) => null(), pja(:) => null()
    complex(KIND=8), pointer :: ddum(:) => null(), pa(:) => null(), pb(:,:) => null(), px(:,:) => null()
    INTEGER*4 maxfct, mnum, mtype, phase, nrhs, error, msglvl, nnz, numIterations, idxIterations, offset
    INTEGER*4 iparm(64)


!     The matrix is given in the CSR format A = (a, ia, ja)
      INTEGER*4 sizeA, nonzerosA !matrix A is of size sizeA x sizeA
      INTEGER*4, TARGET ::  ia(sizeA+1), ja(nonzerosA) 
      COMPLEX*16, TARGET :: a(nonzerosA)

      INTEGER*4, ALLOCATABLE :: perm(:)

      INTEGER*4 sizeInv !size of the inverse matrix
      real*8 :: real_inva(sizeinv,sizeinv) 
      real*8 :: imag_inva(sizeinv,sizeinv) 
      
      INTEGER*4 idx, jdx, maxmem



      DATA maxfct /1/, mnum /1/
      allocate( pt(64), idum(0), ddum(0) )
  
      do idx = 1, 64
        pt(idx)%DUMMY = 0
      end do
      pa => a
      pia => ia
      pja => ja
! initialize error flag
      error = 0


#ifdef DEBUG
      write(0,*)
      write(0,*) 'The input matrix: (getPartialInv)'
      idx = 0
      do jdx = 1, nonzerosa
        if (jdx >= ia(idx+1)) then
          idx=idx+1
        end if

        if (jdx < ia(idx+1)) then  
          !write(0,*) 'row', idx, 'col', ja(jdx), 'element of a = ', a(jdx)
        end if
      end do
      write(0,*) 'The size of the partial inverse to be calculated:', sizeinv
#endif

!
!  .. Setup Pardiso control parameters und initialize the solvers     
!     internal adress pointers. This is only necessary for the FIRST   
!     call of the PARDISO solver.                                     
!       
!     choose the matrix type for pardiso factorization
      mtype     = 3  ! structurally symmetric complex matrix

!  .. PARDISO license check and initialize solver
      call pardisoinit(pt, mtype, iparm)




!     Reordering and Factorization, This step also allocates
!     all memory that is necessary for the factorization
#ifdef DEBUG      
      msglvl    = 1      ! with statistical information
      iparm(27) = 1      ! perform matrix check
#else
      msglvl    = 0      ! without statistical information
      iparm(27) = 0      ! perform matrix check
#endif
      iparm(6) = 0  ! The solution is written into x  
      iparm(28) = 0 ! double real/complex

#ifdef MIC
      iparm(2) = 3 ! parallel OpenMP version of the nested dissection
      iparm(24) = 1 ! Intel MKL PARDISO uses a two-level factorization algorithm
      iparm(25) = 0 ! Intel MKL PARDISO uses a parallel algorithm for the solve step.
      iparm(34) = 0 ! 
      iparm(60) = 0 ! in-core (IC) mode
#else
      iparm(2) = 2 ! parallel OpenMP version of the nested dissection
      iparm(24) = 0 ! Intel MKL PARDISO uses a two-level factorization algorithm
      iparm(25) = 0 ! Intel MKL PARDISO uses a parallel algorithm for the solve step.
      iparm(34) = 0 ! 
      iparm(60) = 0 ! in-core (IC) mode
#endif

      phase     = 11  ! Analysis

      ! setting the appropriate number of right hand sides.
#ifdef MIC
      nrhs = maxnrhs      
      if (nrhs.GT.sizeinv) then
        nrhs = sizeinv
      end if
#else
      nrhs = 1
#endif

#ifdef DEBUG
      write(0,*) 'number of right sides: ', nrhs
#endif

      

      CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, pa, pia, pja, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error) 
      maxmem = (iparm(16)+iparm(17))

#ifdef DEBUG
      write(0,*) 'maxmem: ', maxmem, ' kbytes'
      write(0,*) 'maxmem total: ', maxmem+sizea*nrhs*32/1000, ' kbytes'
#endif

      IF ((maxmem+sizea*nrhs*32/1000).GT.maxmem) then
        nrhs = floor(real(maxmem-maxmem)/real(sizea*32)*1000)
        write(0,*) 'number of right sides: ', nrhs
        IF (nrhs.EQ.0) THEN
          error = -2
          RETURN
        ENDIF
      ENDIF

      phase     = 22  ! factorization  
      CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, pa, pia, pja, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error) 

#ifdef DEBUG
      write(0,*) 'Reordering and factorization completed ... '
#endif
      IF (error .NE. 0) THEN
         write(0,*) 'The following ERROR was detected at the factorization step: ', error
        RETURN
      ENDIF 



!     Set parameters to calculate only the last sizeInv components of the solution
      iparm(31) = 3
      phase = 33 ! solver
      allocate( perm(sizea) )
      perm = 0
      do idx = 1, sizeinv
        perm(sizea-sizeinv+idx) = 1
      end do

!     Calculating the chunk of the solution to be calculated at once
      if (modulo(sizeinv, nrhs).EQ.0) then
        numiterations = sizeinv/nrhs
      else
        numiterations = floor(real(sizeinv)/real(nrhs)) + 1
      end if

#ifdef DEBUG
      write(0,*) 'Number of iterations to obtain the right solutions ', numiterations
#endif

!     Allocating the right hand side
      allocate( pb(sizea, nrhs) )
      allocate( px(sizea, nrhs) )

      do idxiterations = 1, numiterations

        offset = (idxiterations-1)*nrhs

        if ((sizeinv-offset).LT.nrhs) then
          nrhs = sizeinv-offset
          deallocate(pb, px) 
          allocate( pb(sizea, nrhs) )
          allocate( px(sizea, nrhs) )
        end if

!       Creating the right hand side
        pb = (0.d0, 0.d0)
        px = (0.d0, 0.d0)
        do idx = 1, nrhs
          pb(sizea-sizeinv+offset+idx, idx) = (1.d0, 0.d0)
        end do
!       Solver        
        CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, pa, pia, pja, &
                   perm, nrhs, iparm, msglvl, pb, px, error)  

        IF (error .NE. 0) THEN
          exit
        ENDIF 

        do idx = 1, nrhs 
          do jdx = 1,sizeinv
            real_inva(offset+idx,jdx) = real(px(sizea-sizeinv+jdx, idx))
            imag_inva(offset+idx,jdx) = aimag(px(sizea-sizeinv+jdx, idx))
          end do
        end do  

      end do

      deallocate(pb, px, perm) 


      IF (error .NE. 0) THEN
        write(0,*) 'The following ERROR was detected at the solution step: ', error
        RETURN
      ENDIF 

!  
!   Termination and release of memory
      phase     = -1           ! release internal memory
      CALL pardiso (pt, maxfct, mnum, mtype, phase, sizea, ddum, &
                   idum, idum, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error) 

      deallocate(idum, ddum, pt)

      IF (error .NE. 0) THEN
        write(0,*) 'The following ERROR was detected at releasing the memory: ', error
        RETURN
      ENDIF 

  end



end module