MatrixIO.F90

Go to the documentation of this file.

!#include "fintrf.h"
!======================================================================
!     Test program to test the 
!     Copyright (C) 2018 Peter Rakyta, Ph.D.
!
!     This program is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Publi! 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/.
!     
!======================================================================
!     Gateway routine
      PROGRAM matrixio

      use getpartialinv


!     Declarations
      implicit none


      integer::narg,cptarg !#of arg & counter of arg
      character(len=100):: name !Arg name
      logical lookforfilename, lookforinvdim
      logical lookforoutput
      logical fileexist
      character(len=100) :: filename, output

!     Array information:
      integer*4 nonzerosa !number of nonzeros elements in A
      integer*4 sizeam, sizean !size of the matrix A (rows, and cols) 
      integer*4 sizeinv !the size of the partial inverse to be calculated
      character(len=20) sizeinv_string !the size of the partial inverse -- string from the input
      character(len=8) out_status

!     Arguments for computational routine:  
      integer*4, allocatable :: ia(:), ja(:)
#ifdef COMP
      COMPLEX*16, ALLOCATABLE ::  a(:)    
#else
      real*8, ALLOCATABLE ::  a(:)    
#endif       

#ifdef COMP
      COMPLEX*16, allocatable ::  inva(:,:)  
#else
      real*8, allocatable ::  inva(:,:)  
#endif
!      
        
      integer i,j

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

!     https://genomeek.wordpress.com/2012/02/09/fortran-command-line-arguments/
      !Check if any arguments are found
      lookforfilename = .false.
      lookforinvdim = .false.
      lookforoutput = .false.
      filename = ''
      output = ''
      sizeinv = 0
      
      narg=command_argument_count()
      !Loop over the arguments
      if(narg>0)then
      !loop across options
        do cptarg=1,narg
          call get_command_argument(cptarg,name)
          select case(adjustl(name))
          case("--help","-h")
            write(0,*)"This is program calculates a " , &
         "partial inverse of a sparse matirix. Usage: ./MatrixIO ",  &
         "--filename FILENAME --invDim INTEGER --Output FILENAME"
            stop
          case("--filename","-F", "--Filename")
            lookforfilename = .true.
          case("--output","-o", "--Output")
            lookforoutput = .true.
          case("--invDim","-I")
            lookforinvdim = .true.
          case default
            if (lookforfilename) then
              filename=adjustl(name) !assign a value to pedfile
              lookforfilename = .false.
            elseif (lookforoutput) then
              output=adjustl(name) !assign a value to pedfile
              lookforoutput = .false.
            elseif (lookforinvdim) then 
              sizeinv_string = trim(adjustl(name))
              READ( sizeinv_string, * ) sizeinv
              lookforinvdim = .false.
            else
              write(0,*)"Option ",adjustl(name),"unknown"
            end if
          end select
        end do
      end if

      if (output.eq.'') then
        output = filename
      end if

#ifdef DEBUG
      write(0,*) 'Input file name:', trim(filename)
      write(0,*) 'Output file name:', trim(output)
#endif

!     Testing the existance of the input file containing the sparse matrix A    
      inquire(file=trim(filename),exist=fileexist)!check if it exist
      if(.not.fileexist)then
        write(0,*)"file ",trim(filename)," not found!"
        stop
      endif

!     Opening the file to read in data   
      OPEN(unit=6, file=trim(filename),status='OLD',form='UNFORMATTED')

      READ(6,end=999,err=1000) sizean
      READ(6,end=999,err=1000) nonzerosa

      allocate(ia(sizean+1))
      allocate(ja(nonzerosa))
      allocate(a(nonzerosa))

      READ(6,end=999,err=1000) ia(:)
      READ(6,end=999,err=1000) ja(:)
      READ(6,end=999,err=1000) a(:)



!     Closing the input file
      CLOSE(unit=6)

#ifdef DEBUG
      write(0,*) 'Size of the partial inverse:', sizeinv
      write(0,*) 'The size of the input MATRIX:', sizean
      write(0,*) 'Nonzeros in the input matrix:', nonzerosa     
      write(0,*) 'The input MATRIX:'
      i = 0
      do j = 1, size(a)
        if (j >= ia(i+1)) then
          i=i+1
        end if

        if (j < ia(i+1)) then  
          !write(0,*) 'row', i, 'col', ja(j), 'element of A = ', a(j)
        end if
      end do

#endif




!     Calculating the partial inverse
      ALLOCATE(inva(sizeinv, sizeinv))

      CALL dgetpartialinv(sizean, nonzerosa, a, ia, ja, sizeinv, inva)

!     deallocate the fortran matrix A
      deallocate(ia)
      deallocate(ja)
      deallocate(a)
#ifdef DEBUG
      do i=1, sizeinv
        do j=1, sizeinv
          write(0,*) 'row', i, 'col', j, ' element of A = ', inva(i,j)
        end do
      end do
#endif

      
!     Deallocating memory for the partial inverse  

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


      stop

 999  write(*,'(/"End-of-file when i = ",I5)') i
      stop

 1000 write(*,'(/"ERROR reading when i = ",I5)') i
      stop
      end