export_complex_sparse.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

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


C     Function declarations:
      mwsize mxgetnzmax
      mwpointer mxgetjc
      mwpointer mxgetir
      mwpointer mxcreatestring
      mwpointer mxgetstring, getstring_stat 
      integer mxIsChar

C     Arguments for copying data between fortran and MATLAB
      integer*4 complexflag

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

      logical fileExist
      character(len=100) :: filename 
      character(len=7) :: file_status

C     Array information:
      mwsize nonzerosa !number of nonzeros elements in A
      mwsize sizeam, sizean !size of the matrix A (rows, and cols) 
      mwpointer mrows, ncols !size of the input string containing the filename
      integer*4 strlen ! the length of the input string
      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(:)
      complex*16, allocatable ::  a(:) 
        
C     Parameters for pardiso
      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:export_real_sparse:nInput',
     +                           'Two inputs required.')
      elseif(nlhs .gt. 0) then
         call mexerrmsgidandtxt ('MATLAB:export_real_sparse: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:export_real_sparse:NonDouble',
     +                           'Input must be a double.')
      endif
      
C     Check if input matrix must be complex
      if ((mxiscomplex(prhs(1)) .ne. 1) ) then
         call mexerrmsgidandtxt ('MATLAB:convec:NonComplex',
     +                           'Inputs must be complex.')
      endif    

C     The input filename must be a string.
       if(mxischar(prhs(2)) .ne. 1) then
         call mexerrmsgidandtxt ('MATLAB:export_real_sparse:NonString',
     +                 'Input filename must be a string.')
       endif

C     Get the size of the input string.
       mrows = mxgetm(prhs(2))
       ncols = mxgetn(prhs(2))

C     Get the length of the input string and validate.
       strlen = int(mrows*ncols, 4)
       if (strlen .gt. len(filename)) then
         call mexerrmsgidandtxt ('MATLAB:export_real_sparse:NonDouble',
     +             'input filename is greater than max str size.')   
       endif

C     Get the input filename.
      getstring_stat = 
     + mxgetstring(prhs(2), filename, int(len(filename), 8))

#ifdef DEBUG
      write(line,*) 'Input file name: ',filename, new_line('A')
      tmp = mexprintf(line)
#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 mxcopyptrtocomplex16( mxgetpr(prhs(1)), mxgetpi(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     Opening the file to export data 
      inquire(file=trim(filename),exist=fileexist)!check if it exist
      if(.not.fileexist)then
        file_status = 'NEW'
      else
        file_status = 'REPLACE'
      endif   
      OPEN(unit=6, file=trim(filename),status=file_status
     + ,form='UNFORMATTED')

      WRITE(6) int(sizean, 4)
      WRITE(6) int(nonzerosa, 4)
      WRITE(6) ia(:)
      WRITE(6) ja(:)
      WRITE(6) a(:)


C     Closing the file
      CLOSE(unit=6)

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

      return
      

      end