Add sparse matrix construction routine

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
 
 ! ------------------------------------------------------------------------------

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
 ! ------------------------------------------------------------------------------

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

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
 
 ! ------------------------------------------------------------------------------