From 5ee805b29bd67311f9547671c52ea85236c41ac5 Mon Sep 17 00:00:00 2001 From: Jason Christopherson Date: Tue, 30 Jul 2024 06:20:25 -0500 Subject: [PATCH] Add sparse matrix construction routine --- src/linalg.f90 | 10 ++++++ src/linalg_sparse.f90 | 73 +++++++++++++++++++++++++++++++++++++++++++ src/sparskit.f90 | 66 ++++++++++++++++++++++++++++++++++++++ tests/test_sparse.f90 | 22 ++++++++++++- 4 files changed, 170 insertions(+), 1 deletion(-) diff --git a/src/linalg.f90 b/src/linalg.f90 index 71834e5..cb89f1f 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -204,6 +204,7 @@ module linalg public :: msr_to_csr public :: dense_to_msr public :: msr_to_dense + public :: create_csr_matrix public :: matmul public :: operator(+) public :: operator(-) @@ -5738,6 +5739,15 @@ module subroutine csr_pgmres_solver(a, lu, ju, b, x, im, tol, maxits, & real(real64), intent(in), optional :: tol class(errors), intent(inout), optional, target :: err end subroutine + + module function create_csr_matrix(m, n, rows, cols, vals, err) & + result(rst) + integer(int32), intent(in) :: m, n + integer(int32), intent(in), dimension(:) :: rows, cols + real(real64), intent(in), dimension(:) :: vals + class(errors), intent(inout), optional, target :: err + type(csr_matrix) :: rst + end function end interface ! ------------------------------------------------------------------------------ diff --git a/src/linalg_sparse.f90 b/src/linalg_sparse.f90 index e9a6cb9..977a607 100644 --- a/src/linalg_sparse.f90 +++ b/src/linalg_sparse.f90 @@ -1138,6 +1138,79 @@ module subroutine msr_assign_to_csr(csr, msr) csr = msr_to_csr(msr) end subroutine +! ------------------------------------------------------------------------------ +module function create_csr_matrix(m, n, rows, cols, vals, err) result(rst) + ! Arguments + integer(int32), intent(in) :: m, n + integer(int32), intent(in), dimension(:) :: rows, cols + real(real64), intent(in), dimension(:) :: vals + class(errors), intent(inout), optional, target :: err + type(csr_matrix) :: rst + + ! Local Variables + integer(int32) :: i, flag, nnz + integer(int32), allocatable, dimension(:) :: ir + class(errors), pointer :: errmgr + type(errors), target :: deferr + + ! Initialization + if (present(err)) then + errmgr => err + else + errmgr => deferr + end if + nnz = size(rows) + + ! Input Checking + if (m < 0) then + call errmgr%report_error("create_csr_matrix", & + "The number of rows must be a positive value.", & + LA_INVALID_INPUT_ERROR) + return + end if + if (n < 0) then + call errmgr%report_error("create_csr_matrix", & + "The number of columns must be a positive value.", & + LA_INVALID_INPUT_ERROR) + return + end if + if (size(cols) /= nnz .or. size(vals) /= nnz) then + call errmgr%report_error("create_csr_matrix", & + "The size of the input arrays must be the same.", & + LA_ARRAY_SIZE_ERROR) + return + end if + do i = 1, nnz + if (rows(i) < 1 .or. rows(i) > m) then + call errmgr%report_error("create_csr_matrix", & + "All row indices must be within the bounds of the matrix.", & + LA_INVALID_INPUT_ERROR) + return + end if + if (cols(i) < 1 .or. cols(i) > n) then + call errmgr%report_error("create_csr_matrix", & + "All column indices must be within the bounds of the matrix.", & + LA_INVALID_INPUT_ERROR) + return + end if + end do + allocate(ir(nnz), source = rows, stat = flag) + if (flag /= 0) then + call errmgr%report_error("create_csr_matrix", & + "Memory allocation error.", LA_OUT_OF_MEMORY_ERROR) + return + end if + + ! Create an empty matrix + rst = create_empty_csr_matrix(m, n, nnz, errmgr) + if (errmgr%has_error_occurred()) return + + ! Populate the empty matrix + call coocsr(m, nnz, vals, ir, cols, rst%values, rst%column_indices, & + rst%row_indices) + call csort(m, rst%values, rst%column_indices, rst%row_indices, .true.) +end function + ! ****************************************************************************** ! LU PRECONDITIONER ROUTINES ! ------------------------------------------------------------------------------ diff --git a/src/sparskit.f90 b/src/sparskit.f90 index 00f78e8..0ed76bb 100644 --- a/src/sparskit.f90 +++ b/src/sparskit.f90 @@ -195,6 +195,28 @@ subroutine msrcsr(n, a, ja, ao, jao, iao, wk, iwk) real(real64), intent(in) :: a(*) real(real64), intent(out) :: ao(*), wk(n) end subroutine + + !> @brief Converte a matrix stored in coordinate format to CSR format. + !! + !! @param[in] nrow The number of rows in the matrix. + !! @param[in] nnz The number of non-zero elements in the matrix. + !! @param[in] a An NNZ-element array containing the non-zero elements + !! of the matrix. + !! @param[in,out] ir An NNZ-element array containing the row indices of + !! each non-zero element. + !! @param[in] jc An NNZ-element array containing the column indices of + !! each non-zero element. + !! @param[out] ao The non-zero elements of matrix A. + !! @param[out] jao The column indices of matrix A. + !! @param[out] iao The index in A where the requested row starts. + subroutine coocsr(nrow, nnz, a, ir, jc, ao, jao, iao) + use iso_fortran_env, only : int32, real64 + integer(int32), intent(in) :: nrow, nnz, jc(*) + integer(int32), intent(inout) :: ir(*) + real(real64), intent(in) :: a(*) + integer(int32), intent(out) :: jao(*), iao(*) + real(real64) :: ao(*) + end subroutine end interface ! UNARY.F @@ -238,6 +260,50 @@ subroutine getdia(nrow, ncol, job, a, ja, ia, len, diag, idiag, ioff) real(real64), intent(in) :: a(*) real(real64), intent(out) :: diag(*) end subroutine + + !> @brief Sorces the elements of a CSR matrix in increasing order of + !! their column indices within each row. + !! + !! @param[in] n The number of rows in the matrix. + !! @param[in,out] a The non-zero values. + !! @param[in,out] ja An array of column indices of the elements in A. + !! @param[in] ia An array of pointers to the rows. + !! @param[in] values Idicates whether A must also be permuted. If + !! false, A can be a dummy array. + subroutine csort(n, a, ja, ia, values) + use iso_fortran_env, only : int32, real64 + integer(int32), intent(in) :: n + real(real64), intent(inout) :: a(*) + integer(int32), intent(inout) :: ja(*) + integer(int32), intent(in) :: ia(*) + logical, intent(in) :: values + end subroutine + + !> @breif Cleans up a CSR matrix. + !! + !! @param[in] job The job to perform. + !! - 0: Nothing is done + !! - 1: Eliminate duplicate entries and zero entries. + !! - 2: Eliminate duplicate entries and perform partial ordering. + !! - 3: Eliminate duplicate entries and sort the entries in increasing + !! order of column indices. + !! @param[in] value2 0 if the matrix is pattern only (A is not touched), + !! or 1 if the matrix has values. + !! @param[in] nrow The number of rows in the matrix. + !! @param[in,out] a The non-zero values. + !! @param[in,out] ja An array of column indices of the elements in A. + !! @param[in,out] ia An array of pointers to the rows. + !! @param[out] indu An NROW array containing pointers to the beginning + !! of the upper triangular portion if job > 1. + !! @param[out] iwk An NROW+1 element workspace array. + + subroutine clncsr(job, value2, nrow, a, ja, ia, indu, iwk) + use iso_fortran_env, only : int32, real64 + integer(int32), intent(in) :: job, value2, nrow + real(real64), intent(inout) :: a(*) + integer(int32), intent(inout) :: ja(*), ia(*) + integer(int32), intent(inout) :: indu(*), iwk(*) + end subroutine end interface ! ILUT.F diff --git a/tests/test_sparse.f90 b/tests/test_sparse.f90 index 256b132..0acbedb 100644 --- a/tests/test_sparse.f90 +++ b/tests/test_sparse.f90 @@ -15,13 +15,18 @@ function test_csr_1() result(rst) integer(int32), parameter :: n = 6 integer(int32), parameter :: nnz = 8 real(real64) :: dense(m, n), v(nnz), check(m, n) - integer(int32) :: ja(nnz), ia(m + 1) + integer(int32) :: ja(nnz), ia(m + 1), rows(nnz), cols(nnz) type(csr_matrix) :: sparse ! Initialization rst = .true. ! Construct a small sparse matrix, but in dense form + ! + ! | 1 2 0 0 0 0| + ! | 0 3 0 4 0 0| + ! | 0 0 5 6 7 0| + ! | 0 0 0 0 0 8| dense = reshape([ & 1.0d0, 0.0d0, 0.0d0, 0.0d0, & 2.0d0, 3.0d0, 0.0d0, 0.0d0, & @@ -64,6 +69,21 @@ function test_csr_1() result(rst) rst = .false. print "(A)", "Test Failed: test_csr_1 -6" end if + + ! ---------- + ! Construct from coordinate format + rows = [1, 1, 2, 2, 3, 3, 3, 4] + cols = [1, 2, 2, 4, 3, 4, 5, 6] + sparse = create_csr_matrix(m, n, rows, cols, v) + + ! Convert back to dense + check = sparse + + ! Test + if (.not.assert(dense, check)) then + rst = .false. + print "(A)", "Test Failed: test_csr_1 -7" + end if end function ! ------------------------------------------------------------------------------