dgetPartialInv.F

Go to the documentation of this file.

#include "fintrf.h"
C======================================================================
C     Gateway routine to call the gatPartialInv.F subroutine.
C     (See: http://www.netlib.org/lapack/explore-html/d9/d52/dggev_8f.html for details)
C     Copyright (C) 2009-2015 Peter Rakyta, Ph.D.
C
C     This program is free software: you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation, either version 3 of the License, or
C     (at your option) any later version.
C 
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License
C     along with this program.  If not, see http://www.gnu.org/licenses/.
C     
C======================================================================
C     Gateway routine
      subroutine mexfunction(nlhs, plhs, nrhs, prhs)

C     Declarations
      implicit none

      interface
        subroutine getpartialinv 
     1          (asize, a, ia, ja, sizeinv, inva)
          INTEGER*4 asize !matrix A is of size asize x asize
          INTEGER*4, ALLOCATABLE :: ia(:), ja(:) 
          real*8, ALLOCATABLE ::  a(:)    
          INTEGER*4 sizeInv !size of the Schur complement to be calculated
          real*8, ALLOCATABLE ::  inva(:,:)    
        end subroutine
      end interface

C     mexFunction arguments:
      mwpointer plhs(*), prhs(*)
      integer nlhs, nrhs


C     Function declarations:
      mwsize mxgetnzmax
      mwpointer mxgetjc
      mwpointer mxgetir

C     Arguments for mxCreateNumericArray
      integer*4 complexflag
      mwsize ndim

      mwpointer mxgetpr
      mwpointer mxcreatedoublematrix
      mwpointer mxgetm, mxgetn
      integer*4 mxIsDouble
      integer*4 mxIsComplex
      integer*4 mexPrintf

  

C     Array information:
      mwsize nonzerosa !number of nonzeros elements in A
      mwsize sizeam, sizean !size of the matrix A (rows, and cols)  
      integer*4 sizeInv !Size of the Schur complement to be calculated

C     Arguments for computational routine:  
C     fortran representation of the input matrix (CSR format)
      integer*8, allocatable :: ia_int8(:), ja_int8(:) 
      integer*4, allocatable :: ia(:), ja(:)
      real*8, allocatable ::  a(:) 
      real*8, allocatable ::  inva(:,:) 
      mwpointer mxinva
        
      integer i,j

#ifdef DEBUG
C     Debug parameters
      integer*4 tmp
      character*800 line
#endif


C-----------------------------------------------------------------------
C     Check for proper number of arguments. 
      if(nrhs .ne. 2) then
         call mexerrmsgidandtxt ('MATLAB:timestwo:nInput',
     +                           'Two inputs required.')
      elseif(nlhs .gt. 1) then
         call mexerrmsgidandtxt ('MATLAB:timestwo:nOutput',
     +                           'Too many output arguments.')
      endif

C     Validate inputs
      
C     Check that the input is a double.
      if ((mxisdouble(prhs(1)) .eq. 0)) then
         call mexerrmsgidandtxt ('MATLAB:timestwo:NonDouble',
     +                           'Input must be a double.')
      endif
      
C     Check if inputs are complex.
      if ((mxiscomplex(prhs(1)) .ne. 0) ) then
         call mexerrmsgidandtxt ('MATLAB:convec:NonComplex',
     +                           'Inputs must be real.')
      endif    

C     Get the size of the input array. (MATLAB)
      nonzerosa = mxgetnzmax(prhs(1))
      sizeam = mxgetm(prhs(1))
      sizean = mxgetn(prhs(1))


C     Allocating fortran arrays
      allocate(ia_int8(sizean+1))
      allocate(ja_int8(nonzerosa))
      allocate(ia(sizean+1))
      allocate(ja(nonzerosa))
      allocate(a(nonzerosa))
C
C     Check the allocated arrays
      if (.not. ( allocated(ia_int8) .and.
     + allocated(ja_int8) .and.
     + allocated(ia) .and.
     + allocated(ja) .and.
     + allocated(a) )) then
         call mexerrmsgidandtxt ('MATLAB:dgetSchur:BadAllocation',
     +       'The allocation of variables was unsuccesfull.')
      endif
C     Copy data into the fortran arrays
C     MATLAB does CSC format and PARDISO CSR format for the sparse matrices! Need to transpose.   
      call mxcopyptrtointeger8( 
     +                   mxgetjc(prhs(1)), ia_int8, sizean+1 )
      call mxcopyptrtointeger8( 
     +                   mxgetir(prhs(1)), ja_int8, nonzerosa )
      call mxcopyptrtoreal8( mxgetpr(prhs(1)), a, nonzerosa )
C     Fortran starts the indexes from 1
      ia = int(ia_int8, 4) + 1
      ja = int(ja_int8, 4) + 1
      deallocate(ia_int8)
      deallocate(ja_int8)

C     choose the matrix type for pardiso factorization
C     getting the size of the partial inverse
      ndim = 1
      call mxcopyptrtointeger4( mxgetpr(prhs(2)), sizeinv, ndim )


C     Do the pardiso calulations to get the partial inverse
      ALLOCATE(inva(sizeinv, sizeinv))
      CALL getpartialinv(int(sizean, 4), a, ia, ja, sizeinv, inva)

C     deallocate the fortran matrix A
      deallocate(ia)
      deallocate(ja)
      deallocate(a)


C     Allocating memory for the Schur complement  
#ifdef DEBUG
      write(*,*) 'diagonal Elements of the partial solution'
      do i = 1, sizeinv
        write(line,*) inva(i,i), new_line('A')
        tmp = mexprintf( line )
      end do
#endif
 
C     Converting the partial inverse to matlab format
      ndim = int(sizeinv, 8)
      complexflag = 0;
      plhs(1) = mxcreatedoublematrix(ndim, ndim, 
     + complexflag) ! for invA
      mxinva = mxgetpr(plhs(1))
      call mxcopyreal8toptr(inva, mxinva,ndim*ndim) 


      
C.. Deallocating memory for the partial inverse

      if ( allocated(inva) ) then
        DEALLOCATE(inva)
      endif      


      return
      end