! SPDX-Identifier: MIT
#:include "common.fypp"
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
#:set I_KINDS_TYPES = list(zip(INT_KINDS, INT_TYPES, INT_KINDS))
#:set RCI_KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES + I_KINDS_TYPES
submodule (stdlib_io_mm) stdlib_io_mm_load
use stdlib_error, only : error_stop
use stdlib_strings, only : to_string, starts_with
use stdlib_str2num, only: to_num_from_stream
use stdlib_kinds
implicit none
enum, bind(c)
enumerator :: MF_array = 1
enumerator :: MF_coordinate = 2
end enum
enum, bind(c)
enumerator :: MQ_real = 1
enumerator :: MQ_integer = 2
enumerator :: MQ_complex = 3
enumerator :: MQ_pattern = 4
end enum
enum, bind(c)
enumerator :: MS_general = 1
enumerator :: MS_symmetric = 2
enumerator :: MS_skew_symmetric = 3
enumerator :: MS_hermitian = 4
end enum
integer(int8), parameter :: LF = 10, CR = 13, PP=iachar('%')
contains
#:for k, t, s in RCI_KINDS_TYPES
module subroutine load_mm_dense_${s}$(filename, matrix, iostat, iomsg)
!> Name of the Matrix Market file to load from
character(len=*), intent(in) :: filename
!> Matrix to be loaded from the Matrix Market file
${t}$, allocatable, intent(out) :: matrix(:,:)
!> Error status of loading, zero on success
integer, intent(out), optional :: iostat
!> Associated error message in case of non-zero status code
character(len=:), allocatable, intent(out), optional :: iomsg
! Internal variables
type(mm_header_type) :: header
integer :: u , fsze, err, eol_position
integer :: nrows, ncols, i, j
integer(int8) :: stat
character(:), allocatable, target :: ff
character(len=:), pointer :: ffp
#:if t.startswith('complex')
real(${k}$) :: mold, val_r, val_i
#:else
${t}$ :: mold
#:endif
if (present(iostat)) iostat = 0
if (present(iomsg)) iomsg = ''
stat = 0
err = 0
!-----------------------------------------------------------------------------
! Open file for regular reading
open( newunit = u , file=filename, status = 'old' , access='stream', action="read", iostat=err )
if( err /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err,&
message = 'Error opening matrix market file')
return
end if
err = 1
!-----------------------------------------
! Load file in a single string
inquire(unit=u, size=fsze)
allocate(character(fsze) :: ff)
read(u) ff
ffp => ff(1:)
close(u)
!-----------------------------------------
! Read header
call read_mm_header(ffp, header, err)
if( err /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err,&
message = 'Error reading mm header')
return
end if
if( header%format /= MF_array ) then
err = 2
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err, &
message = 'warning: a dense matrix is expected for the current file')
return
end if
!-----------------------------------------
! Skip comments
eol_position = shift_to_eol(ffp)
ffp => ffp(eol_position+1:)
do while( iachar(ffp(1:1))==PP )
eol_position = shift_to_eol(ffp)
ffp => ffp(eol_position+1:)
end do
!-----------------------------------------
! Read matrix dimensions
nrows = to_num_from_stream(ffp, nrows, stat)
if( stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat),&
message = 'Error reading number of rows')
return
end if
ncols = to_num_from_stream(ffp, ncols, stat)
if( stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat),&
message = 'Error reading number of columns')
return
end if
!-----------------------------------------
! Read actual matrix data
allocate(matrix(nrows, ncols), stat=err)
matrix = 0
if( err /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err,&
message = 'Error allocating matrix')
return
end if
read_vals: block
if(header%symmetry==MS_general) then
do j = 1, ncols
do i = 1, nrows
#:if t.startswith('complex')
val_r = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
val_i = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
matrix(i,j) = cmplx(val_r, val_i, kind = ${k}$)
#:else
matrix(i,j) = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
#:endif
end do
end do
else
do j = 1, ncols
do i = j, nrows
! Keep diagonal elements as zero incase of skew-symmetric cases
if(header%symmetry==MS_skew_symmetric .and. i==j) cycle
#:if t.startswith('complex')
val_r = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
val_i = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
matrix(i,j) = cmplx(val_r, val_i, kind = ${k}$)
#:else
matrix(i,j) = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
#:endif
! Assign transpose of the current element
matrix(j, i) = matrix(i, j)
if(header%symmetry==MS_skew_symmetric) matrix(j, i) = -matrix(j, i)
#:if t.startswith('complex')
if(header%symmetry==MS_hermitian) matrix(j, i) = conjg(matrix(j, i))
#:endif
end do
end do
end if
end block read_vals
if( stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat),&
message = 'Error reading matrix value')
return
end if
end subroutine
#:endfor
#:for k, t, s in RCI_KINDS_TYPES
module subroutine load_mm_coo_${s}$(filename, index, data, iostat, iomsg)
!> Name of the Matrix Market file to load from
character(len=*), intent(in) :: filename
!> Matrix indices in COO format read from the file:
!> index(1, :) = row indices, index(2, :) = column indices
integer, allocatable, intent(out) :: index(:, :)
!> Nonzero matrix values corresponding to each (row, column) index pair
${t}$, allocatable, intent(out) :: data(:)
!> Error status of loading, zero on success
integer, intent(out), optional :: iostat
!> Associated error message in case of non-zero status code
character(len=:), allocatable, intent(out), optional :: iomsg
! Internal variables
type(mm_header_type) :: header
integer :: u , fsze, err, eol_position
integer :: i, nnz, adr
integer(int8) :: stat
character(:), allocatable, target :: ff
character(len=:), pointer :: ffp
integer, allocatable :: rows(:), cols(:)
${t}$, allocatable :: vals(:)
integer :: nrows, ncols, n_diag
#:if t.startswith('complex')
real(${k}$) :: mold, val_r, val_i
#:else
${t}$ :: mold
#:endif
if (present(iostat)) iostat = 0
if (present(iomsg)) iomsg = ''
stat = 0
err = 0
!-----------------------------------------------------------------------------
! Open file for regular reading
open( newunit = u , file=filename, status = 'old' , access='stream', action="read", iostat=err )
if( err /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err,&
message = 'Error opening matrix market file')
return
end if
err = 1
!-----------------------------------------
! Load file in a single string
inquire(unit=u, size=fsze)
allocate(character(fsze) :: ff)
read(u) ff
ffp => ff(1:)
close(u)
!-----------------------------------------
! Read header
call read_mm_header(ffp, header, err)
if( err /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err,&
message = 'Error reading mm header')
return
end if
if( header%format /= MF_coordinate ) then
err = 2
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = err, &
message = 'warning: a coordinate matrix is expected for the current file')
return
end if
!-----------------------------------------
! Skip comments
eol_position = shift_to_eol(ffp)
ffp => ffp(eol_position+1:)
do while( iachar(ffp(1:1))==PP )
eol_position = shift_to_eol(ffp)
ffp => ffp(eol_position+1:)
end do
!-----------------------------------------
! Read matrix dimensions
nrows = to_num_from_stream(ffp, nrows, stat)
if( stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat),&
message = 'Error reading nrows')
return
end if
ncols = to_num_from_stream(ffp, ncols, stat)
if( stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat),&
message = 'Error reading ncols')
return
end if
nnz = to_num_from_stream(ffp, nnz, stat)
if( stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat),&
message = 'Error reading nnz')
return
end if
!-----------------------------------------
! Allocate temporary arrays to hold the file data
allocate(rows(nnz))
allocate(cols(nnz))
if(header%qualifier/=MQ_pattern) allocate(vals(nnz))
!-----------------------------------------
! Read actual matrix data and store inside temporary arrays
n_diag = 0
read_vals: block
if(header%qualifier==MQ_pattern) then
do i = 1, nnz ! read entries from file
rows(i) = to_num_from_stream(ffp, rows(i), stat)
if(stat/=0) exit read_vals
cols(i) = to_num_from_stream(ffp, cols(i), stat)
if(stat/=0) exit read_vals
if(rows(i) == cols(i)) n_diag = n_diag + 1
if(stat/=0) exit read_vals
end do
else
do i = 1, nnz ! read entries from file
rows(i) = to_num_from_stream(ffp, rows(i), stat)
if(stat/=0) exit read_vals
cols(i) = to_num_from_stream(ffp, cols(i), stat)
if(stat/=0) exit read_vals
if(rows(i) == cols(i)) n_diag = n_diag + 1
#:if t.startswith('complex')
val_r = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
val_i = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
vals(i) = cmplx(val_r, val_i, kind = ${k}$)
#:else
vals(i) = to_num_from_stream(ffp, mold, stat)
if(stat/=0) exit read_vals
#:endif
end do
end if
end block read_vals
if(stat /= 0 ) then
call mm_fail_process(iostat = iostat, iomsg = iomsg, code = int(stat), &
message = 'Error reading the Matrix Market coordinate data')
return
end if
!-----------------------------------------
! check storage hypothesis
if(header%symmetry == MS_symmetric .or. header%symmetry == MS_hermitian) then
allocate(index(2, 2*nnz-n_diag))
if(header%qualifier/=MQ_pattern) allocate(data(2*nnz-n_diag))
else if(header%symmetry == MS_skew_symmetric) then
allocate(index(2, 2*nnz))
if(header%qualifier/=MQ_pattern) allocate(data(2*nnz))
else
allocate(index(2, nnz))
if(header%qualifier/=MQ_pattern) allocate(data(nnz))
end if
!-----------------------------------------
! Fill in matrix entries from temporary arrays
if(header%qualifier==MQ_pattern) then
do i = 1, nnz
index(1,i) = rows(i)
index(2,i) = cols(i)
end do
else
do i = 1, nnz
index(1,i) = rows(i)
index(2,i) = cols(i)
data(i) = vals(i)
end do
end if
if(allocated(rows)) deallocate(rows)
if(allocated(cols)) deallocate(cols)
if(allocated(vals)) deallocate(vals)
!-----------------------------------------
! Fill in symmetric entries if needed
if(header%symmetry==MS_general) return
adr = 1
if(header%qualifier==MQ_pattern) then
do i = 1, nnz
if(index(1,i)==index(2,i)) cycle
index(1,nnz+adr) = index(2,i)
index(2,nnz+adr) = index(1,i)
adr = adr + 1
end do
else
do i = 1, nnz
if(index(1,i)==index(2,i)) cycle
index(1,nnz+adr) = index(2,i)
index(2,nnz+adr) = index(1,i)
data(nnz+adr) = data(i)
if(header%symmetry==MS_skew_symmetric) data(nnz+adr) = -data(i)
#:if t.startswith('complex')
if(header%symmetry==MS_hermitian) data(nnz+adr) = conjg(data(i))
#:endif
adr = adr + 1
end do
end if
end subroutine
#:endfor
subroutine read_mm_header(ffp, header, err)
character(len=:), intent(inout), pointer :: ffp
type(mm_header_type), intent(out) :: header
integer, intent(out) :: err
!----------------------------------------------
err = 0
if( .not. starts_with(ffp, "%%MatrixMarket ") ) then
err = 1
return
end if
ffp => ffp(16:)
! Read object type: matrix
if( .not. starts_with(ffp, "matrix ") ) then
err = 1
return
end if
ffp => ffp(8:)
header%object = 1 ! matrix
! Read format type: coordinate or array
if( starts_with(ffp, "arr") ) then
ffp => ffp(7:) ! array
header%format = MF_array
else if( starts_with(ffp, "coo") ) then
ffp => ffp(12:) ! coordinate
header%format = MF_coordinate
else
err = 1
return
end if
! Read first qualifier: real, complex, integer, pattern (sparse)
if( starts_with(ffp, "real") ) then
ffp => ffp(6:) ! real
header%qualifier = MQ_real
else if( starts_with(ffp, "complex") ) then
ffp => ffp(9:) ! complex
header%qualifier = MQ_complex
else if( starts_with(ffp, "integer") ) then
ffp => ffp(9:) ! integer
header%qualifier = MQ_integer
else if( starts_with(ffp, "pattern") ) then
ffp => ffp(9:) ! pattern
header%qualifier = MQ_pattern
else
err = 1
return
end if
! Read second qualifier: general, symmetric, skew-symmetric, hermitian
if( starts_with(ffp, "general") ) then
ffp => ffp(9:) ! general
header%symmetry = MS_general
else if( starts_with(ffp, "symmetric") ) then
ffp => ffp(11:) ! symmetric
header%symmetry = MS_symmetric
else if( starts_with(ffp, "skew-symmetric") ) then
ffp => ffp(16:) ! skew-symmetric
header%symmetry = MS_skew_symmetric
else if( starts_with(ffp, "hermitian") ) then
ffp => ffp(11:) ! hermitian
header%symmetry = MS_hermitian
else
err = 1
return
end if
end subroutine
elemental function shift_to_eol(s) result(p)
!! move string to position of the next end-of-line character
character(*),intent(in) :: s !! character chain
integer :: p !! position
!----------------------------------------------
p = 1
do while( p<len(s) .and. .not.(iachar(s(p:p))==LF .or. iachar(s(p:p))==CR) )
p = p + 1
end do
! If CRLF, move to LF
if (p < len(s)) then
if (iachar(s(p:p)) == CR .and. iachar(s(p+1:p+1)) == LF) then
p = p + 1
end if
end if
end function
subroutine mm_fail_process(code, message, iostat, iomsg)
integer, intent(out), optional :: iostat
character(len=:), allocatable, intent(out), optional :: iomsg
integer, intent(in) :: code
character(*), intent(in) :: message
if (present(iostat)) iostat = code
if (present(iomsg)) then
iomsg = message
else
call error_stop(message)
end if
end subroutine
end submodule stdlib_io_mm_load
