diff --git a/doc/specs/stdlib_io.md b/doc/specs/stdlib_io.md index 46b6952e1..307baa89a 100644 --- a/doc/specs/stdlib_io.md +++ b/doc/specs/stdlib_io.md @@ -321,3 +321,136 @@ Exceptions trigger an `error stop` unless the optional `err` argument is provide {!example/io/example_get_file.f90!} ``` +## Matrix Market Format I/O + +### Status + +Experimental + +### Description + +The Matrix Market I/O module provides support for reading and writing matrices in the Matrix Market format, a simple ASCII format for sparse and dense matrices developed at NIST. The format supports real, complex, and integer matrices with various symmetry properties. + +### `load_mm` - load a matrix from Matrix Market file + +Loads a 2D matrix from a Matrix Market format file. Symmetric matrices are expanded to full storage. + +#### Syntax + +- To load a matrix of `array` format: + +`call ` [[stdlib_io_mm(module):load_mm(interface)]] `(filename, matrix [, iostat, iomsg])` + +- To load a matrix of `coordinate` format: + +`call ` [[stdlib_io_mm(module):load_mm(interface)]] `(filename, index, data [, iostat, iomsg])` + +#### Arguments + +**`Array` format** + +`call ` [[stdlib_io_mm(module):load_mm(interface)]] `(filename, matrix [, iostat, iomsg])` + +- `filename`: Shall be a character expression containing the Matrix Market file name to read from. + +- `matrix`: Shall be an allocatable rank-2 array of type `real`, `complex`, or `integer` that will contain the loaded matrix. + +- `iostat` (optional): Shall be a scalar of type `integer` that receives the error status. Zero indicates success. + +- `iomsg` (optional): Shall be an allocatable character string that receives the error message if iostat is non-zero. + +**`Coordinate` format** + +`call ` [[stdlib_io_mm(module):load_mm(interface)]] `(filename, index, data [, iostat, iomsg])` + +- `filename`: Shall be a character expression containing the Matrix Market file name to read from. + +- `index`: Shall be an allocatable rank-2 array of type `integer` that will contain the indices of the loaded matrix. + +- `data`: Shall be an allocatable rank-1 array of type `real`, `complex`, or `integer` that will contain the values of the loaded matrix. + +- `iostat` (optional): Shall be a scalar of type `integer` that receives the error status. Zero indicates success. + +- `iomsg` (optional): Shall be an allocatable character string that receives the error message if iostat is non-zero. + +### `save_mm` - save a matrix to Matrix Market file + +Saves a 2D matrix to Matrix Market format file. + +#### Syntax + +- To save a matrix of `array` format: + +`call ` [[stdlib_io_mm(module):save_mm(interface)]] `(filename, matrix [, comment, format, symmetry, iostat, iomsg])` + +- To save a matrix of `coordinate` format: + +`call ` [[stdlib_io_mm(module):save_mm(interface)]] `(filename, index, data [, comment, format, symmetry, iostat, iomsg])` + +#### Arguments + +**`Array` format** + +`call ` [[stdlib_io_mm(module):save_mm(interface)]] `(filename, matrix [, comment, format, symmetry, iostat, iomsg])` + +- `filename`: Shall be a character expression containing the Matrix Market file name to write to. + +- `matrix`: Shall be a rank-2 array of type `real`, `complex`, or `integer` to save. + +- `comment` (optional): Shall be a character expression containing additional comments for the file header. + +- `format` (optional): Shall be a character expression specifying how entries are written. + +- `symmetry` (optional): Shall be a character expression defining the symmetry of the matrix. Allowed values: `auto`, `general`, `symmetric`, `skew-symmetric`, `hermitian`. + - `auto`: Detects the symmetry automatically, falls back to `general` if no symmetry is found. + - `general`: all entries are stored + - `symmetric` / `hermitian`: only the lower triangle (including diagonal) is stored + - `skew-symmetric`: only the strictly lower triangle (excluding diagonal) is stored + + Default: `general` + +- `iostat` (optional): Shall be a scalar of type `integer` that receives the error status. Zero indicates success. + +- `iomsg` (optional): Shall be an allocatable character string that receives the error message if iostat is non-zero. + +**`Coordinate` format** + +`call ` [[stdlib_io_mm(module):save_mm(interface)]] `(filename, index, data [, comment, format, symmetry, iostat, iomsg])` + +- `filename`: Shall be a character expression containing the Matrix Market file name to write to. + +- `index`: Shall be a rank-2 `integer` array of shape `(2, n)` specifying the indices of the entries. + - index(1,:) shall contain row indices + - index(2,:) shall contain column indices + +- `data`: Shall be a rank-1 array of type `real`, `complex`, or `integer` to save. + - If `size(data) == n`, the values for each index are written. + - If `size(data) == 1`, a `pattern` matrix is written (no explicit values) + +- `comment` (optional): Shall be a character expression containing additional comments for the file header. + +- `format` (optional): Shall be a character expression specifying how entries are written. + +- `symmetry` (optional): Shall be a character expression defining the symmetry of the matrix. Allowed values: `general`, `symmetric`, `skew-symmetric`, `hermitian`. + - `general`: all entries are stored + - `symmetric` / `hermitian`: only the lower triangle (including diagonal) is stored + - `skew-symmetric`: only the strictly lower triangle (excluding diagonal) is stored + +- `iostat` (optional): Shall be a scalar of type `integer` that receives the error status. Zero indicates success. + +- `iomsg` (optional): Shall be an allocatable character string that receives the error message if iostat is non-zero. + +### Matrix Market Format Details + +The Matrix Market format supports: + +- **Object types**: Currently only `matrix` is supported +- **Formats**: `coordinate` (sparse) and `array` (dense) +- **Data types**: `real`, `complex`, `integer`, `pattern` +- **Symmetry**: `general`, `symmetric`, `skew-symmetric`, `hermitian`, `auto` (Currently only supported for `array` format) + +### Example + +```fortran +{!example/io/example_matrix_market.f90!} +``` \ No newline at end of file diff --git a/example/io/CMakeLists.txt b/example/io/CMakeLists.txt index db663f537..fca9bb4f4 100644 --- a/example/io/CMakeLists.txt +++ b/example/io/CMakeLists.txt @@ -3,6 +3,7 @@ ADD_EXAMPLE(fmt_constants) ADD_EXAMPLE(get_file) ADD_EXAMPLE(loadnpy) ADD_EXAMPLE(loadtxt) +ADD_EXAMPLE(matrix_market) ADD_EXAMPLE(open) ADD_EXAMPLE(savenpy) ADD_EXAMPLE(savetxt) diff --git a/example/io/example_matrix_market.f90 b/example/io/example_matrix_market.f90 new file mode 100644 index 000000000..c19dd2c96 --- /dev/null +++ b/example/io/example_matrix_market.f90 @@ -0,0 +1,88 @@ +program example_matrix_market + use stdlib_io_mm, only : load_mm, save_mm + use stdlib_kinds, only : dp + implicit none + + real(dp), allocatable :: matrix(:,:), matrix2(:,:) + integer, allocatable :: index(:,:) + complex(dp), allocatable :: data(:) + character(len=*), parameter :: dense_filename = "example_dense.mtx" + character(len=*), parameter :: sparse_filename = "example_sparse.mtx" + integer :: iostat, i + character(len=:), allocatable :: iomsg + + ! Create a test dense matrix + allocate(matrix(3,3)) + matrix = reshape([1.0_dp, 2.0_dp, 3.0_dp, & + 4.0_dp, 5.0_dp, 6.0_dp, & + 7.0_dp, 8.0_dp, 9.0_dp], [3,3]) + + print *, "=== Dense Matrix Example ===" + print *, "Original dense matrix:" + call print_matrix(matrix) + + ! Save dense matrix to Matrix Market file + call save_mm(dense_filename, matrix, format="ES24.15E2", symmetry="general", iostat=iostat, iomsg=iomsg) + if (iostat /= 0) then + print *, "Error saving dense matrix: ", iomsg + stop 1 + end if + + print *, "Dense matrix saved to ", dense_filename + + ! Load dense matrix from Matrix Market file + call load_mm(dense_filename, matrix2, iostat=iostat, iomsg=iomsg) + if (iostat /= 0) then + print *, "Error loading dense matrix: ", iomsg + stop 1 + end if + + print *, "Loaded dense matrix:" + call print_matrix(matrix2) + + print *, "=== Sparse Matrix Example ===" + ! Create a test sparse matrix + allocate(index(2,6)) + allocate(data(6)) + index(:,1) = [1,1]; data(1) = (10.0_dp, -1.5_dp) + index(:,2) = [2,2]; data(2) = (20.0_dp, 2.0_dp) + index(:,3) = [3,3]; data(3) = (30.0_dp, 3.0_dp) + index(:,4) = [4,4]; data(4) = (40.0_dp, -4.0_dp) + index(:,5) = [1,4]; data(5) = ( 5.0_dp, -7.5_dp) + index(:,6) = [3,1]; data(6) = (15.0_dp, 25.0_dp) + + ! Save sparse matrix to Matrix Market file + call save_mm(sparse_filename, index, data, format="ES24.15E2", symmetry="general", iostat=iostat, iomsg=iomsg) + if (iostat /= 0) then + print *, "Error saving sparse matrix: ", iomsg + stop 1 + end if + + print *, "Sparse matrix saved to ", sparse_filename + + ! Load sparse matrix from Matrix Market file + call load_mm(sparse_filename, index, data, iostat=iostat, iomsg=iomsg) + if (iostat /= 0) then + print *, "Error loading sparse matrix: ", iomsg + stop 1 + end if + + print *, "Loaded sparse matrix (COO format):" + print *, "Data (row, col, value):" + do i = 1, size(data) + print *, index(1,i), index(2,i), data(i) + end do + +contains + + subroutine print_matrix(mat) + real(dp), intent(in) :: mat(:,:) + integer :: i + + do i = 1, size(mat, 1) + print *, mat(i, :) + end do + print * + end subroutine print_matrix + +end program example_matrix_market \ No newline at end of file diff --git a/src/io/CMakeLists.txt b/src/io/CMakeLists.txt index d25526a3e..4587bc698 100644 --- a/src/io/CMakeLists.txt +++ b/src/io/CMakeLists.txt @@ -3,6 +3,9 @@ set(io_fppFiles stdlib_io_npy.fypp stdlib_io_npy_load.fypp stdlib_io_npy_save.fypp + stdlib_io_mm.fypp + stdlib_io_mm_load.fypp + stdlib_io_mm_save.fypp ) set(io_cppFiles diff --git a/src/io/stdlib_io_mm.fypp b/src/io/stdlib_io_mm.fypp new file mode 100644 index 000000000..fe2ae1799 --- /dev/null +++ b/src/io/stdlib_io_mm.fypp @@ -0,0 +1,94 @@ +! 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 + +!> The Matrix Market (MM) format is a simple, human-readable, ASCII format for sparse +!> and dense matrices. The format was developed at NIST (National Institute of Standards +!> and Technology) for the Matrix Market, a repository of test matrices for use in +!> comparative studies of algorithms for numerical linear algebra. +!> +!> For more information, see: https://math.nist.gov/MatrixMarket/formats.html +module stdlib_io_mm + use stdlib_kinds, only : int8, int16, int32, int64, sp, dp, xdp, qp + implicit none + private + + type, public :: mm_header_type + integer :: object + integer :: format + integer :: qualifier + integer :: symmetry + character(len=1024), allocatable :: comments(:) + end type mm_header_type + + !> Version: experimental + !> + !> Load a matrix from a Matrix Market file + !> ([Specification](../page/specs/stdlib_io.html#load_mm)) + interface load_mm + #: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 + 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 to be read from the Matrix Market file + integer, allocatable, intent(out) :: index(:,:) + !> Matrix data to be read from the Matrix Market file + ${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 + end subroutine + #:endfor + end interface + public :: load_mm + + !> Version: experimental + !> + !> Save a matrix to a Matrix Market file + !> ([Specification](../page/specs/stdlib_io.html#save_mm)) + interface save_mm + #:for k, t, s in RCI_KINDS_TYPES + module subroutine save_mm_dense_${s}$(filename, matrix, comment, format, symmetry, iostat, iomsg) + character(len=*), intent(in) :: filename + ${t}$, intent(in) :: matrix(:,:) + character(len=*), intent(in), optional :: comment + character(len=*), intent(in), optional :: format + character(len=*), intent(in), optional :: symmetry + integer, intent(out), optional :: iostat + character(len=:), allocatable, intent(out), optional :: iomsg + end subroutine + #:endfor + + #:for k, t, s in RCI_KINDS_TYPES + module subroutine save_mm_coo_${s}$(filename, index, data, comment, format, symmetry, iostat, iomsg) + character(len=*), intent(in) :: filename + integer, intent(in) :: index(:,:) + ${t}$, intent(in) :: data(:) + character(len=*), intent(in), optional :: comment + character(len=*), intent(in), optional :: format + character(len=*), intent(in), optional :: symmetry + integer, intent(out), optional :: iostat + character(len=:), allocatable, intent(out), optional :: iomsg + end subroutine + #:endfor + end interface save_mm + public :: save_mm + +end module stdlib_io_mm \ No newline at end of file diff --git a/src/io/stdlib_io_mm_load.fypp b/src/io/stdlib_io_mm_load.fypp new file mode 100644 index 000000000..c81b13b51 --- /dev/null +++ b/src/io/stdlib_io_mm_load.fypp @@ -0,0 +1,483 @@ +! 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 Implementation for saving multidimensional arrays to Matrix Market files +submodule (stdlib_io_mm) stdlib_io_mm_save + use stdlib_error, only : error_stop + use stdlib_strings, only : to_string + use stdlib_io, only : open + use stdlib_ascii, only : to_lower + implicit none + + ! Matrix Market format constants + character(len=*), parameter :: MM_BANNER = "%%MatrixMarket" + character(len=*), parameter :: MM_COMMENT_CHAR = "%" + + ! Matrix Market object types + character(len=*), parameter :: & + MM_MATRIX = "matrix", & + MM_VECTOR = "vector" + + ! Matrix Market format types + character(len=*), parameter :: & + MM_COORDINATE = "coordinate", & + MM_ARRAY = "array" + + ! Matrix Market data types + character(len=*), parameter :: & + MM_REAL = "real", & + MM_COMPLEX = "complex", & + MM_INTEGER = "integer", & + MM_PATTERN = "pattern" + + ! Matrix Market storage schemes + character(len=*), parameter :: & + MM_GENERAL = "general", & + MM_SYMMETRIC = "symmetric", & + MM_SKEW_SYMMETRIC = "skew-symmetric", & + MM_HERMITIAN = "hermitian", & + MM_AUTO = "auto" + +contains + + #:for k, t, s in RCI_KINDS_TYPES + module subroutine save_mm_dense_${s}$(filename, matrix, comment, format, symmetry, iostat, iomsg) + !> Name of the Matrix Market file to save to + character(len=*), intent(in) :: filename + !> Matrix to be saved to the Matrix Market file + ${t}$, intent(in) :: matrix(:,:) + !> Optional comment information + character(len=*), intent(in), optional :: comment + !> Format in which matrix data needs to be stored + character(len=*), intent(in), optional :: format + !> Symmetry type of the matrix (general, symmetric, skew-symmetric, hermitian) + character(len=*), intent(in), optional :: symmetry + !> Error status of saving, 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 + + integer :: io, stat, i, j, nnz, nrows, ncols + character(len=:), allocatable :: msg + character(len=:), allocatable :: field_type + character(len=:), allocatable :: fmt_ + character(len=:), allocatable :: symmetry_ + #:if t.startswith('complex') + real(${k}$) :: real_part, imag_part + #:endif + if(present(iostat)) iostat = 0 + if(present(iomsg)) iomsg = '' + stat = 0 + + #:if t.startswith('integer') + fmt_ = "I0" + #:else + fmt_ = "ES24.16E3" + #:endif + if(present(format)) fmt_ = format + + #:if t.startswith('real') + fmt_ = '(' // fmt_ //')' + #:elif t.startswith('complex') + fmt_ = '(' // fmt_ //',1X,'// fmt_ //')' + #:elif t.startswith('integer') + fmt_ = '('// fmt_ //')' + #:endif + + io = open(filename, "w", iostat=stat) + if (stat /= 0) then + call mm_fail_process(iostat = iostat, iomsg = iomsg, code = stat,& + message = "Could not create file: " // filename) + return + end if + + nrows = size(matrix, 1) + ncols = size(matrix, 2) + ! Determine symmetry type + symmetry_ = MM_GENERAL + if (present(symmetry)) then + symmetry_ = to_lower(trim(symmetry)) + if (symmetry_ == MM_AUTO) then + + symmetry_block: block + ! Non-square matrices cannot be symmetric/skew/hermitian + if (nrows /= ncols) then + symmetry_ = MM_GENERAL + exit symmetry_block + end if + + ! Try symmetric + symmetry_ = MM_SYMMETRIC + do j = 1, ncols + do i = j+1, nrows + if (matrix(i,j) /= matrix(j,i)) then + symmetry_ = MM_GENERAL + exit + end if + end do + if (symmetry_ == MM_GENERAL) exit + end do + if (symmetry_ == MM_SYMMETRIC) exit symmetry_block + + ! Try skew-symmetric + symmetry_ = MM_SKEW_SYMMETRIC + do j = 1, ncols + do i = j, nrows + if (matrix(i,j) /= -matrix(j,i)) then + symmetry_ = MM_GENERAL + exit + end if + end do + if (symmetry_ == MM_GENERAL) exit + end do + if (symmetry_ == MM_SKEW_SYMMETRIC) exit symmetry_block + + #:if t.startswith('complex') + ! Try hermitian + symmetry_ = MM_HERMITIAN + do j = 1, ncols + do i = j, nrows + if (matrix(i,j) /= conjg(matrix(j,i))) then + symmetry_ = MM_GENERAL + exit + end if + end do + if (symmetry_ == MM_GENERAL) exit + end do + if (symmetry_ == MM_HERMITIAN) exit symmetry_block + #:endif + + symmetry_ = MM_GENERAL + end block symmetry_block + end if + end if + + ! Determine field type based on matrix type + #:if t.startswith('real') + field_type = MM_REAL + #:elif t.startswith('complex') + field_type = MM_COMPLEX + #:elif t.startswith('integer') + field_type = MM_INTEGER + #:endif + + catch: block + ! Write header + call write_mm_header(io, MM_ARRAY, field_type, symmetry_, & + nrows, ncols, comment=comment, iostat=stat, iomsg=msg) + if (stat /= 0) exit catch + + ! Write array format (column-major order) + if(symmetry_ == MM_GENERAL) then + do j = 1, ncols + do i = 1, nrows + #:if t.startswith('real') + write(io, fmt=fmt_, iostat=stat) matrix(i, j) + #:elif t.startswith('complex') + real_part = real(matrix(i, j), kind=${k}$) + imag_part = aimag(matrix(i, j)) + write(io, fmt=fmt_, iostat=stat) real_part, imag_part + #:elif t.startswith('integer') + write(io, fmt=fmt_, iostat=stat) matrix(i, j) + #:endif + if (stat /= 0) then + msg = "Error writing array element (" // & + to_string(i) // "," // to_string(j) // ")" + exit catch + end if + end do + end do + else + ! For symmetric and hermitian matrices, only the lower triangle + ! (including the diagonal) is written. + ! For skew-symmetric matrices, only the strictly lower triangle is written + ! (the diagonal is omitted and assumed zero). + do j = 1, ncols + do i = j, nrows + if(symmetry_ == MM_SKEW_SYMMETRIC .and. i == j) cycle + #:if t.startswith('real') + write(io, fmt=fmt_, iostat=stat) matrix(i, j) + #:elif t.startswith('complex') + real_part = real(matrix(i, j), kind=${k}$) + imag_part = aimag(matrix(i, j)) + if(i==j .and. symmetry_ == MM_HERMITIAN) imag_part = 0 + write(io, fmt=fmt_, iostat=stat) real_part, imag_part + #:elif t.startswith('integer') + write(io, fmt=fmt_, iostat=stat) matrix(i, j) + #:endif + if (stat /= 0) then + msg = "Error writing array element (" // & + to_string(i) // "," // to_string(j) // ")" + exit catch + end if + end do + end do + end if + end block catch + + close(io) + if(stat/=0) then + call mm_fail_process(iostat = iostat, iomsg = iomsg, code = stat,& + message = "Failed to save Matrix Market file '" // filename // "': " // msg) + return + end if + + if (present(iomsg) .and. allocated(msg)) call move_alloc(msg, iomsg) + end subroutine + #:endfor + + #:for k, t, s in RCI_KINDS_TYPES + module subroutine save_mm_coo_${s}$(filename, index, data, comment, format, symmetry, iostat, iomsg) + !> Name of the Matrix Market file to save to + character(len=*), intent(in) :: filename + !> Matrix indices in COO format to be written to the file: + !> index(1, :) = row indices, index(2, :) = column indices + integer, intent(in) :: index(:, :) + !> Nonzero matrix values corresponding to each (row, column) index pair + ${t}$, intent(in) :: data(:) + !> Optional comment information + character(len=*), intent(in), optional :: comment + !> Format in which matrix data needs to be stored + character(len=*), intent(in), optional :: format + !> Symmetry type of the matrix (general, symmetric, skew-symmetric, hermitian) + character(len=*), intent(in), optional :: symmetry + !> Error status of saving, 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 + + integer :: io, stat, i, nnz_to_write + character(len=:), allocatable :: msg + character(len=:), allocatable :: field_type + character(len=:), allocatable :: fmt_ + character(len=:), allocatable :: symmetry_ + logical :: is_pattern + #:if t.startswith('complex') + real(${k}$) :: real_part, imag_part + #:endif + if(present(iostat)) iostat = 0 + if(present(iomsg)) iomsg = '' + stat = 0 + + if(size(data)==1 .and. size(data)/=size(index,dim=2)) then + is_pattern = .true. + else + is_pattern = .false. + end if + + if(is_pattern) then + fmt_ = '(I0,1X,I0)' + else + #:if t.startswith('integer') + fmt_ = "I0" + #:else + fmt_ = "ES24.16E3" + #:endif + if(present(format)) fmt_ = format + #:if t.startswith('real') + fmt_ = '(I0,1X,I0,1X,' // fmt_ //')' + #:elif t.startswith('complex') + fmt_ = '(I0,1X,I0,1X,' // fmt_//',1X,'//fmt_//')' + #:elif t.startswith('integer') + fmt_ = '(I0,1X,I0,1X,'// fmt_ //')' + #:endif + end if + + io = open(filename, "w", iostat=stat) + if (stat /= 0) then + call mm_fail_process(iostat = iostat, iomsg = iomsg, code = stat,& + message = "Could not create file: " // filename) + return + end if + + if (.not. is_pattern .and. size(data) /= size(index,dim=2)) then + call mm_fail_process(iostat=iostat, iomsg=iomsg, code=1, & + message="Invalid COO data size") + return + end if + + if(size(index, dim=1)/=2) then + call mm_fail_process(iostat = iostat, iomsg = iomsg, code = stat,& + message = "Invalid index dimensions: first dimension must be 2") + return + end if + + ! Determine symmetry type + symmetry_ = MM_GENERAL + if (present(symmetry)) then + symmetry_ = to_lower(trim(symmetry)) + end if + + ! Determine field type based on matrix type + if(is_pattern) then + field_type = MM_PATTERN + else + #:if t.startswith('real') + field_type = MM_REAL + #:elif t.startswith('complex') + field_type = MM_COMPLEX + #:elif t.startswith('integer') + field_type = MM_INTEGER + #:endif + end if + + catch: block + ! Calculate the nnz to write inside mtx file + if (symmetry_ == MM_GENERAL) then + nnz_to_write = size(index, dim=2) + else if(symmetry_ == MM_SKEW_SYMMETRIC) then + nnz_to_write = count(index(1,:) > index(2,:)) + else + nnz_to_write = count(index(1,:) >= index(2,:)) + end if + ! Write header + if(nnz_to_write == 0) then + call write_mm_header(io, MM_COORDINATE, field_type, symmetry_, & + 0, 0,& + nnz=nnz_to_write, comment=comment, iostat=stat, iomsg=msg) + exit catch + end if + call write_mm_header(io, MM_COORDINATE, field_type, symmetry_, & + maxval(index(1,:)), maxval(index(2,:)),& + nnz=nnz_to_write, comment=comment, iostat=stat, iomsg=msg) + if (stat /= 0) exit catch + + ! Write coordinate format (row, column, value) + if(symmetry_ == MM_GENERAL) then + if(is_pattern) then + do i = 1, nnz_to_write + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i) + if (stat /= 0) then + msg = "Error writing array element (" // to_string(i) // ")" + exit catch + end if + end do + else + do i = 1, nnz_to_write + #:if t.startswith('real') + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i), data(i) + #:elif t.startswith('complex') + real_part = real(data(i), kind=${k}$) + imag_part = aimag(data(i)) + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i), real_part, imag_part + #:elif t.startswith('integer') + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i), data(i) + #:endif + if (stat /= 0) then + msg = "Error writing array element (" // to_string(i) // ")" + exit catch + end if + end do + end if + else + ! For symmetric and hermitian matrices, only the lower triangle + ! (including the diagonal) is written. + ! For skew-symmetric matrices, only the strictly lower triangle is written + ! (the diagonal is omitted and assumed zero). + if(is_pattern) then + do i = 1, size(index, dim=2) + if(index(1,i) < index(2,i)) cycle + if(symmetry_ == MM_SKEW_SYMMETRIC .and. index(1,i) == index(2,i)) cycle + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i) + if (stat /= 0) then + msg = "Error writing array element (" // to_string(i) // ")" + exit catch + end if + end do + else + do i = 1, size(index, dim=2) + if(index(1,i) < index(2,i)) cycle + if(symmetry_ == MM_SKEW_SYMMETRIC .and. index(1,i) == index(2,i)) cycle + #:if t.startswith('real') + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i), data(i) + #:elif t.startswith('complex') + real_part = real(data(i), kind=${k}$) + imag_part = aimag(data(i)) + if(index(1,i)==index(2,i) .and. symmetry_ == MM_HERMITIAN) imag_part = 0 + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i), real_part, imag_part + #:elif t.startswith('integer') + write(io, fmt=fmt_, iostat=stat) & + index(1,i), index(2,i), data(i) + #:endif + if (stat /= 0) then + msg = "Error writing array element (" // to_string(i) // ")" + exit catch + end if + end do + end if + end if + end block catch + + close(io) + + if(stat/=0) then + call mm_fail_process(iostat = iostat, iomsg = iomsg, code = stat,& + message = "Failed to save Matrix Market file '" // filename // "': " // msg) + return + end if + + if (present(iomsg) .and. allocated(msg)) call move_alloc(msg, iomsg) + end subroutine + #:endfor + + !> Write Matrix Market header + subroutine write_mm_header(io, format, field, symmetry, nrows, ncols, nnz, & + comment, iostat, iomsg) + integer, intent(in) :: io + character(len=*), intent(in) :: format, field, symmetry + integer, intent(in) :: nrows, ncols + integer, intent(in), optional :: nnz + character(len=*), intent(in), optional :: comment + integer, intent(out) :: iostat + character(len=:), allocatable, intent(out) :: iomsg + + integer :: stat + character(len=*), parameter :: iso_date_fmt = '(I4.4,"-",I2.2,"-",I2.2)' + integer :: date_values(8) + + iostat = 0 + + ! Write banner line + write(io, '(A)', iostat=stat) MM_BANNER // " " // MM_MATRIX // " " // & + format // " " // field // " " // symmetry + if (stat /= 0) then + iostat = stat + iomsg = "Error writing Matrix Market banner" + return + end if + + ! Write comments (including optional header_info and generation info) + call date_and_time(values=date_values) + write(io, '(A)', iostat=stat) "% Generated by Fortran stdlib on " // & + to_string(date_values(1)) // "-" // & + to_string(date_values(2)) // "-" // & + to_string(date_values(3)) + if (stat /= 0) then + iostat = stat + iomsg = "Error writing comment line" + return + end if + + if (present(comment)) then + if(len_trim(comment) > 0) then + write(io, '(A)', iostat=stat) "% " // trim(comment) + if (stat /= 0) then + iostat = stat + iomsg = "Error writing header info" + return + end if + end if + end if + + ! Write size line + if (format == MM_COORDINATE) then + write(io, '(I0,1X,I0,1X,I0)', iostat=stat) nrows, ncols, nnz + else + write(io, '(I0,1X,I0)', iostat=stat) nrows, ncols + end if + + if (stat /= 0) then + iostat = stat + iomsg = "Error writing matrix dimensions" + return + end if + end subroutine write_mm_header + + 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_save \ No newline at end of file diff --git a/test/data/matrix_market/README.md b/test/data/matrix_market/README.md new file mode 100644 index 000000000..511d9e309 --- /dev/null +++ b/test/data/matrix_market/README.md @@ -0,0 +1,76 @@ +# Test Data: Matrix Market + +This directory contains matrices used for testing purposes. These matrices are obtained from the SuiteSparse Matrix Collection (formerly the University of Florida Sparse Matrix Collection): + +https://sparse.tamu.edu + +## License + +The matrices in this directory are distributed under the CC-BY 4.0 license: + +https://creativecommons.org/licenses/by/4.0/ + +## Attribution + +- Kolodziej et al., (2019). The SuiteSparse Matrix Collection Website Interface. Journal of Open Source Software, 4(35), 1244. DOI: https://doi.org/10.21105/joss.01244 + +- Timothy A. Davis and Yifan Hu. 2011. The University of Florida Sparse Matrix Collection. ACM Transactions on Mathematical Software 38, 1, Article 1 (December 2011), 25 pages. DOI: https://doi.org/10.1145/2049662.2049663 + +## Matrix Market Metadata + +The matrices stored in the Matrix Market (`.mtx`) files include metadata in their headers, which may contain additional citations specific to individual matrices. + +These headers have been preserved and must not be removed. + +## Matrix Sources + +The matrix market files used are listed below: +- https://sparse.tamu.edu/HB/ash85 +- https://sparse.tamu.edu/HB/bcsstk01 +- https://sparse.tamu.edu/Mallya/lhr01 + +## About .npy files + +The `.npy` files in this directory are derived from the corresponding Matrix Market (`.mtx`) files using SciPy and NumPy. + +### Storage format + +For coordinate (COO) matrices: + +- Non-pattern matrices: + - `*_data.npy` contains the nonzero values + - `*_indices.npy` contains the indices in a flattened format: + [row_1, row_2, ..., row_n, col_1, col_2, ..., col_n] + +- Pattern matrices: + - `*_indices.npy` contains the indices in the same flattened format + +For array-type matrices: +- `*_data.npy` contains the matrix entries + +### Generation of .npy files + +The following Python snippet was used to generate the `.npy` files: + +```python +import numpy as np +from scipy.io import mmread + +FILE_NAME = "mm_file_name.mtx" + +loaded = mmread(FILE_NAME) +indices = np.concatenate((loaded.row, loaded.col)) # not applicable for array-type matrices + +np.save(FILE_NAME + "_indices.npy", indices) # not applicable for array-type matrices +np.save(FILE_NAME + "_data.npy", loaded.data) # not applicable for pattern matrices +``` + +## Notes + +- The matrices included are solely for testing purposes. +- The matrices in the Matrix Market (`.mtx`) files are included without modification from their original source. +- The stdlib project itself remains licensed under the MIT License. +- The CC-BY 4.0 license applies only to the Matrix Market (`.mtx`) files in this directory. +- No warranties are provided by the original authors. + +© The original authors of the Matrix Market (`.mtx`) files, as listed in the SuiteSparse Matrix Collection. \ No newline at end of file diff --git a/test/data/matrix_market/ash85.mtx b/test/data/matrix_market/ash85.mtx new file mode 100644 index 000000000..293634394 --- /dev/null +++ b/test/data/matrix_market/ash85.mtx @@ -0,0 +1,318 @@ +%%MatrixMarket matrix coordinate pattern symmetric +%------------------------------------------------------------------------------- +% UF Sparse Matrix Collection, Tim Davis +% http://www.cise.ufl.edu/research/sparse/matrices/HB/ash85 +% name: HB/ash85 +% [SYMMETRIC PATTERN OF NORMAL MATRIX OF HOLLAND SURVEY. ASHKENAZI, 1974] +% id: 11 +% date: 1974 +% author: V. Askenazi +% ed: A. Curtis, I. Duff, J. Reid +% fields: title A name id date author ed kind +% kind: least squares problem +%------------------------------------------------------------------------------- +85 85 304 +1 1 +2 1 +6 1 +7 1 +8 1 +2 2 +3 2 +8 2 +9 2 +10 2 +3 3 +4 3 +10 3 +4 4 +5 4 +10 4 +11 4 +12 4 +5 5 +12 5 +23 5 +6 6 +7 6 +13 6 +14 6 +45 6 +7 7 +8 7 +14 7 +15 7 +16 7 +8 8 +9 8 +16 8 +18 8 +9 9 +10 9 +11 9 +18 9 +19 9 +20 9 +10 10 +11 10 +11 11 +12 11 +20 11 +21 11 +22 11 +12 12 +22 12 +23 12 +13 13 +44 13 +45 13 +14 14 +15 14 +45 14 +46 14 +47 14 +15 15 +16 15 +17 15 +31 15 +47 15 +48 15 +16 16 +17 16 +18 16 +17 17 +18 17 +29 17 +30 17 +31 17 +18 18 +19 18 +28 18 +29 18 +19 19 +20 19 +27 19 +28 19 +20 20 +21 20 +27 20 +21 21 +22 21 +24 21 +25 21 +26 21 +27 21 +22 22 +23 22 +24 22 +23 23 +24 23 +32 23 +24 24 +25 24 +32 24 +33 24 +34 24 +25 25 +26 25 +34 25 +26 26 +27 26 +34 26 +36 26 +37 26 +27 27 +28 27 +37 27 +39 27 +28 28 +29 28 +39 28 +29 29 +30 29 +39 29 +40 29 +30 30 +31 30 +40 30 +42 30 +43 30 +31 31 +43 31 +48 31 +49 31 +32 32 +33 32 +64 32 +33 33 +34 33 +63 33 +64 33 +34 34 +35 34 +36 34 +64 34 +65 34 +35 35 +36 35 +38 35 +65 35 +66 35 +36 36 +37 36 +38 36 +37 37 +38 37 +39 37 +38 38 +39 38 +66 38 +67 38 +69 38 +39 39 +40 39 +41 39 +69 39 +70 39 +40 40 +41 40 +42 40 +71 40 +41 41 +70 41 +71 41 +42 42 +43 42 +50 42 +71 42 +72 42 +43 43 +49 43 +50 43 +51 43 +44 44 +45 44 +52 44 +53 44 +54 44 +45 45 +46 45 +52 45 +57 45 +46 46 +47 46 +57 46 +58 46 +47 47 +48 47 +58 47 +48 48 +49 48 +58 48 +59 48 +49 49 +51 49 +59 49 +60 49 +50 50 +51 50 +72 50 +73 50 +51 51 +60 51 +61 51 +73 51 +52 52 +54 52 +57 52 +81 52 +53 53 +54 53 +55 53 +85 53 +54 54 +55 54 +56 54 +81 54 +55 55 +56 55 +83 55 +84 55 +85 55 +56 56 +81 56 +82 56 +83 56 +57 57 +58 57 +80 57 +81 57 +58 58 +59 58 +80 58 +59 59 +60 59 +61 59 +79 59 +80 59 +60 60 +61 60 +61 61 +62 61 +73 61 +76 61 +78 61 +79 61 +62 62 +73 62 +74 62 +76 62 +78 62 +63 63 +64 63 +64 64 +65 64 +65 65 +66 65 +66 66 +67 66 +67 67 +68 67 +69 67 +68 68 +69 68 +69 69 +70 69 +70 70 +71 70 +71 71 +72 71 +72 72 +73 72 +73 73 +74 73 +74 74 +75 74 +76 74 +75 75 +76 75 +76 76 +77 76 +78 76 +77 77 +78 77 +78 78 +79 78 +79 79 +80 79 +80 80 +81 80 +81 81 +82 81 +82 82 +83 82 +83 83 +84 83 +84 84 +85 84 +85 85 diff --git a/test/data/matrix_market/ash85.mtx_indices.npy b/test/data/matrix_market/ash85.mtx_indices.npy new file mode 100644 index 000000000..2251eb6d5 Binary files /dev/null and b/test/data/matrix_market/ash85.mtx_indices.npy differ diff --git a/test/data/matrix_market/bcsstk01.mtx b/test/data/matrix_market/bcsstk01.mtx new file mode 100644 index 000000000..728073cb1 --- /dev/null +++ b/test/data/matrix_market/bcsstk01.mtx @@ -0,0 +1,238 @@ +%%MatrixMarket matrix coordinate real symmetric +%------------------------------------------------------------------------------- +% UF Sparse Matrix Collection, Tim Davis +% http://www.cise.ufl.edu/research/sparse/matrices/HB/bcsstk01 +% name: HB/bcsstk01 +% [SYMMETRIC STIFFNESS MATRIX SMALL GENERALIZED EIGENVALUE PROBLEM] +% id: 23 +% date: 1982 +% author: J. Lewis +% ed: I. Duff, R. Grimes, J. Lewis +% fields: title A name id date author ed kind +% kind: structural problem +%------------------------------------------------------------------------------- +48 48 224 +1 1 2832268.51852 +5 1 1e6 +6 1 2083333.33333 +7 1 -3333.33333333 +11 1 1e6 +19 1 -2.8e6 +25 1 -28935.1851852 +30 1 2083333.33333 +2 2 1635447.53086 +4 2 -2e6 +6 2 5555555.55555 +8 2 -6666.66666667 +10 2 -2e6 +20 2 -30864.1975309 +24 2 5555555.55555 +26 2 -1597916.66667 +3 3 1724367.28395 +4 3 -2083333.33333 +5 3 -2777777.77778 +9 3 -1.68e6 +21 3 -15432.0987654 +23 3 -2777777.77778 +27 3 -28935.1851852 +28 3 -2083333.33333 +4 4 1003333333.33 +8 4 2e6 +10 4 4e8 +22 4 -3333333.33333 +27 4 2083333.33333 +28 4 1e8 +5 5 1.0675e9 +7 5 -1e6 +11 5 2e8 +21 5 2777777.77778 +23 5 333333333.333 +29 5 -833333.333333 +6 6 1535333333.33 +12 6 -2e6 +20 6 -5555555.55555 +24 6 666666666.667 +25 6 -2083333.33333 +30 6 1e8 +7 7 2832268.51852 +11 7 -1e6 +12 7 2083333.33333 +13 7 -2.8e6 +31 7 -28935.1851852 +36 7 2083333.33333 +8 8 1635447.53086 +10 8 2e6 +12 8 5555555.55555 +14 8 -30864.1975309 +18 8 5555555.55555 +32 8 -1597916.66667 +9 9 1724367.28395 +10 9 -2083333.33333 +11 9 -2777777.77778 +15 9 -15432.0987654 +17 9 -2777777.77778 +33 9 -28935.1851852 +34 9 -2083333.33333 +10 10 1003333333.33 +16 10 -3333333.33333 +33 10 2083333.33333 +34 10 1e8 +11 11 1.0675e9 +15 11 2777777.77778 +17 11 333333333.333 +35 11 -833333.333333 +12 12 1535333333.33 +14 12 -5555555.55555 +18 12 666666666.667 +31 12 -2083333.33333 +36 12 1e8 +13 13 2836099.4695 +17 13 -2149285.29451 +18 13 2359161.80402 +19 13 -3333.33333333 +23 13 -1e6 +37 13 -28935.1851852 +42 13 2083333.33333 +43 13 -3830.95098171 +47 13 -1149285.29451 +48 13 275828.470683 +14 14 1767410.74446 +15 14 517922.131816 +16 14 4298570.58902 +18 14 -5555555.55555 +20 14 -6666.66666667 +22 14 2e6 +38 14 -1597916.66667 +44 14 -131963.213599 +45 14 -517922.131816 +46 14 2298570.58902 +15 15 3890038.06848 +16 15 -2634990.2747 +17 15 2777777.77778 +21 15 -1.68e6 +39 15 -28935.1851852 +40 15 -2083333.33333 +44 15 -517922.131816 +45 15 -2165670.78453 +46 15 -551656.941367 +16 16 1975720635.31 +20 16 -2e6 +22 16 4e8 +39 16 2083333.33333 +40 16 1e8 +44 16 -2298570.58902 +45 16 551656.941366 +46 16 486193650.99 +17 17 1527346515.47 +18 17 -109779731.332 +19 17 1e6 +23 17 2e8 +41 17 -833333.333333 +43 17 1149285.29451 +47 17 229724661.236 +48 17 -55717351.0779 +18 18 1564111437.11 +24 18 -2e6 +37 18 -2083333.33333 +42 18 1e8 +43 18 -275828.470683 +47 18 -55717351.0779 +48 18 10941196.0038 +19 19 2832268.51852 +23 19 1e6 +24 19 2083333.33333 +43 19 -28935.1851852 +48 19 2083333.33333 +20 20 1635447.53086 +22 20 -2e6 +24 20 -5555555.55555 +44 20 -1597916.66667 +21 21 1724367.28395 +22 21 -2083333.33333 +23 21 2777777.77778 +45 21 -28935.1851852 +46 21 -2083333.33333 +22 22 1003333333.33 +45 22 2083333.33333 +46 22 1e8 +23 23 1.0675e9 +47 23 -833333.333333 +24 24 1535333333.33 +43 24 -2083333.33333 +48 24 1e8 +25 25 60879.6296296 +29 25 1.25e6 +30 25 416666.666667 +31 25 -4166.66666667 +35 25 1.25e6 +26 26 3372916.66667 +28 26 -2.5e6 +32 26 -8333.33333333 +34 26 -2.5e6 +27 27 2411712.96296 +28 27 -416666.666667 +33 27 -2.355e6 +28 28 1.5e9 +32 28 2.5e6 +34 28 5e8 +29 29 501833333.333 +31 29 -1.25e6 +35 29 2.5e8 +30 30 5.025e8 +36 30 -2.5e6 +31 31 3985879.62963 +35 31 -1.25e6 +36 31 416666.666667 +37 31 -3.925e6 +32 32 3411496.91358 +34 32 2.5e6 +36 32 6944444.44444 +38 32 -38580.2469136 +42 32 6944444.44445 +33 33 2431003.08642 +34 33 -416666.666667 +35 33 -3472222.22222 +39 33 -19290.1234568 +41 33 -3472222.22222 +34 34 1504166666.67 +40 34 -4166666.66667 +35 35 1335166666.67 +39 35 3472222.22222 +41 35 416666666.667 +36 36 2169166666.67 +38 36 -6944444.44444 +42 36 833333333.333 +37 37 3985879.62963 +41 37 -1.25e6 +42 37 416666.666667 +43 37 -4166.66666667 +47 37 -1.25e6 +38 38 3411496.91358 +40 38 2.5e6 +42 38 -6944444.44445 +44 38 -8333.33333333 +46 38 2.5e6 +39 39 2431003.08642 +40 39 -416666.666667 +41 39 3472222.22222 +45 39 -2.355e6 +40 40 1504166666.67 +44 40 -2.5e6 +46 40 5e8 +41 41 1335166666.67 +43 41 1.25e6 +47 41 2.5e8 +42 42 2169166666.67 +48 42 -2.5e6 +43 43 64710.5806113 +47 43 2399285.29451 +48 43 140838.195984 +44 44 3504879.88027 +45 44 517922.131816 +46 44 -4798570.58902 +45 45 4577383.74749 +46 45 134990.2747 +46 46 2472387301.98 +47 47 961679848.804 +48 47 -109779731.332 +48 48 531278103.775 diff --git a/test/data/matrix_market/bcsstk01.mtx_data.npy b/test/data/matrix_market/bcsstk01.mtx_data.npy new file mode 100644 index 000000000..587b025e4 Binary files /dev/null and b/test/data/matrix_market/bcsstk01.mtx_data.npy differ diff --git a/test/data/matrix_market/bcsstk01.mtx_indices.npy b/test/data/matrix_market/bcsstk01.mtx_indices.npy new file mode 100644 index 000000000..f7f8268ff Binary files /dev/null and b/test/data/matrix_market/bcsstk01.mtx_indices.npy differ diff --git a/test/data/matrix_market/lhr01_b.mtx b/test/data/matrix_market/lhr01_b.mtx new file mode 100644 index 000000000..6f2859660 --- /dev/null +++ b/test/data/matrix_market/lhr01_b.mtx @@ -0,0 +1,1484 @@ +%%MatrixMarket matrix array real general +%------------------------------------------------------------------------------- +% UF Sparse Matrix Collection, Tim Davis +% http://www.cise.ufl.edu/research/sparse/matrices/Mallya/lhr01 +% name: Mallya/lhr01 : b matrix +%------------------------------------------------------------------------------- +1477 1 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +-3.18512309264e-7 +6.90905177952e-8 +.0350651168595 +-3.2326875953e-5 +-7.63939169701e-8 +1.02136254458e-9 +2.21098162001e-9 +3.33102434524e-10 +1.36139988172e-9 +5.09544406668e-10 +4.15741396864e-9 +3.73887587557e-10 +4.2811620915e-10 +6.47617071081e-10 +8.66094751473e-10 +2.16829221245e-10 +6.96871893524e-10 +9.02076635612e-10 +4.08221012549e-10 +6.15159478912e-10 +2.70767455674e-6 +0 +0 +0 +.0350651398076 +0 +6.26187102171e-8 +-2.52008981061 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +2.70767455674e-6 +0 +0 +-364.600326385 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +2.70767455674e-6 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +3.8224999595e-6 +0 +3.8224999595e-6 +0 +0 +-1.56319401867e-13 +20.4025355751 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +2.70767455674e-6 +0 +1.90024820768e-6 +324.78054912 +-2.04163370308e-7 +2.20984475163e-9 +7.71160557633e-9 +1.84983946383e-7 +7.34326022211e-9 +0 +-4.6554760047e-9 +-1.39516487252e-9 +-4.43014869234e-9 +-5.24221377418e-9 +-1.24538246382e-8 +-2.14183160097e-9 +-3.64997276847e-9 +-1.21076482174e-10 +-1.0762249758e-8 +-9.66767288446e-9 +-1.91752747014e-9 +-2.67269228971e-9 +-5.02788566337e-9 +-2.96603275274e-9 +-1.86297253522e-7 +0 +0 +0 +.0302915736684 +0 +-.000374576199368 +0 +.00458407276047 +0 +-.0102632853654 +0 +.0043951391887 +0 +.00393091191425 +0 +.0018186909557 +0 +.00120820496899 +0 +.000346859101932 +0 +.000228796488123 +0 +8.31426559778e-6 +0 +4.57459893725e-5 +0 +1.88149782718e-6 +0 +3.80871684787e-7 +0 +7.54859273006e-8 +0 +1.56610488198e-8 +0 +6.67079536493e-5 +0 +5.86948567891e-5 +0 +1.46027004981e-5 +0 +1.31468630813e-5 +0 +3.07391891621e-6 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +-1.86297253522e-7 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +-1.64188804774e-8 +0 +-1.64188804774e-8 +0 +0 +-1.56319401867e-13 +-268.228351898 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +-1.86299072511e-7 +0 +1.51572749019e-6 +265.182511195 +-1.06501829352e-8 +3.99249522332e-11 +1.36734001899e-9 +5.55110091227e-9 +1.27329258248e-9 +3.0786395655e-10 +4.63501237391e-10 +1.68256519828e-11 +-2.18619788939e-10 +-1.45746525959e-10 +-1.69598024513e-9 +-5.79575498705e-10 +-9.38626953939e-10 +-2.50105358646e-9 +-2.67314703706e-9 +-3.24840243593e-9 +-1.28210331241e-10 +-1.23918653117e-11 +-5.98191718382e-10 +-2.58893351202e-10 +4.10191205185e-7 +0 +0 +0 +-.00446066943379 +0 +.00151899737709 +0 +.00556025129786 +0 +-.117019994474 +0 +.00536197452727 +0 +.0036778339236 +0 +.0014990049338 +0 +.000972914431133 +0 +.000261602339309 +0 +.00017309635193 +0 +6.11115900213e-6 +0 +1.61370171795e-5 +0 +6.45993071336e-7 +0 +1.25830014052e-7 +0 +2.30023337414e-8 +0 +4.50507332985e-9 +0 +2.17875012349e-5 +0 +1.98706263134e-5 +0 +5.11218715431e-6 +0 +4.57586896816e-6 +0 +1.0650755504e-6 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +5.68434188608e-14 +0 +0 +0 +0 +0 +0 +4.10191205185e-7 +0 +-3.72616341338e-6 +.00446443757392 +3.26836016029e-8 +1.48525636234e-12 +-5.8164175698e-9 +-1.16710907605e-9 +-4.68241978289e-9 +-5.1579718247e-9 +-9.39849087011e-10 +3.40651240549e-9 +0 +1.09309894469e-9 +1.22179244499e-9 +8.83062512003e-11 +-1.92414972844e-11 +-2.05176320378e-10 +-3.21733750752e-11 +-1.84772375178e-9 +6.19877482677e-10 +8.47933279147e-10 +4.43975523012e-10 +2.54658516496e-10 +-6.45158689599e-8 +0 +-3.4514028942e-7 +0 +0 +0 +-.199238497502 +0 +-.000388266893176 +0 +-.796715267608 +0 +-.750142755376 +0 +-.409748449039 +0 +-.0595268675034 +0 +.0102318681196 +0 +.0138692563655 +0 +.0169840331029 +0 +.0160228238482 +0 +.00604078515611 +0 +.0105525683032 +0 +.00335846304061 +0 +.00177488283806 +0 +.000915876150734 +0 +.000476534789599 +0 +.0114864472596 +0 +.0114162932182 +0 +.00727897785613 +0 +.00696134830986 +0 +.00403577504429 +4.54747350886e-13 +0 +4.16193298331e-9 +-6.47446538265e-10 +2.81261236523e-10 +2.29874785873e-10 +-2.40731878876e-10 +-8.52310222399e-10 +2.84217094304e-14 +-8.97699692359e-11 +-1.14486198299e-10 +-5.17879072959e-11 +0 +2.22044604925e-16 +-1.33715261086e-12 +0 +-7.03153091308e-11 +-7.45004058444e-12 +-2.61994870243e-11 +-2.66524580184e-11 +3.24266683179e-8 +0 +3.00270807942e-8 +0 +6.70752342558e-8 +6.71774165508 +-4.45538717031e-9 +-9.56768518579e-12 +1.87660327322e-8 +-4.22300416858e-10 +-3.12255110657e-9 +3.48910588512e-9 +-1.6066792341e-9 +-6.46537046123e-10 +1.27897692437e-13 +-1.6567014427e-10 +-7.83089148731e-11 +-1.7763568394e-15 +-1.56193946665e-11 +-3.77761294202e-11 +-1.4157338496e-11 +-1.17973400146e-11 +-1.8065549046e-12 +-7.37010452667e-12 +-1.93698390564e-11 +-2.37663222435e-11 +-1.54411648082e-7 +0 +3.86219493183e-7 +0 +1.68980279265e-8 +0 +0 +0 +7.00383433023 +0 +1.14379308875e-5 +0 +-.0350545438584 +0 +-.0457023291204 +0 +-.0120731857962 +0 +.00682215509769 +0 +-.0216563098352 +0 +-.0298777781017 +0 +-.021474752021 +0 +-.0134831926787 +0 +-.000451448109188 +0 +.00236842175716 +0 +-.00133626318982 +0 +-.000497079969079 +0 +-.000140086761205 +0 +-3.78699910774e-5 +0 +.00266336244546 +0 +.00210335314361 +0 +.000417225107643 +0 +.00019883433047 +0 +-.00144218970873 +1.86255419582e-8 +1.51712148413e-12 +5.81344386213e-9 +5.05587388445e-9 +-2.48562059824e-10 +3.26542703578e-9 +1.46246747885e-9 +-2.49730192081e-9 +1.84741111298e-13 +-1.720152909e-10 +-4.74342787006e-12 +7.015277248e-12 +-7.62175738105e-12 +-4.1294155595e-12 +-7.53628550624e-13 +-2.49403713909e-13 +4.97104579722e-11 +-4.38848957174e-12 +5.41239275122e-12 +3.71802588717e-12 +2.23576882837e-7 +0 +-1.76339563041e-6 +0 +0 +0 +-.147643457071 +0 +2.91166479263e-5 +0 +-.0346500261325 +0 +-.0460476552 +0 +.0235049440078 +0 +.0108523870103 +0 +-.030249844952 +0 +-.0381378930894 +0 +-.0225889127188 +0 +-.0113848129854 +0 +-.00118129765668 +0 +.00605624943093 +0 +-.000182476234811 +0 +-2.66988381178e-5 +0 +-3.0854303241e-6 +0 +-3.53773532156e-7 +0 +.00859019572513 +0 +.00669164657482 +0 +.000294024951462 +0 +8.17557370022e-5 +0 +-.000270025770921 +1.07030473281e-8 +-8.51235748556e-13 +1.48542435922e-8 +-5.66444668948e-10 +-8.0077597886e-9 +1.13686837722e-13 +-8.52651282912e-14 +4.8577248832e-9 +-7.1054273576e-15 +-1.1932499433e-10 +-1.2122386428e-12 +1.46578305049e-11 +-5.58222088345e-13 +-3.95757577642e-14 +1.41288104372e-13 +7.09390093325e-21 +5.10900211026e-11 +-4.4408920985e-15 +9.40116039855e-13 +-3.916728053e-13 +7.80850991815e-8 +0 +5.65440475765e-7 +0 +0 +0 +-.176268530282 +0 +5.15047490787e-5 +0 +-.0180771100413 +0 +-.0389615009067 +0 +.0880599284593 +0 +.0188883626166 +0 +-.053803328564 +0 +-.0669397969019 +0 +-.0410674692195 +0 +-.020614586099 +0 +-.000285064934862 +0 +.0042732452155 +0 +-1.68867418714e-5 +0 +-1.00123522572e-6 +0 +-4.79837391018e-8 +0 +-2.34942099658e-9 +0 +.0118213585034 +0 +.00703422333679 +0 +7.35162426614e-5 +0 +1.16563684542e-5 +0 +-3.30451819915e-5 +1.03743786749e-8 +-1.64172495043e-13 +-3.69643871068e-10 +1.89913861369e-9 +-6.51068887692e-9 +3.37283267509e-9 +2.89659851515e-10 +-1.34065203383e-9 +2.9157121148e-10 +-2.7917224088e-11 +-1.51482992816e-13 +.00217874841656 +-2.01873294122e-14 +1.89050045614e-15 +-6.32222545006e-15 +3.8153705754e-23 +.00473591209159 +.00364647802045 +1.09482735276e-12 +5.74193470548e-14 +-.0333629605056 +0 +-.0223287145287 +0 +0 +0 +.220063585884 +0 +3.11446718569e-5 +0 +-.0128435925255 +0 +-.0275924955831 +0 +.0547109732548 +0 +.016938085487 +0 +-.0542819584328 +0 +-.0665829188399 +0 +-.0324581358704 +0 +-.00909817520297 +0 +2.26894837596e-5 +0 +-.00228410147734 +0 +-1.11561018202e-6 +0 +-2.89459728313e-8 +0 +-5.75874605189e-10 +0 +-1.18863165264e-11 +0 +-.00762509316763 +0 +-.00435957109097 +0 +9.63133133656e-5 +0 +6.24483009272e-5 +0 +-2.59826548116e-6 +-3.9794656459e-9 +-9.28175070894e-13 +1.44027180001e-8 +1.85267977003e-9 +-1.56694568432e-9 +-2.73189470892e-10 +1.13722364811e-10 +-3.68565622466e-10 +1.33618449575e-10 +5.93867177656e-11 +2.93044668003e-14 +-.0028507329893 +2.54053255045e-15 +-4.88706735198e-15 +1.74224916519e-23 +-2.42741966867e-25 +-.00778534876956 +-.00525394932293 +-7.56863664833e-6 +-4.55957798127e-6 +-.00621782676966 +0 +-.0115453676665 +0 +0 +0 +.0927779772604 +0 +-2.71089192304e-7 +0 +-.0448312690439 +0 +-.0464251064025 +0 +-.0245759551802 +0 +.0165876576612 +0 +-.00755238721016 +0 +-.0135251754058 +0 +-.00136469791203 +0 +.00778301173358 +0 +7.92439311327e-6 +0 +.000589716679691 +0 +-2.99774918264e-8 +0 +-3.4428484455e-10 +0 +-2.80354645834e-12 +0 +-2.39719817631e-14 +0 +.00184430120023 +0 +.00122663612949 +0 +-1.53813462076e-5 +0 +-9.47748499332e-6 +0 +-8.05977930824e-8 +-1.25965016196e-9 +1.44428900823e-12 +-5.98512883698e-9 +1.01513236952e-8 +-1.13686837722e-13 +3.55839802069e-11 +-4.1467274059e-11 +0 +3.5704772455e-12 +0 +8.86595484528e-15 +-3.63657515772e-5 +4.23516473627e-22 +6.61744490042e-24 +-1.03397576569e-25 +0 +-.000102773220837 +-5.78928589156e-5 +-3.06808280978e-7 +-1.03723500902e-7 +-.0055056092155 +0 +0 +0 +.0146945375861 +0 +-6.21794366383e-7 +0 +-.0180767901915 +0 +-.0182759162447 +0 +-.00985069030895 +0 +.0193188175897 +0 +-.00546646046677 +0 +-.0122040778633 +0 +-.00282610806789 +0 +.00606251002741 +0 +1.04985594005e-6 +0 +.000134681649983 +0 +-3.71874192744e-9 +0 +-1.6970293946e-11 +0 +-5.6232710224e-14 +0 +-1.9774396841e-16 +0 +.000557044192498 +0 +.000355597557505 +0 +2.6092196329e-6 +0 +1.89088380673e-6 +0 +-2.47110719026e-8 +-4.45470504928e-9 +1.6095458271e-13 +-1.23578303146e-9 +1.26672716249e-9 +-3.46972228726e-10 +3.4890490417e-9 +-1.13686837722e-13 +-3.23257154378e-9 +5.81508174416e-10 +-3.87444742955e-10 +-2.76486388798e-10 +2.16289208765e-10 +-1.84229520528e-10 +-1.74452452484e-10 +1.19467799711e-9 +-8.58335624798e-12 +-1.62714286159e-10 +-8.98978668312e-11 +-4.71430890391e-10 +2.71342059932e-10 +-2.53838150252e-7 +0 +6.63931132294e-8 +3.08655374273 +-1.07393134251e-8 +6.19643224789e-14 +1.2219842874e-9 +-1.98471639123e-9 +3.57772478165e-10 +1.85991666513e-9 +5.32850208401e-10 +3.40355654771e-9 +2.00975591724e-9 +3.08921244141e-9 +1.71964984474e-8 +-1.47340983858e-9 +1.16119736049e-9 +8.83858319867e-10 +4.19402113039e-9 +7.91663978816e-9 +2.7088447041e-9 +-8.53646042742e-10 +2.30386376643e-10 +3.02253511109e-9 +1.28198839344e-7 +0 +-2.05967350681e-6 +0 +0 +0 +-4.98707650155 +0 +-3.5988931178e-5 +0 +.140754465959 +0 +.110811093732 +0 +-.0912597597893 +0 +-.0312612509738 +0 +.00571147787014 +0 +.0148550322744 +0 +.0173726431684 +0 +.015166781601 +0 +.00420368145115 +0 +.00680235078504 +0 +.002746162188 +0 +.00171368548263 +0 +.0010213156521 +0 +.000607732014185 +0 +.00761312672467 +0 +.00766704833551 +0 +.00483310265065 +0 +.00468610062209 +0 +.00303533561571 +-3.41970007438e-10 +-1.16507844896e-12 +-4.88419971134e-9 +1.36239464155e-10 +-1.53477230924e-12 +-2.84217094304e-11 +3.71983333025e-10 +-1.79829839908e-9 +-9.63495949691e-10 +4.35363745055e-9 +7.61281171435e-9 +2.89048784907e-10 +-2.28462226914e-9 +2.04522621061e-10 +4.37239577877e-10 +-1.72173031387e-9 +-1.79369408215e-10 +3.87672116631e-10 +1.63532831721e-9 +-3.9797782847e-9 +-2.08933570541e-6 +0 +9.675864233e-5 +0 +0 +0 +-4.1491748948 +0 +2.7328070284e-5 +0 +.246459733582 +0 +.420957929396 +0 +.3930336947 +0 +.0539735410346 +0 +.019058927647 +0 +.0254248603778 +0 +.0162921083797 +0 +.0115318656192 +0 +-.00135112589203 +0 +-.00128190319611 +0 +-.000324334488554 +0 +6.70758396411e-5 +0 +.000154085448913 +0 +.000153192566331 +0 +-.000620918863371 +0 +-.000560434909542 +0 +-.00144104008013 +0 +-.00132151998453 +0 +-.000711121905939 +4.86652424572e-8 +-2.06456379635e-13 +-4.34300061959e-9 +-6.70915767387e-10 +6.85759005137e-10 +2.90515345114e-9 +-6.32326191408e-10 +-1.63504410011e-9 +1.40357769851e-9 +-5.51642642677e-9 +-2.06536014957e-8 +2.11682049667e-9 +-1.52203938342e-9 +8.17180989543e-10 +-5.93732352172e-9 +1.55677071234e-9 +-2.3580355446e-9 +3.4860931919e-9 +-8.16868350739e-10 +8.29714963402e-10 +1.02081187132e-6 +0 +-.000364611633729 +0 +0 +0 +-6.30642090981 +0 +2.64008013241e-5 +0 +.30189530113 +0 +.496790342051 +0 +.749877612579 +0 +.182834101855 +0 +.0404650727786 +0 +.0389542791709 +0 +.0130570226493 +0 +.00564082569494 +0 +-.00842503754363 +0 +-.0088888204679 +0 +-.00365670228084 +0 +-.00183467960423 +0 +-.000913900768568 +0 +-.000445487785651 +0 +-.00815518803669 +0 +-.00816858497938 +0 +-.00769180849061 +0 +-.00733325794912 +0 +-.0046707888683 +-2.02726369025e-8 +6.46050486134e-13 +-1.2742589206e-9 +1.47515777092e-9 +-2.39995756601e-9 +-4.59863258584e-9 +-3.84829945688e-10 +-1.5768364392e-9 +-2.23735696636e-9 +-2.96040525427e-10 +5.36147126695e-10 +-7.12276460035e-10 +1.38368250191e-9 +-6.4218852458e-10 +2.00122940441e-9 +-4.08718392464e-9 +1.4647127955e-9 +-3.30450689034e-9 +-2.9333762086e-9 +1.83430870493e-9 +-3.59836979485e-5 +0 +-.000323218017758 +0 +0 +0 +-7.63030547062 +0 +1.20245744923e-5 +0 +.282013791978 +0 +.349248292016 +0 +.79072318683 +0 +.206664937022 +0 +.0388690388754 +0 +.0380231814981 +0 +.0132592477786 +0 +.00573256118972 +0 +-.00869434340951 +0 +-.00974392260966 +0 +-.00394837318776 +0 +-.00199253650574 +0 +-.000968145205227 +0 +-.000471701055732 +0 +-.00884162881156 +0 +-.00883002323106 +0 +-.00845553256716 +0 +-.00798867992571 +0 +-.00508592823672 +-1.65609890246e-8 +2.76635029352e-12 +2.95635159657e-9 +-4.41751524249e-10 +2.65600874627e-10 +-9.41327015305e-11 +-5.92194737692e-10 +-9.18589648791e-10 +6.85531631461e-11 +-2.64549271378e-9 +-1.93938376469e-9 +-1.16227738545e-9 +6.98349822414e-10 +1.30634703055e-9 +-2.45708520197e-9 +-3.75391095986e-9 +-8.43130010253e-10 +-4.17003320763e-10 +9.01820840227e-10 +-2.94002688861e-9 +-6.44761098537e-5 +0 +.0226512106376 +0 +0 +0 +-6.28553778981 +0 +1.87262827964e-6 +0 +.118424213469 +0 +.0970308455209 +0 +.249771478472 +0 +-.0130361977898 +0 +-.0136352424274 +0 +.00403448455496 +0 +.0174502037485 +0 +.0157140902669 +0 +.00532388780604 +0 +.00595127967624 +0 +.00385768965776 +0 +.0027935724324 +0 +.00188108025294 +0 +.00121597445374 +0 +.00618547985039 +0 +.0070603404165 +0 +.00466800032895 +0 +.00513297627261 +0 +.00357462607304 +1.54359440785e-8 +-6.4061928505e-14 +1.2051799532e-9 +1.42108547152e-14 +-2.66481947619e-10 +3.41060513165e-13 +4.6304649004e-10 +-3.16049408866e-11 +-2.81033862848e-10 +-2.14413375943e-10 +-7.90691956354e-9 +-6.01090732744e-10 +-2.11883843804e-9 +-5.32872945769e-9 +-2.87516854769e-9 +-3.86935994356e-9 +-1.00351371657e-9 +-8.90992168934e-10 +-1.76865455614e-9 +-3.93356458517e-10 +.0812075537354 +0 +0 +0 +-.176933364035 +0 +6.77018307841e-7 +0 +.0264257505448 +0 +.0211922949515 +0 +.0562170636732 +0 +-.0508116234836 +0 +-.0126279574678 +0 +.00291646273862 +0 +.0248782142522 +0 +.0159386900522 +0 +.0118214849555 +0 +.00521387266521 +0 +.00284095476618 +0 +.00228060259598 +0 +.00146836355507 +0 +.000773694130742 +0 +.00409857200446 +0 +.00552844123009 +0 +.00285630982869 +0 +.00434207704071 +0 +.00292300739687 +4.03179001296e-9 +0 +-7.53210827042e-11 +-8.10480571545e-11 +-2.47212028626e-10 +-2.04295247386e-10 +9.47466105572e-10 +1.66460267792e-9 +1.36139988172e-9 +5.09544406668e-10 +7.48332240619e-9 +3.73859165848e-10 +4.28059365731e-10 +-2.15862883038e-10 +8.66094751473e-10 +-6.50430820058e-10 +-6.96815048891e-10 +1.62364699463e-9 +4.08192590839e-10 +1.43526790453e-9 +.140260541331 +.0350651398076 +2.03287907341e-20 +4.51954029757e-11 +5.40523181769e-12 +-2.90094961243e-11 +1.56251936421e-12 +-6.05431700421e-12 +2.9002716262e-15 +-3.84868123994e-12 +4.76719305862e-12 +4.89981359736e-12 +8.89872897962e-12 +-1.43646653834e-11 +-4.8109273975e-12 +1.92684747871e-11 +0 +4.27260505669e-12 +4.36321837338e-12 +-1.40071935177e-11 +-4.68085562293e-12 +4.8882532447e-12 +.00101508435232 +0 +2.56709111127 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +.0350651398076 +0 +0 +-2.06786550942 +-4.79310938158e-7 +-4.57907786184e-9 +-1.54977894694e-9 +7.71251507103e-8 +3.03771230392e-8 +3.08318703901e-10 +1.79511516762e-10 +3.18323145621e-12 +-7.03437308402e-13 +-7.3399064604e-12 +-1.67466041034e-12 +-1.05360165037e-12 +-3.97737398572e-14 +-1.28213412109e-14 +-3.87927537315e-16 +-6.22657307658e-16 +-4.88942220045e-13 +-5.40012479178e-13 +-1.88737914186e-13 +-3.02480263059e-13 +-3.16280113845e-6 +0 +0 +49.2012580097 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +-3.16280113845e-6 +0 +0 +-1.56319401867e-13 +-33.8844050154 diff --git a/test/data/matrix_market/lhr01_b.mtx_data.npy b/test/data/matrix_market/lhr01_b.mtx_data.npy new file mode 100644 index 000000000..b4f5eb652 Binary files /dev/null and b/test/data/matrix_market/lhr01_b.mtx_data.npy differ diff --git a/test/io/CMakeLists.txt b/test/io/CMakeLists.txt index 4e19b5fbe..f6f78ca82 100644 --- a/test/io/CMakeLists.txt +++ b/test/io/CMakeLists.txt @@ -2,6 +2,7 @@ set( fppFiles "test_loadtxt_qp.fypp" "test_savetxt_qp.fypp" + "test_io_mm.fypp" ) fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) @@ -14,6 +15,7 @@ set_tests_properties(loadtxt_qp PROPERTIES LABELS quadruple_precision) set_tests_properties(savetxt_qp PROPERTIES LABELS quadruple_precision) ADDTEST(get_line) +ADDTEST(io_mm) ADDTEST(npy) ADDTEST(open) ADDTEST(parse_mode) diff --git a/test/io/test_io_mm.fypp b/test/io/test_io_mm.fypp new file mode 100644 index 000000000..c111f9388 --- /dev/null +++ b/test/io/test_io_mm.fypp @@ -0,0 +1,597 @@ +#: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 KINDS_TYPES = R_KINDS_TYPES + C_KINDS_TYPES + I_KINDS_TYPES +module test_io_mm + use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test + use stdlib_kinds + use stdlib_math, only: all_close + use stdlib_io_npy, only: load_npy + use stdlib_io_mm + implicit none + + integer, parameter :: MS_general = 1 + integer, parameter :: MS_symmetric = 2 + integer, parameter :: MS_skew_symmetric = 3 + integer, parameter :: MS_hermitian = 4 + +contains + + + !> Collect all exported unit tests + subroutine collect_suite(testsuite) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest('io_mm_array', test_io_mm_array), & + new_unittest('io_mm_coordinate', test_io_mm_coordinate) & + ] + end subroutine + + #:for k, t, s in KINDS_TYPES + subroutine generate_random_${s}$_dense_matrix(A) + ${t}$, intent(out) :: A(:, :) + + ! Internal variables + #:if t.startswith('complex') + real(${k}$), allocatable :: R(:, :, :) + #:endif + #:if t.startswith('integer') + real :: rnd + integer :: i, j + #:endif + integer :: n + + n = size(A,dim=1) + + #:if t.startswith('real') + call random_number(A) + #:endif + #:if t.startswith('complex') + allocate(R(n,n,2)) + call random_number(R(:,:,1)) + call random_number(R(:,:,2)) + A = cmplx(R(:,:,1), R(:,:,2), kind=${k}$) + #:endif + #:if t.startswith('integer') + do j = 1, n + do i = 1,n + call random_number(rnd) + A(i,j) = int(rnd * 100 - 50, kind=${k}$) + end do + end do + #:endif + end subroutine + pure function compare_dense_${s}$(A, B) result(result) + ${t}$, intent(in) :: A(:, :), B(:, :) + logical :: result + + #:if t.startswith('integer') + result = all(A==B) + #:else + result = all_close(A, B) + #:endif + end function + #:endfor + + subroutine test_io_mm_array(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + integer :: u, stat + #:for k, t, s in (KINDS_TYPES) + block + integer, parameter :: n = 5 + ${t}$, allocatable :: matrix_save(:, :), matrix_load(:, :), A(:, :) + logical :: result + allocate(matrix_save(n,n)) + allocate(A(n,n)) + + call random_seed() + ! General matrix + call generate_random_${s}$_dense_matrix(matrix_save) + call save_mm("test_mmio_dense.mtx", matrix_save, format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=general, symmetry_arg=unspecified, type=${t}$") + if(allocated(error)) return + ! Check if symmetry = auto + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "auto", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=general, symmetry_arg=auto, type=${t}$") + if(allocated(error)) return + + ! Symmetric matrix + call generate_random_${s}$_dense_matrix(A) + ! Construct symmetric matrix using (A + A.T) + matrix_save = A + transpose(A) + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "symmetric", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=symmetric, symmetry_arg=symmetric, type=${t}$") + if(allocated(error)) return + ! Check if symmetry = auto + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "auto", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=symmetric, symmetry_arg=auto, type=${t}$") + if(allocated(error)) return + + ! Skew-symmetric matrix + call generate_random_${s}$_dense_matrix(A) + ! Construct symmetric matrix using (A - A.T) + matrix_save = A - transpose(A) + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "skew-symmetric", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=skew-symmetric, symmetry_arg=skew-symmetric, type=${t}$") + if(allocated(error)) return + ! Check if symmetry = auto + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "auto", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=skew-symmetric, symmetry_arg=auto, type=${t}$") + if(allocated(error)) return + + #:if t.startswith('complex') + ! Hermitian matrix + call generate_random_${s}$_dense_matrix(A) + ! Construct symmetric matrix using (A + A.H) + matrix_save = A + transpose(conjg(A)) + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "hermitian", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=hermitian, symmetry_arg=hermitian, type=${t}$") + if(allocated(error)) return + ! Check if symmetry = auto + call save_mm("test_mmio_dense.mtx", matrix_save, symmetry = "auto", format = "G0") + call load_mm("test_mmio_dense.mtx", matrix_load) + result = compare_dense_${s}$(matrix_save, matrix_load) + call check(error, result, .true.,& + "MM array test failed: matrix=hermitian, symmetry_arg=auto, type=${t}$") + if(allocated(error)) return + #:endif + end block + #:endfor + + stat = 0 + open(newunit=u, file="test_mmio_dense.mtx", status="old", iostat=stat) + if (stat == 0) then + close(u, status="delete") + end if + end subroutine + + subroutine generate_random_positions(pos, t_entries) + integer, intent(out) :: pos(:) + integer, intent(in) :: t_entries + + ! Internal variables + integer :: i, j, temp + real(dp) :: rnd + + pos = [(i,i=1,t_entries)] + do i = t_entries, 2, -1 + call random_number(rnd) + j = int(i*rnd) + 1 + temp = pos(i) + pos(i) = pos(j) + pos(j) = temp + end do + end subroutine + + subroutine fill_first_half_indices(index_save, pos, nnz, nrows, ncols, symmetry, j) + integer, intent(out) :: index_save(:, :) + integer, intent(in) :: pos(:), nnz, nrows, ncols, symmetry + integer, intent(out) :: j + + ! Internal variables + integer :: i, row, col + + if(symmetry == MS_symmetric .or. symmetry == MS_hermitian) then + j = 1 + do i = 1, nnz + row = mod(pos(i) - 1,nrows) + 1 + col = (pos(i) - 1)/nrows + 1 + if(row < col) cycle + index_save(1,j) = row + index_save(2,j) = col + j = j + 1 + end do + else + j = 1 + do i = 1, nnz + row = mod(pos(i) - 1,nrows) + 1 + col = (pos(i) - 1)/nrows + 1 + if(row <= col) cycle + index_save(1,j) = row + index_save(2,j) = col + j = j + 1 + end do + end if + end subroutine + + #:for k, t, s in KINDS_TYPES + subroutine generate_random_data_for_${s}$_coo(A, nnz_to_write) + ${t}$, intent(out) :: A(:) + integer, intent(in) :: nnz_to_write + + ! Internal variables + #:if t.startswith('complex') + real(${k}$), allocatable :: R(:, :) + #:endif + #:if t.startswith('integer') + real :: rnd + integer :: i + #:endif + + #:if t.startswith('real') + call random_number(A(1:nnz_to_write)) + #:endif + #:if t.startswith('complex') + allocate(R(nnz_to_write, 2)) + call random_number(R(:,1)) + call random_number(R(:,2)) + A(1:nnz_to_write) = cmplx(R(1:nnz_to_write, 1), R(1:nnz_to_write, 2), kind=${k}$) + #:endif + #:if t.startswith('integer') + do i = 1, nnz_to_write + call random_number(rnd) + A(i) = int(rnd * 100 - 50, kind=${k}$) + end do + #:endif + end subroutine + + pure function compare_coo_${s}$(index_save, index_load, data_save, data_load) result(result) + ${t}$, intent(in) :: data_save(:), data_load(:) + integer, intent(in) :: index_save(:, :), index_load(:,:) + logical :: result + + #:if t.startswith('integer') + result = all(index_save == index_load) .and. all(data_save == data_load) + #:else + result = all(index_save == index_load) .and. all_close(data_save, data_load) + #:endif + end function + + subroutine fill_other_half_${s}$(index_save, data_save, j, half_nnz, symmetry, is_pattern) + integer, intent(inout) :: index_save(:, :) + ${t}$, intent(inout) :: data_save(:) + integer, intent(in) :: half_nnz, symmetry + integer, intent(inout) :: j + logical, intent(in) :: is_pattern + + ! Internal variables. + integer :: i + + if(.not. is_pattern) then + if(symmetry == MS_symmetric) then + do i = 1, half_nnz + if(index_save(1,i) == index_save(2,i)) cycle + index_save(1,j) = index_save(2,i) + index_save(2,j) = index_save(1,i) + data_save(j) = data_save(i) + j=j+1 + end do + #:if t.startswith('complex') + else if(symmetry == MS_hermitian) then + do i = 1, half_nnz + if(index_save(1,i) == index_save(2,i)) cycle + index_save(1,j) = index_save(2,i) + index_save(2,j) = index_save(1,i) + data_save(j) = conjg(data_save(i)) + j=j+1 + end do + #:endif + else + do i = 1, half_nnz + if(index_save(1,i) == index_save(2,i)) cycle + index_save(1,j) = index_save(2,i) + index_save(2,j) = index_save(1,i) + data_save(j) = -data_save(i) + j=j+1 + end do + end if + else + if(symmetry == MS_symmetric) then + do i = 1, half_nnz + if(index_save(1,i) == index_save(2,i)) cycle + index_save(1,j) = index_save(2,i) + index_save(2,j) = index_save(1,i) + j=j+1 + end do + #:if t.startswith('complex') + else if(symmetry == MS_hermitian) then + do i = 1, half_nnz + if(index_save(1,i) == index_save(2,i)) cycle + index_save(1,j) = index_save(2,i) + index_save(2,j) = index_save(1,i) + j=j+1 + end do + #:endif + else + do i = 1, half_nnz + if(index_save(1,i) == index_save(2,i)) cycle + index_save(1,j) = index_save(2,i) + index_save(2,j) = index_save(1,i) + j=j+1 + end do + end if + end if + end subroutine + + subroutine generate_random_${s}$_coo_matrix(index_save, data_save, nrows, ncols, symmetry, is_pattern) + ${t}$, allocatable, intent(out) :: data_save(:) + integer, allocatable, intent(out) :: index_save(:, :) + integer, intent(in) :: nrows, ncols, symmetry + logical, intent(in) :: is_pattern + + ! Internal variables + integer, allocatable :: pos(:) + integer :: nnz, nnz_lower, nnz_diag, i, j + #:if t.startswith('integer') + real(dp) :: density + #:else + real(${k}$) :: density + #:endif + + allocate(pos(nrows * ncols)) + call generate_random_positions(pos, nrows * ncols) + call random_number(density) + nnz = ceiling(density*nrows*ncols) + + if(symmetry == MS_general) then + allocate(index_save(2, nnz)) + if(.not. is_pattern) allocate(data_save(nnz)) + do i = 1, nnz + index_save(1,i) = mod(pos(i) - 1,nrows) + 1 + index_save(2,i) = (pos(i) - 1)/nrows + 1 + end do + if(.not. is_pattern) call generate_random_data_for_${s}$_coo(data_save, nnz) + else if(symmetry == MS_symmetric) then + nnz_lower = count(mod(pos(1:nnz) - 1,nrows) > (pos(1:nnz) - 1)/nrows) !! lower triangular part + nnz_diag = count(mod(pos(1:nnz) - 1,nrows) == (pos(1:nnz) - 1)/nrows) !! diagonal + allocate(index_save(2, 2*nnz_lower + nnz_diag)) + if(.not. is_pattern) allocate(data_save(2*nnz_lower + nnz_diag)) + call fill_first_half_indices(index_save, pos, nnz, nrows, ncols, MS_symmetric, j) + if(.not. is_pattern) call generate_random_data_for_${s}$_coo(data_save, nnz_lower + nnz_diag) + call fill_other_half_${s}$(index_save, data_save, j, nnz_lower+nnz_diag, MS_symmetric, is_pattern) + #:if t.startswith('complex') + else if(symmetry == MS_hermitian) then + nnz_lower = count(mod(pos(1:nnz) - 1,nrows) > (pos(1:nnz) - 1)/nrows) !! lower triangular part + nnz_diag = count(mod(pos(1:nnz) - 1,nrows) == (pos(1:nnz) - 1)/nrows) !! diagonal + allocate(index_save(2, 2*nnz_lower + nnz_diag)) + if(.not. is_pattern) allocate(data_save(2*nnz_lower + nnz_diag)) + call fill_first_half_indices(index_save, pos, nnz, nrows, ncols, MS_hermitian, j) + if(.not. is_pattern) call generate_random_data_for_${s}$_coo(data_save, nnz_lower+nnz_diag) + if(.not. is_pattern) then + do i = 1, nnz_lower + nnz_diag + if(index_save(1, i) == index_save(2,i)) data_save(i) = real(data_save(i)) + end do + end if + call fill_other_half_${s}$(index_save, data_save, j, nnz_lower+nnz_diag, MS_hermitian, is_pattern) + #:endif + else + nnz_lower = count(mod(pos(1:nnz) - 1,nrows) > (pos(1:nnz) - 1)/nrows) !! lower triangular part + allocate(index_save(2, 2*nnz_lower)) + if(.not. is_pattern) allocate(data_save(2*nnz_lower)) + call fill_first_half_indices(index_save, pos, nnz, nrows, ncols, MS_skew_symmetric, j) + if(.not. is_pattern) call generate_random_data_for_${s}$_coo(data_save, nnz_lower) + call fill_other_half_${s}$(index_save, data_save, j, nnz_lower, MS_skew_symmetric, is_pattern) + end if + if(allocated(pos)) deallocate(pos) + end subroutine + #:endfor + + pure function compare_coo_pattern(index_save, index_load) result(result) + integer, intent(in) :: index_save(:, :), index_load(:,:) + logical :: result + + result = all(index_save == index_load) + end function + + subroutine test_io_mm_coordinate(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + integer :: u, stat + #:for k, t, s in KINDS_TYPES + block + integer :: nrows, ncols + ${t}$, allocatable :: data_save(:), data_load(:) + integer, allocatable :: index_save(:, :), index_load(:,:) + logical :: result + + nrows = 5 + ncols = 5 + + call random_seed() + ! General matrix + call generate_random_${s}$_coo_matrix(index_save, data_save, nrows, ncols, MS_general, .false.) + call save_mm("test_mmio_sparse.mtx", index_save, data_save, format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, data_load) + result = compare_coo_${s}$(index_save, index_load, data_save, data_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=unspecified, type=${t}$") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + if(allocated(data_save)) deallocate(data_save) + + ! Symmetric matrix + call generate_random_${s}$_coo_matrix(index_save, data_save, nrows, ncols, MS_symmetric, .false.) + call save_mm("test_mmio_sparse.mtx", index_save, data_save, symmetry = "symmetric", format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, data_load) + result = compare_coo_${s}$(index_save, index_load, data_save, data_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=symmetric, type=${t}$") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + if(allocated(data_save)) deallocate(data_save) + + ! Skew-symmetric matrix + call generate_random_${s}$_coo_matrix(index_save, data_save, nrows, ncols, MS_skew_symmetric, .false.) + call save_mm("test_mmio_sparse.mtx", index_save, data_save, symmetry = "skew-symmetric", format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, data_load) + result = compare_coo_${s}$(index_save, index_load, data_save, data_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=skew-symmetric, type=${t}$") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + if(allocated(data_save)) deallocate(data_save) + + #:if t.startswith('complex') + ! Hermitian matrix + call generate_random_${s}$_coo_matrix(index_save, data_save, nrows, ncols, MS_hermitian, .false.) + call save_mm("test_mmio_sparse.mtx", index_save, data_save, symmetry = "hermitian", format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, data_load) + result = compare_coo_${s}$(index_save, index_load, data_save, data_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=hermitian, type=${t}$") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + if(allocated(data_save)) deallocate(data_save) + #:endif + end block + #:endfor + + ! Pattern tests + block + integer :: nrows, ncols + real(sp), allocatable :: dummy(:) ! Dummy data matrix + integer, allocatable :: index_save(:, :), index_load(:,:) + logical :: result + + nrows = 5 + ncols = 5 + + call random_seed() + ! General matrix + call generate_random_sp_coo_matrix(index_save, dummy, nrows, ncols, MS_general, .true.) + call save_mm("test_mmio_sparse.mtx", index_save, [0], format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, dummy) + result = compare_coo_pattern(index_save, index_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=unspecified, type=pattern") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + + ! Symmetric matrix + call generate_random_sp_coo_matrix(index_save, dummy, nrows, ncols, MS_symmetric, .true.) + call save_mm("test_mmio_sparse.mtx", index_save, [0], symmetry = "symmetric", format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, dummy) + result = compare_coo_pattern(index_save, index_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=symmetric, type=pattern") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + + ! Skew-symmetric matrix + call generate_random_sp_coo_matrix(index_save, dummy, nrows, ncols, MS_skew_symmetric, .true.) + call save_mm("test_mmio_sparse.mtx", index_save, [0], symmetry = "skew-symmetric", format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, dummy) + result = compare_coo_pattern(index_save, index_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=skew-symmetric, type=pattern") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + + ! Hermitian matrix + call generate_random_sp_coo_matrix(index_save, dummy, nrows, ncols, MS_hermitian, .true.) + call save_mm("test_mmio_sparse.mtx", index_save, [0], symmetry = "hermitian", format = "G0") + call load_mm("test_mmio_sparse.mtx", index_load, dummy) + result = compare_coo_pattern(index_save, index_load) + call check(error, result, .true.,& + "MM coordinate test failed: symmetry_arg=hermitian, type=pattern") + if(allocated(error)) return + if(allocated(index_save)) deallocate(index_save) + end block + + ! Validate given matrices inside https://sparse.tamu.edu/ + block + real(dp), allocatable :: data_f_load_1(:), data_s_load_1(:), data_f_load_2(:,:), data_s_load_2(:,:) + integer, allocatable :: index_f_load(:,:), index_s_load(:) + integer :: nrows, nnz, stat, i + character(len=:), allocatable :: path, str, filename + + str = "${_FILE_}$" + i = len_trim(str) + do while(i > 0) + if(str(i:i) == '/' .or. str(i:i) == '\') exit + i = i - 1 + end do + path = str(:i) // "../data/matrix_market/" + + ! coordinate real symmetric + filename = "bcsstk01" + call load_mm(path // filename // ".mtx", index_f_load, data_f_load_1, iostat=stat) + nnz = size(index_f_load, dim=2) + call load_npy(path // filename // ".mtx_indices.npy", index_s_load, iostat=stat) + index_s_load = index_s_load + 1 + call load_npy(path // filename //".mtx_data.npy", data_s_load_1, iostat=stat) + call check(error, all(index_f_load(1,:)==index_s_load(1:nnz)), .true.,& + "MM coordinate test failed: bcsstk01.mtx rows not matched") + call check(error, all(index_f_load(2,:)==index_s_load(nnz+1:2*nnz)), .true.,& + "MM coordinate test failed: bcsstk01.mtx cols not matched") + call check(error, all_close(data_f_load_1, data_s_load_1), .true.,& + "MM coordinate test failed: bcsstk01.mtx data not matched") + + ! coordinate pattern symmetric + filename = "ash85" + call load_mm(path // filename // ".mtx", index_f_load, data_f_load_1, iostat=stat) + nnz = size(index_f_load, dim=2) + call load_npy(path // filename // ".mtx_indices.npy", index_s_load, iostat=stat) + index_s_load = index_s_load + 1 + call check(error, all(index_f_load(1,:)==index_s_load(1:nnz)), .true.,& + "MM coordinate test failed: bcsstk01.mtx rows not matched") + call check(error, all(index_f_load(2,:)==index_s_load(nnz+1:2*nnz)), .true.,& + "MM coordinate test failed: bcsstk01.mtx cols not matched") + + ! array real general + filename = "lhr01_b" + call load_mm(path // filename // ".mtx", data_f_load_2, iostat=stat) + call load_npy(path // filename // ".mtx_data.npy", data_s_load_2, iostat=stat) + call check(error, all_close(data_f_load_2, transpose(data_s_load_2)), .true.,& + "MM coordinate test failed: bcsstk01.mtx data not matched") + end block + stat = 0 + open(newunit=u, file="test_mmio_sparse.mtx", status="old", iostat=stat) + if (stat == 0) then + close(u, status="delete") + end if + end subroutine + +end module + + +program tester + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_io_mm, only : collect_suite + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("io_mm", collect_suite) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program