diff --git a/src/conversion/fsparse_conversions.f90 b/src/conversion/fsparse_conversions.f90 index 6c1fd3e..8e36db1 100644 --- a/src/conversion/fsparse_conversions.f90 +++ b/src/conversion/fsparse_conversions.f90 @@ -8,6 +8,11 @@ module fsparse_conversions use fsparse_constants use fsparse_matrix_gallery implicit none + private + public :: dense2coo, dense2diagonal + public :: coo2dense, coo2diagonal + public :: coo2csr, csr2coo, csr2diagonal, csr2sellc + public :: csr_block_expansion interface dense2coo module procedure dense2coo_sp @@ -67,6 +72,13 @@ module fsparse_conversions module procedure csr2sellc_cdp end interface + interface csr_block_expansion + module procedure csr_block_expansion_sp + module procedure csr_block_expansion_dp + module procedure csr_block_expansion_csp + module procedure csr_block_expansion_cdp + end interface + contains subroutine dense2coo_sp(dense,COO) integer, parameter :: wp = sp @@ -942,5 +954,374 @@ subroutine csr2diagonal_cdp(CSR,diagonal) end select end subroutine + + subroutine csr_block_expansion_sp(CSR,num_dof) + !! This function enables expanding an ordered CSR matrix to include contigous blocks + !! Warning: this operation does not preserve the data buffer + type(CSR_sp), intent(inout) :: CSR + integer, intent(in) :: num_dof + real(sp), parameter :: zero = zero_sp + integer, allocatable :: rowptr_expn(:), col_expn(:) + integer :: i, j, p, q, adr1, adr2, block_nnz + + select case(CSR%sym) + case(k_NOSYMMETRY) + CSR%NNZ = CSR%NNZ * num_dof ** 2 + case(k_SYMTRISUP,k_SYMTRIINF) + block_nnz = num_dof + num_dof * (num_dof-1) / 2 + CSR%NNZ = CSR%nrows * block_nnz + (CSR%NNZ-CSR%nrows) * num_dof ** 2 + end select + + allocate( col_expn(CSR%NNZ) ) + allocate( rowptr_expn(CSR%nrows*num_dof+1) ) + + rowptr_expn(1) = 1 + select case(CSR%sym) + case(k_NOSYMMETRY) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRISUP) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - p + 1 + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + j = CSR%rowptr(i) + adr2 = adr1 - p + 1 + do q = p, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + do j = CSR%rowptr(i)+1, CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) - p + 1 + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRIINF) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - num_dof + p + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-2 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + j = CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, p + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end select + + CSR%nrows = CSR%nrows * num_dof + CSR%ncols = CSR%ncols * num_dof + + call move_alloc( rowptr_expn , CSR%rowptr ) + call move_alloc( col_expn , CSR%col ) + + if( allocated(CSR%data) ) deallocate( CSR%data ) + allocate( CSR%data(CSR%NNZ) , source = zero ) + end subroutine + + subroutine csr_block_expansion_dp(CSR,num_dof) + !! This function enables expanding an ordered CSR matrix to include contigous blocks + !! Warning: this operation does not preserve the data buffer + type(CSR_dp), intent(inout) :: CSR + integer, intent(in) :: num_dof + real(dp), parameter :: zero = zero_dp + integer, allocatable :: rowptr_expn(:), col_expn(:) + integer :: i, j, p, q, adr1, adr2, block_nnz + + select case(CSR%sym) + case(k_NOSYMMETRY) + CSR%NNZ = CSR%NNZ * num_dof ** 2 + case(k_SYMTRISUP,k_SYMTRIINF) + block_nnz = num_dof + num_dof * (num_dof-1) / 2 + CSR%NNZ = CSR%nrows * block_nnz + (CSR%NNZ-CSR%nrows) * num_dof ** 2 + end select + + allocate( col_expn(CSR%NNZ) ) + allocate( rowptr_expn(CSR%nrows*num_dof+1) ) + + rowptr_expn(1) = 1 + select case(CSR%sym) + case(k_NOSYMMETRY) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRISUP) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - p + 1 + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + j = CSR%rowptr(i) + adr2 = adr1 - p + 1 + do q = p, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + do j = CSR%rowptr(i)+1, CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) - p + 1 + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRIINF) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - num_dof + p + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-2 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + j = CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, p + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end select + + CSR%nrows = CSR%nrows * num_dof + CSR%ncols = CSR%ncols * num_dof + + call move_alloc( rowptr_expn , CSR%rowptr ) + call move_alloc( col_expn , CSR%col ) + + if( allocated(CSR%data) ) deallocate( CSR%data ) + allocate( CSR%data(CSR%NNZ) , source = zero ) + end subroutine + + subroutine csr_block_expansion_csp(CSR,num_dof) + !! This function enables expanding an ordered CSR matrix to include contigous blocks + !! Warning: this operation does not preserve the data buffer + type(CSR_csp), intent(inout) :: CSR + integer, intent(in) :: num_dof + complex(sp), parameter :: zero = zero_csp + integer, allocatable :: rowptr_expn(:), col_expn(:) + integer :: i, j, p, q, adr1, adr2, block_nnz + + select case(CSR%sym) + case(k_NOSYMMETRY) + CSR%NNZ = CSR%NNZ * num_dof ** 2 + case(k_SYMTRISUP,k_SYMTRIINF) + block_nnz = num_dof + num_dof * (num_dof-1) / 2 + CSR%NNZ = CSR%nrows * block_nnz + (CSR%NNZ-CSR%nrows) * num_dof ** 2 + end select + + allocate( col_expn(CSR%NNZ) ) + allocate( rowptr_expn(CSR%nrows*num_dof+1) ) + + rowptr_expn(1) = 1 + select case(CSR%sym) + case(k_NOSYMMETRY) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRISUP) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - p + 1 + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + j = CSR%rowptr(i) + adr2 = adr1 - p + 1 + do q = p, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + do j = CSR%rowptr(i)+1, CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) - p + 1 + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRIINF) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - num_dof + p + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-2 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + j = CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, p + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end select + + CSR%nrows = CSR%nrows * num_dof + CSR%ncols = CSR%ncols * num_dof + + call move_alloc( rowptr_expn , CSR%rowptr ) + call move_alloc( col_expn , CSR%col ) + + if( allocated(CSR%data) ) deallocate( CSR%data ) + allocate( CSR%data(CSR%NNZ) , source = zero ) + end subroutine + + subroutine csr_block_expansion_cdp(CSR,num_dof) + !! This function enables expanding an ordered CSR matrix to include contigous blocks + !! Warning: this operation does not preserve the data buffer + type(CSR_cdp), intent(inout) :: CSR + integer, intent(in) :: num_dof + complex(dp), parameter :: zero = zero_cdp + integer, allocatable :: rowptr_expn(:), col_expn(:) + integer :: i, j, p, q, adr1, adr2, block_nnz + + select case(CSR%sym) + case(k_NOSYMMETRY) + CSR%NNZ = CSR%NNZ * num_dof ** 2 + case(k_SYMTRISUP,k_SYMTRIINF) + block_nnz = num_dof + num_dof * (num_dof-1) / 2 + CSR%NNZ = CSR%nrows * block_nnz + (CSR%NNZ-CSR%nrows) * num_dof ** 2 + end select + + allocate( col_expn(CSR%NNZ) ) + allocate( rowptr_expn(CSR%nrows*num_dof+1) ) + + rowptr_expn(1) = 1 + select case(CSR%sym) + case(k_NOSYMMETRY) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRISUP) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - p + 1 + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + j = CSR%rowptr(i) + adr2 = adr1 - p + 1 + do q = p, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + do j = CSR%rowptr(i)+1, CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) - p + 1 + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRIINF) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - num_dof + p + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-2 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + j = CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, p + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end select + + CSR%nrows = CSR%nrows * num_dof + CSR%ncols = CSR%ncols * num_dof + + call move_alloc( rowptr_expn , CSR%rowptr ) + call move_alloc( col_expn , CSR%col ) + + if( allocated(CSR%data) ) deallocate( CSR%data ) + allocate( CSR%data(CSR%NNZ) , source = zero ) + end subroutine + end module fsparse_conversions \ No newline at end of file diff --git a/src/conversion/fsparse_conversions.fypp b/src/conversion/fsparse_conversions.fypp index e7ada8d..c533f8f 100644 --- a/src/conversion/fsparse_conversions.fypp +++ b/src/conversion/fsparse_conversions.fypp @@ -10,6 +10,11 @@ module fsparse_conversions use fsparse_constants use fsparse_matrix_gallery implicit none + private + public :: dense2coo, dense2diagonal + public :: coo2dense, coo2diagonal + public :: coo2csr, csr2coo, csr2diagonal, csr2sellc + public :: csr_block_expansion interface dense2coo #:for k1, t1, s1 in (KINDS_TYPES) @@ -61,6 +66,12 @@ module fsparse_conversions #:endfor end interface + interface csr_block_expansion + #:for k1, t1, s1 in (KINDS_TYPES) + module procedure csr_block_expansion_${s1}$ + #:endfor + end interface + contains #:for k1, t1, s1 in (KINDS_TYPES) subroutine dense2coo_${s1}$(dense,COO) @@ -337,5 +348,100 @@ contains end subroutine #:endfor + + #:for k1, t1, s1 in (KINDS_TYPES) + subroutine csr_block_expansion_${s1}$(CSR,num_dof) + !! This function enables expanding an ordered CSR matrix to include contigous blocks + !! Warning: this operation does not preserve the data buffer + type(CSR_${s1}$), intent(inout) :: CSR + integer, intent(in) :: num_dof + ${t1}$, parameter :: zero = zero_${s1}$ + integer, allocatable :: rowptr_expn(:), col_expn(:) + integer :: i, j, p, q, adr1, adr2, block_nnz + + select case(CSR%sym) + case(k_NOSYMMETRY) + CSR%NNZ = CSR%NNZ * num_dof ** 2 + case(k_SYMTRISUP,k_SYMTRIINF) + block_nnz = num_dof + num_dof * (num_dof-1) / 2 + CSR%NNZ = CSR%nrows * block_nnz + (CSR%NNZ-CSR%nrows) * num_dof ** 2 + end select + + allocate( col_expn(CSR%NNZ) ) + allocate( rowptr_expn(CSR%nrows*num_dof+1) ) + + rowptr_expn(1) = 1 + select case(CSR%sym) + case(k_NOSYMMETRY) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRISUP) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - p + 1 + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + j = CSR%rowptr(i) + adr2 = adr1 - p + 1 + do q = p, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + do j = CSR%rowptr(i)+1, CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) - p + 1 + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end do + case(k_SYMTRIINF) + do i = 1, CSR%nrows + do p = 1, num_dof + rowptr_expn(num_dof*(i-1)+p+1) = rowptr_expn(num_dof*(i-1)+p) & + & + num_dof*(CSR%rowptr(i+1)-CSR%rowptr(i)) - num_dof + p + end do + do p = 1, num_dof + adr1 = rowptr_expn(num_dof*(i-1)+p) + do j = CSR%rowptr(i), CSR%rowptr(i+1)-2 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, num_dof + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + j = CSR%rowptr(i+1)-1 + adr2 = adr1 + num_dof*(j-CSR%rowptr(i)) + do q = 1, p + col_expn(adr2+q-1) = num_dof*(CSR%col(j)-1)+q + end do + end do + end do + end select + + CSR%nrows = CSR%nrows * num_dof + CSR%ncols = CSR%ncols * num_dof + + call move_alloc( rowptr_expn , CSR%rowptr ) + call move_alloc( col_expn , CSR%col ) + + if( allocated(CSR%data) ) deallocate( CSR%data ) + allocate( CSR%data(CSR%NNZ) , source = zero ) + end subroutine + + #:endfor end module fsparse_conversions \ No newline at end of file diff --git a/src/krylov/fsparse_krylov.f90 b/src/krylov/fsparse_krylov.f90 index 5829a78..bf37acc 100644 --- a/src/krylov/fsparse_krylov.f90 +++ b/src/krylov/fsparse_krylov.f90 @@ -53,24 +53,6 @@ real(dp) function abs_dot_dp(x,y) result(r) public :: cgsolve interface cgsolve - module subroutine cgsolve_COO_sp(A,b,x,di,rtol,maxiter,restart) - type(COO_sp), intent(in), target :: A - real(sp), intent(in) :: b(:) - real(sp), intent(out) :: x(:) - logical, intent(in), optional, target :: di(:) !! Dirichlet conditions for blocked d.o.f - real(sp), intent(in), optional :: rtol - integer, intent(in), optional :: maxiter - logical, intent(in), optional :: restart !! Set X=0 or use it as starting solution - end subroutine - module subroutine cgsolve_COO_dp(A,b,x,di,rtol,maxiter,restart) - type(COO_dp), intent(in), target :: A - real(dp), intent(in) :: b(:) - real(dp), intent(out) :: x(:) - logical, intent(in), optional, target :: di(:) !! Dirichlet conditions for blocked d.o.f - real(dp), intent(in), optional :: rtol - integer, intent(in), optional :: maxiter - logical, intent(in), optional :: restart !! Set X=0 or use it as starting solution - end subroutine module subroutine cgsolve_CSR_sp(A,b,x,di,rtol,maxiter,restart) type(CSR_sp), intent(in), target :: A real(sp), intent(in) :: b(:) @@ -189,6 +171,12 @@ subroutine loc_matvec_sp(op,x,y) y = zero_sp call matvec(op%ptr_CSR_sp,x,y) end subroutine + real(sp) function dot_sp(x,y) result(r) + real(sp), intent(in) :: x(:) + real(sp), intent(in) :: y(:) + r = dot_product(x,y) + end function + subroutine loc_precond_sp(op,x,y) class(linop_sp) :: op real(sp), intent(in) :: x(:) @@ -196,11 +184,6 @@ subroutine loc_precond_sp(op,x,y) y = zero_sp call matvec(op%ptr_pcd_sp,x,y) end subroutine - real(sp) function dot_sp(x,y) result(r) - real(sp), intent(in) :: x(:) - real(sp), intent(in) :: y(:) - r = dot_product(x,y) - end function subroutine factorization_jacobi_sp(op) class(linop_sp) :: op @@ -226,6 +209,12 @@ subroutine loc_matvec_dp(op,x,y) y = zero_dp call matvec(op%ptr_CSR_dp,x,y) end subroutine + real(dp) function dot_dp(x,y) result(r) + real(dp), intent(in) :: x(:) + real(dp), intent(in) :: y(:) + r = dot_product(x,y) + end function + subroutine loc_precond_dp(op,x,y) class(linop_dp) :: op real(dp), intent(in) :: x(:) @@ -233,11 +222,6 @@ subroutine loc_precond_dp(op,x,y) y = zero_dp call matvec(op%ptr_pcd_dp,x,y) end subroutine - real(dp) function dot_dp(x,y) result(r) - real(dp), intent(in) :: x(:) - real(dp), intent(in) :: y(:) - r = dot_product(x,y) - end function subroutine factorization_jacobi_dp(op) class(linop_dp) :: op @@ -383,6 +367,12 @@ subroutine loc_matvec_sp(op,x,y) y = zero_sp call matvec(op%ptr_dense_sp,x,y) end subroutine + real(sp) function dot_sp(x,y) result(r) + real(sp), intent(in) :: x(:) + real(sp), intent(in) :: y(:) + r = dot_product(x,y) + end function + subroutine loc_precond_sp(op,x,y) class(linop_sp) :: op real(sp), intent(in) :: x(:) @@ -390,11 +380,6 @@ subroutine loc_precond_sp(op,x,y) y = zero_sp call matvec(op%ptr_pcd_sp,x,y) end subroutine - real(sp) function dot_sp(x,y) result(r) - real(sp), intent(in) :: x(:) - real(sp), intent(in) :: y(:) - r = dot_product(x,y) - end function subroutine factorization_jacobi_sp(op) class(linop_sp) :: op @@ -420,6 +405,12 @@ subroutine loc_matvec_dp(op,x,y) y = zero_dp call matvec(op%ptr_dense_dp,x,y) end subroutine + real(dp) function dot_dp(x,y) result(r) + real(dp), intent(in) :: x(:) + real(dp), intent(in) :: y(:) + r = dot_product(x,y) + end function + subroutine loc_precond_dp(op,x,y) class(linop_dp) :: op real(dp), intent(in) :: x(:) @@ -427,11 +418,6 @@ subroutine loc_precond_dp(op,x,y) y = zero_dp call matvec(op%ptr_pcd_dp,x,y) end subroutine - real(dp) function dot_dp(x,y) result(r) - real(dp), intent(in) :: x(:) - real(dp), intent(in) :: y(:) - r = dot_product(x,y) - end function subroutine factorization_jacobi_dp(op) class(linop_dp) :: op @@ -482,6 +468,8 @@ module subroutine cgsolve_dense_sp(A,b,x,di,rtol,maxiter,restart) !----------------------------------------------------- ! set matrix pointer op%ptr_dense_sp%data => A(:,:) + op%ptr_dense_sp%nrows = size(A,dim=1) + op%ptr_dense_sp%ncols = size(A,dim=2) !----------------------------------------------------- ! set preconditioner @@ -528,6 +516,8 @@ module subroutine cgsolve_dense_dp(A,b,x,di,rtol,maxiter,restart) !----------------------------------------------------- ! set matrix pointer op%ptr_dense_dp%data => A(:,:) + op%ptr_dense_dp%nrows = size(A,dim=1) + op%ptr_dense_dp%ncols = size(A,dim=2) !----------------------------------------------------- ! set preconditioner diff --git a/src/krylov/fsparse_krylov.fypp b/src/krylov/fsparse_krylov.fypp index dc1154d..125ce15 100644 --- a/src/krylov/fsparse_krylov.fypp +++ b/src/krylov/fsparse_krylov.fypp @@ -44,7 +44,7 @@ module fsparse_krylov interface cgsolve #:for ST1 in (SPARSE_TYPES) - #:if ST1=="COO" or ST1=="CSR" or ST1=="dense" + #:if ST1=="CSR" or ST1=="dense" #:for k1, t1, s1 in (KINDS_TYPES) module subroutine cgsolve_${ST1}$_${s1}$(A,b,x,di,rtol,maxiter,restart) #:if ST1 == "dense" @@ -132,7 +132,7 @@ submodule (fsparse_krylov) fsparse_krylov_${ST1}$ ${t1}$, intent(in) :: y(:) r = dot_product(x,y) end function - + subroutine loc_precond_${s1}$(op,x,y) class(linop_${s1}$) :: op ${t1}$, intent(in) :: x(:) @@ -201,6 +201,8 @@ submodule (fsparse_krylov) fsparse_krylov_${ST1}$ ! set matrix pointer #:if ST1 == "dense" op%ptr_${ST1}$_${s1}$%data => A(:,:) + op%ptr_${ST1}$_${s1}$%nrows = size(A,dim=1) + op%ptr_${ST1}$_${s1}$%ncols = size(A,dim=2) #:else op%ptr_${ST1}$_${s1}$ => A #:endif diff --git a/test/test_matrices.f90 b/test/test_matrices.f90 index c26f71f..fe2e9c7 100644 --- a/test/test_matrices.f90 +++ b/test/test_matrices.f90 @@ -7,6 +7,7 @@ module test_fsparse use testdrive, only: new_unittest, unittest_type, error_type, check use fsparse + use fsparse_constants implicit none contains @@ -406,16 +407,38 @@ subroutine test_symmetries(error) subroutine test_cells2sparse(error) !> Error handling type(error_type), allocatable, intent(out) :: error - type(COO_dp) :: COO - integer :: cells(4,2) + type(COO_dp) :: COO, COO_n + type(CSR_dp) :: CSR, CSR_n1 + real(8), allocatable :: dense(:,:), dense_n(:,:) + integer :: i, k, dof + integer, allocatable :: cells(:,:), cells_n(:,:) + allocate(cells(4,2)) cells(1:4,1) = [2,5,3,4] cells(1:4,2) = [4,1,3,2] - call coo_from_cells(COO,cells,num_points=5,num_cells=2,selfloop=.true.,symtype=k_SYMTRISUP) - + call coo_from_cells(COO,cells,num_points=5,num_cells=2,selfloop=.true.,symtype=k_SYMTRIINF) + call coo2csr(COO,CSR) + call check(error, size(COO%data) == 14 .and. size( COO%index , dim=2 ) == 14 ) if (allocated(error)) return + dof = 4 + allocate(cells_n(4*dof,2)) + do i =1 , 4 + do k = 1, dof + cells_n(dof*(i-1)+k,1) = dof*(cells(i,1)-1)+k + cells_n(dof*(i-1)+k,2) = dof*(cells(i,2)-1)+k + end do + end do + + call coo_from_cells(COO_n,cells_n,num_points=5*dof,num_cells=2,selfloop=.true.,symtype=k_SYMTRIINF) + call coo2csr(COO_n,CSR_n1) + + call csr_block_expansion(CSR,dof) + + call check(error, all(CSR%rowptr==CSR_n1%rowptr) .and. all(CSR%col==CSR_n1%col) ) + if (allocated(error)) return + end subroutine end module test_fsparse \ No newline at end of file diff --git a/test/test_matrices.fypp b/test/test_matrices.fypp index 980cf32..c4e79b3 100644 --- a/test/test_matrices.fypp +++ b/test/test_matrices.fypp @@ -9,6 +9,7 @@ module test_fsparse use testdrive, only: new_unittest, unittest_type, error_type, check use fsparse + use fsparse_constants implicit none contains @@ -266,16 +267,38 @@ module test_fsparse subroutine test_cells2sparse(error) !> Error handling type(error_type), allocatable, intent(out) :: error - type(COO_dp) :: COO - integer :: cells(4,2) + type(COO_dp) :: COO, COO_n + type(CSR_dp) :: CSR, CSR_n1 + real(8), allocatable :: dense(:,:), dense_n(:,:) + integer :: i, k, dof + integer, allocatable :: cells(:,:), cells_n(:,:) + allocate(cells(4,2)) cells(1:4,1) = [2,5,3,4] cells(1:4,2) = [4,1,3,2] - call coo_from_cells(COO,cells,num_points=5,num_cells=2,selfloop=.true.,symtype=k_SYMTRISUP) - + call coo_from_cells(COO,cells,num_points=5,num_cells=2,selfloop=.true.,symtype=k_SYMTRIINF) + call coo2csr(COO,CSR) + call check(error, size(COO%data) == 14 .and. size( COO%index , dim=2 ) == 14 ) if (allocated(error)) return + dof = 4 + allocate(cells_n(4*dof,2)) + do i =1 , 4 + do k = 1, dof + cells_n(dof*(i-1)+k,1) = dof*(cells(i,1)-1)+k + cells_n(dof*(i-1)+k,2) = dof*(cells(i,2)-1)+k + end do + end do + + call coo_from_cells(COO_n,cells_n,num_points=5*dof,num_cells=2,selfloop=.true.,symtype=k_SYMTRIINF) + call coo2csr(COO_n,CSR_n1) + + call csr_block_expansion(CSR,dof) + + call check(error, all(CSR%rowptr==CSR_n1%rowptr) .and. all(CSR%col==CSR_n1%col) ) + if (allocated(error)) return + end subroutine end module test_fsparse \ No newline at end of file diff --git a/test/test_solvers.f90 b/test/test_solvers.f90 index 2e94077..5b9590d 100644 --- a/test/test_solvers.f90 +++ b/test/test_solvers.f90 @@ -28,7 +28,6 @@ subroutine test_pccg(error) real :: A(2,2) = reshape([4.0,1.0,1.0,3.0],shape=[2,2]) real :: b(2) = [1.0,2.0] real :: x(2) - real :: tol = 0.0001 call cgsolve(A,b,x)