dgetSchur.F

Go to the documentation of this file.

#include "fintrf.h"
C======================================================================
C     Gateway routine to call the zggev.f function from LAPACK package.
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 getschur 
     1          (asize, a, ia, ja, ssize, s, is, js, mtype)
          INTEGER*4 asize !matrix A is of size asize x asize
          INTEGER*4, ALLOCATABLE :: ia(:), ja(:) 
          real*8, ALLOCATABLE ::  a(:)           
          INTEGER*4 ssize !size of the Schur complement to be calculated
          INTEGER, ALLOCATABLE :: iS(:), jS(:) ! row and culomn indices of the Schur complement
          real*8, ALLOCATABLE :: s(:) !nonzeros elements in the Schur complement
          INTEGER mtype
        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 nonzerosS !number of nonzeros elements in S
      integer*4 sizeS !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(:) 
C     fortran representation of the Schur complement matrix (CSR format)
      integer, allocatable :: iS(:), jS(:) 
      real*8, allocatable ::  s(:) 
      real*8, allocatable ::  s_full(:,:) 
      mwpointer mxs_full
        
C     Parameters for pardiso
      INTEGER mtype
      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
      mtype     = 1  ! structurally symmetric real matrix
C     getting the size of the Schur complement
      ndim = 1
      call mxcopyptrtointeger4( mxgetpr(prhs(2)), sizes, ndim )


C     Do the pardiso calulations to get the Schur complement
C     Passing deallocated arrays might cause troubles
      ALLOCATE(is(0))
      ALLOCATE(js(0))
      ALLOCATE(s(0))
      CALL getschur(int(sizeam, 4), a, ia, ja, 
     +   sizes, s, is, js, mtype)

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


#ifdef DEBUG
      write(line,*) new_line('A')
      tmp = mexprintf(line)
      write(line,*) 'The Calculated Schur complement:', new_line('A')
      tmp = mexprintf(line)
      i = 0
      do j = 1, size(s)
        if (j >= is(i+1)) then
          i=i+1
        end if

        if (j < is(i+1)) then  
          write(line,*) 'row', i, 'col', js(j),
     + 'element of S = ', s(j), new_line('A')
          tmp = mexprintf(line)
        end if
      end do

#endif

C     creating full S matrix
      ALLOCATE(s_full(sizes, sizes))
      s_full = 0.d0
      i = 0
      do j = 1, size(s)
        if (j >= is(i+1)) then
          i=i+1
        end if

        if (j < is(i+1)) then 
C     MATLAB does CSC format and PARDISO CSR format for the sparse matrices! Need to transpose back.          
          s_full(js(j), i) = s(j)
        end if
      end do
      ndim = int(sizes, 8)
      complexflag = 0;
      plhs(1) = mxcreatedoublematrix(ndim, ndim, 
     + complexflag) ! for S_full
      mxs_full = mxgetpr(plhs(1))
      call mxcopyreal8toptr(s_full, mxs_full,ndim*ndim) 




      
C.. Deallocating memory for the Schur complement
      if ( allocated(is) ) then
        DEALLOCATE(is)
      endif

      if ( allocated(js) ) then
        DEALLOCATE(js)
      endif

      if ( allocated(s) ) then
        DEALLOCATE(s)
      endif     

      if ( allocated(s_full) ) then
        DEALLOCATE(s_full)
      endif      


      return
      end