From 716b3c5535c5937535d95e647639554599ab70ff Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 2 Mar 2025 22:26:52 +0100 Subject: [PATCH 01/26] start iterative solvers --- src/CMakeLists.txt | 2 + src/stdlib_linalg_iterative_solvers.fypp | 77 ++++++++++++++ src/stdlib_linalg_iterative_solvers_cg.fypp | 105 ++++++++++++++++++++ src/stdlib_sparse_constants.fypp | 11 +- 4 files changed, 185 insertions(+), 10 deletions(-) create mode 100644 src/stdlib_linalg_iterative_solvers.fypp create mode 100644 src/stdlib_linalg_iterative_solvers_cg.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 933da34de..b922109b3 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -92,6 +92,8 @@ set(cppFiles stdlib_linalg_blas.fypp stdlib_linalg_lapack.fypp + stdlib_linalg_iterative_solvers.fypp + stdlib_linalg_iterative_solvers_cg.fypp ) add_subdirectory(blas) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp new file mode 100644 index 000000000..c3aca6258 --- /dev/null +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -0,0 +1,77 @@ +#: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 MATRIX_TYPES = ["dense", "CSR"] +#:set RANKS = range(1, 2+1) + +module stdlib_linalg_iterative_solvers + use stdlib_kinds + use stdlib_sparse + implicit none + private + + #:for k, t, s in R_KINDS_TYPES + type, public :: linop_${s}$ + procedure(vector_sub_${s}$), nopass, pointer :: matvec => null() + procedure(vector_sub_${s}$), nopass, pointer :: preconditionner => null() + procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => null() + end type + #:endfor + + #:for k, t, s in R_KINDS_TYPES + type, abstract, public :: solver_workspace_${s}$ + end type + type, public, extends(solver_workspace_${s}$) :: cg_workspace_${s}$ + ${t}$, allocatable :: r(:), p(:), Ap(:) + end type + + #:endfor + + abstract interface + #:for k, t, s in R_KINDS_TYPES + subroutine vector_sub_${s}$(x,y) + import :: ${k}$ + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + end subroutine + pure ${t}$ function reduction_sub_${s}$(x,y) result(r) + import :: ${k}$ + ${t}$, intent(in) :: x(:) + ${t}$, intent(in) :: y(:) + end function + #:endfor + end interface + + interface solve_cg_generic + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) + class(linop_${s}$), intent(in) :: A + ${t}$, intent(in) :: b(:), tol + ${t}$, intent(inout) :: x(:) + integer, intent(in) :: maxiter + type(cg_workspace_${s}$), intent(inout) :: workspace + end subroutine + #:endfor + end interface + public :: solve_cg_generic + + interface solve_cg + #:for matrix in MATRIX_TYPES + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_cg_${matrix}$_${s}$(A,b,x,tol,maxiter,workspace) + #:if matrix == "dense" + ${t}$, intent(in) :: A(:,:) + #:else + type(${matrix}$_${s}$_type), intent(in) :: A + #:endif + ${t}$, intent(in) :: b(:), tol + ${t}$, intent(inout) :: x(:) + integer, intent(in), optional :: maxiter + type(cg_workspace_${s}$), optional, intent(inout), target :: workspace + end subroutine + #:endfor + #:endfor + end interface + public :: solve_cg + +end module stdlib_linalg_iterative_solvers diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp new file mode 100644 index 000000000..a00cbe43c --- /dev/null +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -0,0 +1,105 @@ +#: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 MATRIX_TYPES = ["dense", "CSR"] +#:set RANKS = range(1, 2+1) + +submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg + use stdlib_kinds + use stdlib_sparse + use stdlib_constants + use stdlib_intrinsics, only: dot_product => stdlib_dot_product + use stdlib_linalg_iterative_solvers + implicit none + +contains + + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) + class(linop_${s}$), intent(in) :: A + ${t}$, intent(in) :: b(:), tol + ${t}$, intent(inout) :: x(:) + integer, intent(in) :: maxiter + type(cg_workspace_${s}$), intent(inout) :: workspace + !------------------------- + integer :: iter + ${t}$ :: rtr, rtrold, alpha, beta, norm0_sq + !------------------------- + associate( p => workspace%p, & + r => workspace%r, & + Ap => workspace%Ap) + x = zero_${s}$ + rtr = A%inner_product(r, r) + norm0_sq = A%inner_product(b, b) + p = b + beta = zero_${s}$ + iter = 1 + do while( rtr > tol**2 * norm0_sq .and. iter < maxiter) + p = r + beta * p + call A%matvec(p,Ap) + alpha = rtr / A%inner_product(p, Ap) + x = x + alpha * p + r = r - alpha * Ap + rtrold = rtr + rtr = A%inner_product(r, r) + beta = rtr / rtrold + iter = iter + 1 + end do + end associate + end subroutine + #:endfor + + #:for matrix in MATRIX_TYPES + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_cg_${matrix}$_${s}$(A,b,x,tol,maxiter,workspace) + #:if matrix == "dense" + ${t}$, intent(in) :: A(:,:) + #:else + type(${matrix}$_${s}$_type), intent(in) :: A + #:endif + ${t}$, intent(in) :: b(:), tol + ${t}$, intent(inout) :: x(:) + integer, intent(in), optional :: maxiter + type(cg_workspace_${s}$), optional, intent(inout), target :: workspace + !------------------------- + type(linop_${s}$) :: op + type(cg_workspace_${s}$), pointer :: workspace_ + integer :: n, maxiter_ + !------------------------- + n = size(b) + op%matvec => default_matvec + op%inner_product => default_dot + + maxiter_ = n + if(present(maxiter)) maxiter_ = maxiter + if(present(workspace)) then + workspace_ => workspace + else + allocate( workspace_%r(n), & + workspace_%p(n), & + workspace_%Ap(n)) + + end if + call solve_cg_generic(op,b,x,tol,maxiter_,workspace_) + contains + + subroutine default_matvec(x,y) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + #:if matrix == "dense" + y = matmul(A,x) + #:else + call spmv( A , x, y ) + #:endif + end subroutine + pure ${t}$ function default_dot(x,y) result(r) + ${t}$, intent(in) :: x(:) + ${t}$, intent(in) :: y(:) + r = dot_product(x,y) + end function + end subroutine + + #:endfor + #:endfor + +end submodule stdlib_linalg_iterative_cg \ No newline at end of file diff --git a/src/stdlib_sparse_constants.fypp b/src/stdlib_sparse_constants.fypp index 1e8374bd9..7b5dc4497 100644 --- a/src/stdlib_sparse_constants.fypp +++ b/src/stdlib_sparse_constants.fypp @@ -3,7 +3,7 @@ #:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX)) module stdlib_sparse_constants use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp - + use stdlib_constants implicit none public @@ -20,13 +20,4 @@ module stdlib_sparse_constants ! Integer size support for ILP64 builds should be done here integer, parameter :: ilp = int32 - #:for k1, t1, s1 in (R_KINDS_TYPES) - ${t1}$, parameter :: zero_${s1}$ = 0._${k1}$ - ${t1}$, parameter :: one_${s1}$ = 1._${k1}$ - #:endfor - #:for k1, t1, s1 in (C_KINDS_TYPES) - ${t1}$, parameter :: zero_${s1}$ = (0._${k1}$,0._${k1}$) - ${t1}$, parameter :: one_${s1}$ = (1._${k1}$,1._${k1}$) - #:endfor - end module stdlib_sparse_constants From 9ed419f5e164b7ce51fee7af7cd0209833519cd4 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Mon, 3 Mar 2025 22:59:56 +0100 Subject: [PATCH 02/26] simplify workspace --- src/stdlib_linalg_iterative_solvers.fypp | 10 ++++------ src/stdlib_linalg_iterative_solvers_cg.fypp | 17 +++++++---------- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index c3aca6258..6ce80aa4f 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -19,10 +19,8 @@ module stdlib_linalg_iterative_solvers #:endfor #:for k, t, s in R_KINDS_TYPES - type, abstract, public :: solver_workspace_${s}$ - end type - type, public, extends(solver_workspace_${s}$) :: cg_workspace_${s}$ - ${t}$, allocatable :: r(:), p(:), Ap(:) + type, public :: solver_workspace_${s}$ + ${t}$, allocatable :: tmp(:,:) end type #:endfor @@ -49,7 +47,7 @@ module stdlib_linalg_iterative_solvers ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) integer, intent(in) :: maxiter - type(cg_workspace_${s}$), intent(inout) :: workspace + type(solver_workspace_${s}$), intent(inout) :: workspace end subroutine #:endfor end interface @@ -67,7 +65,7 @@ module stdlib_linalg_iterative_solvers ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) integer, intent(in), optional :: maxiter - type(cg_workspace_${s}$), optional, intent(inout), target :: workspace + type(solver_workspace_${s}$), optional, intent(inout), target :: workspace end subroutine #:endfor #:endfor diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index a00cbe43c..4770076f8 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -20,14 +20,14 @@ contains ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) integer, intent(in) :: maxiter - type(cg_workspace_${s}$), intent(inout) :: workspace + type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter ${t}$ :: rtr, rtrold, alpha, beta, norm0_sq !------------------------- - associate( p => workspace%p, & - r => workspace%r, & - Ap => workspace%Ap) + associate( p => workspace%tmp(:,1), & + r => workspace%tmp(:,2), & + Ap => workspace%tmp(:,3)) x = zero_${s}$ rtr = A%inner_product(r, r) norm0_sq = A%inner_product(b, b) @@ -60,10 +60,10 @@ contains ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) integer, intent(in), optional :: maxiter - type(cg_workspace_${s}$), optional, intent(inout), target :: workspace + type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !------------------------- type(linop_${s}$) :: op - type(cg_workspace_${s}$), pointer :: workspace_ + type(solver_workspace_${s}$), pointer :: workspace_ integer :: n, maxiter_ !------------------------- n = size(b) @@ -75,10 +75,7 @@ contains if(present(workspace)) then workspace_ => workspace else - allocate( workspace_%r(n), & - workspace_%p(n), & - workspace_%Ap(n)) - + allocate( workspace_%tmp(n,3) ) end if call solve_cg_generic(op,b,x,tol,maxiter_,workspace_) contains From 9106676fc85b48e87ef31dc2eb2f049ca1ba91e4 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Sat, 8 Mar 2025 15:21:47 +0100 Subject: [PATCH 03/26] add pccg solver and example --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_solve_pccg.f90 | 32 ++++ src/CMakeLists.txt | 1 + src/stdlib_linalg_iterative_solvers.fypp | 44 ++++- src/stdlib_linalg_iterative_solvers_cg.fypp | 8 +- src/stdlib_linalg_iterative_solvers_pccg.fypp | 161 ++++++++++++++++++ 6 files changed, 244 insertions(+), 3 deletions(-) create mode 100644 example/linalg/example_solve_pccg.f90 create mode 100644 src/stdlib_linalg_iterative_solvers_pccg.fypp diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 10f982a02..d2db30042 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -41,6 +41,7 @@ ADD_EXAMPLE(get_norm) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) +ADD_EXAMPLE(solve_pccg) ADD_EXAMPLE(sparse_from_ijv) ADD_EXAMPLE(sparse_data_accessors) ADD_EXAMPLE(sparse_spmv) diff --git a/example/linalg/example_solve_pccg.f90 b/example/linalg/example_solve_pccg.f90 new file mode 100644 index 000000000..2826e9349 --- /dev/null +++ b/example/linalg/example_solve_pccg.f90 @@ -0,0 +1,32 @@ +program example_solve_pccg + use stdlib_kinds, only: dp + use stdlib_sparse + use stdlib_linalg_iterative_solvers, only: solve_pccg + + type(CSR_dp_type) :: laplacian_csr + type(COO_dp_type) :: COO + real(dp) :: laplacian(5,5) + real(dp) :: x(5), load(5) + logical(1) :: dirichlet(5) + + laplacian = reshape( [1, -1, 0, 0, 0,& + -1, 2, -1, 0, 0,& + 0, -1, 2, -1, 0,& + 0, 0, -1, 2, -1,& + 0, 0, 0, -1, 1] , [5,5]) + call dense2coo(laplacian,COO) + call coo2csr(COO,laplacian_csr) + + x = 0._dp + load = dble( [0,0,5,0,0] ) + + dirichlet = .false._1 + dirichlet([1,5]) = .true._1 + + call solve_pccg(laplacian, load, x, tol=1.d-6, di=dirichlet) + print *, x + x = 0._dp + + call solve_pccg(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) + print *, x +end program \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b922109b3..88c82608d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -94,6 +94,7 @@ set(cppFiles stdlib_linalg_lapack.fypp stdlib_linalg_iterative_solvers.fypp stdlib_linalg_iterative_solvers_cg.fypp + stdlib_linalg_iterative_solvers_pccg.fypp ) add_subdirectory(blas) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 6ce80aa4f..ab196dfe4 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -44,8 +44,9 @@ module stdlib_linalg_iterative_solvers #:for k, t, s in R_KINDS_TYPES module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A - ${t}$, intent(in) :: b(:), tol + ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) + ${t}$, intent(in) :: tol integer, intent(in) :: maxiter type(solver_workspace_${s}$), intent(inout) :: workspace end subroutine @@ -62,8 +63,9 @@ module stdlib_linalg_iterative_solvers #:else type(${matrix}$_${s}$_type), intent(in) :: A #:endif - ${t}$, intent(in) :: b(:), tol + ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) + ${t}$, intent(in), optional :: tol integer, intent(in), optional :: maxiter type(solver_workspace_${s}$), optional, intent(inout), target :: workspace end subroutine @@ -71,5 +73,43 @@ module stdlib_linalg_iterative_solvers #:endfor end interface public :: solve_cg + + interface solve_pccg_generic + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_pccg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + class(linop_${s}$), intent(in) :: A + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$, intent(in) :: tol + logical(1), intent(in) :: di(:) + integer, intent(in) :: maxiter + logical, intent(in) :: restart + type(solver_workspace_${s}$), intent(inout) :: workspace + end subroutine + #:endfor + end interface + public :: solve_pccg_generic + + interface solve_pccg + #:for matrix in MATRIX_TYPES + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + #:if matrix == "dense" + ${t}$, intent(in) :: A(:,:) + #:else + type(${matrix}$_${s}$_type), intent(in) :: A + #:endif + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$, intent(in), optional :: tol + logical(1), intent(in), optional, target :: di(:) + integer, intent(in), optional :: maxiter + logical, intent(in), optional :: restart + type(solver_workspace_${s}$), optional, intent(inout), target :: workspace + end subroutine + #:endfor + #:endfor + end interface + public :: solve_pccg end module stdlib_linalg_iterative_solvers diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 4770076f8..66c226eb8 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -57,14 +57,16 @@ contains #:else type(${matrix}$_${s}$_type), intent(in) :: A #:endif - ${t}$, intent(in) :: b(:), tol + ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) + ${t}$, intent(in), optional :: tol integer, intent(in), optional :: maxiter type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !------------------------- type(linop_${s}$) :: op type(solver_workspace_${s}$), pointer :: workspace_ integer :: n, maxiter_ + ${t}$ :: tol_ !------------------------- n = size(b) op%matvec => default_matvec @@ -72,7 +74,11 @@ contains maxiter_ = n if(present(maxiter)) maxiter_ = maxiter + tol_ = 1.e-4_${s}$ + if(present(tol)) tol_ = tol + if(present(workspace)) then + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,3) ) workspace_ => workspace else allocate( workspace_%tmp(n,3) ) diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp new file mode 100644 index 000000000..e71653973 --- /dev/null +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -0,0 +1,161 @@ +#: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 MATRIX_TYPES = ["dense", "CSR"] +#:set RANKS = range(1, 2+1) + +submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg + use stdlib_kinds + use stdlib_sparse + use stdlib_constants + use stdlib_intrinsics, only: dot_product => stdlib_dot_product + use stdlib_linalg_iterative_solvers + implicit none + +contains + + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_pccg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + class(linop_${s}$), intent(in) :: A + ${t}$, intent(in) :: b(:), tol + ${t}$, intent(inout) :: x(:) + logical(1), intent(in) :: di(:) + integer, intent(in) :: maxiter + logical, intent(in) :: restart + type(solver_workspace_${s}$), intent(inout) :: workspace + !------------------------- + integer :: iter + ${t}$ :: norm_sq, norm_sq0, norm_sq_old, residual0 + ${t}$ :: zr1, zr2, zv2, alpha, beta, tolsq + !------------------------- + associate( R => workspace%tmp(:,1), & + S => workspace%tmp(:,2), & + P => workspace%tmp(:,3), & + Q => workspace%tmp(:,4)) + norm_sq = A%inner_product( b, b ) + norm_sq0 = norm_sq + residual0 = sqrt(norm_sq0) + if ( norm_sq0 > zero_${s}$ ) then + if(restart) X = zero_${s}$ + where( di ) X = B + + call A%matvec(X, R) + R = B - R + where( di ) R = zero_${s}$ + + call A%preconditionner(R,P) + where( di ) P = zero_${s}$ + + tolsq = tol*tol + iter = 0 + zr1 = zero_${s}$ + zr2 = one_${s}$ + do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) ) + + call A%preconditionner(R,S) + where ( di ) S = zero_${s}$ + zr2 = A%inner_product( R, S ) + + if (iter>0) then + beta = zr2 / zr1 + P = S + beta * P + end if + + call A%matvec(P, Q) + where( di ) Q = zero_${s}$ + zv2 = A%inner_product( P, Q ) + + alpha = zr2 / zv2 + + X = X + alpha * P + R = R - alpha * Q + norm_sq = A%inner_product( R, R ) + norm_sq_old = norm_sq + zr1 = zr2 + iter = iter + 1 + end do + end if + end associate + end subroutine + #:endfor + + #:for matrix in MATRIX_TYPES + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + #:if matrix == "dense" + ${t}$, intent(in) :: A(:,:) + #:else + type(${matrix}$_${s}$_type), intent(in) :: A + #:endif + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$, intent(in), optional :: tol + logical(1), intent(in), optional, target :: di(:) + integer, intent(in), optional :: maxiter + logical, intent(in), optional :: restart + type(solver_workspace_${s}$), optional, intent(inout), target :: workspace + !------------------------- + type(linop_${s}$) :: op + type(solver_workspace_${s}$), pointer :: workspace_ + integer :: n, maxiter_ + ${t}$ :: tol_ + logical :: restart_ + logical(1), pointer :: di_(:) + !------------------------- + n = size(b) + + op%matvec => default_matvec + op%inner_product => default_dot + op%preconditionner => default_preconditionner + + maxiter_ = n + if(present(maxiter)) maxiter_ = maxiter + restart_ = .true. + if(present(restart)) restart_ = restart + tol_ = 1.e-4_${s}$ + if(present(tol)) tol_ = tol + + if(present(di))then + di_ => di + else + allocate(di_(n),source=.false._1) + end if + + if(present(workspace)) then + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) ) + workspace_ => workspace + else + allocate( workspace_ ) + allocate( workspace_%tmp(n,4) ) + end if + + call solve_pccg_generic(op,b,x,di_,tol,maxiter_,restart_,workspace_) + + if(.not.present(di)) deallocate(di_) + contains + + subroutine default_matvec(x,y) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + #:if matrix == "dense" + y = matmul(A,x) + #:else + call spmv( A , x, y ) + #:endif + end subroutine + subroutine default_preconditionner(x,y) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + y = x + end subroutine + pure ${t}$ function default_dot(x,y) result(r) + ${t}$, intent(in) :: x(:) + ${t}$, intent(in) :: y(:) + r = dot_product(x,y) + end function + end subroutine + + #:endfor + #:endfor + +end submodule stdlib_linalg_iterative_pccg \ No newline at end of file From e551a5d052a821f4f6943ceb50cb4f7fe6ea6786 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Mon, 21 Apr 2025 21:26:05 +0200 Subject: [PATCH 04/26] complete cg with dirichlet flag, add example, fix di filter --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_solve_cg.f90 | 17 ++++ example/linalg/example_solve_pccg.f90 | 4 +- src/stdlib_linalg_iterative_solvers.fypp | 8 +- src/stdlib_linalg_iterative_solvers_cg.fypp | 99 ++++++++++++++----- src/stdlib_linalg_iterative_solvers_pccg.fypp | 52 ++++++---- 6 files changed, 130 insertions(+), 51 deletions(-) create mode 100644 example/linalg/example_solve_cg.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index d2db30042..1621919cc 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -41,6 +41,7 @@ ADD_EXAMPLE(get_norm) ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) +ADD_EXAMPLE(solve_cg) ADD_EXAMPLE(solve_pccg) ADD_EXAMPLE(sparse_from_ijv) ADD_EXAMPLE(sparse_data_accessors) diff --git a/example/linalg/example_solve_cg.f90 b/example/linalg/example_solve_cg.f90 new file mode 100644 index 000000000..626cf1184 --- /dev/null +++ b/example/linalg/example_solve_cg.f90 @@ -0,0 +1,17 @@ +program example_solve_pccg + use stdlib_kinds, only: dp + use stdlib_linalg_iterative_solvers, only: solve_cg + + real(dp) :: matrix(2,2) + real(dp) :: x(2), load(2) + + matrix = reshape( [4, 1,& + 1, 3] , [2,2]) + + x = dble( [2,1] ) !> initial guess + load = dble( [1,2] ) !> load vector + + call solve_cg(matrix, load, x, restart=.false.) + print *, x !> solution: [0.0909, 0.6364] + +end program \ No newline at end of file diff --git a/example/linalg/example_solve_pccg.f90 b/example/linalg/example_solve_pccg.f90 index 2826e9349..6c5b47b40 100644 --- a/example/linalg/example_solve_pccg.f90 +++ b/example/linalg/example_solve_pccg.f90 @@ -24,9 +24,9 @@ program example_solve_pccg dirichlet([1,5]) = .true._1 call solve_pccg(laplacian, load, x, tol=1.d-6, di=dirichlet) - print *, x + print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0] x = 0._dp call solve_pccg(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) - print *, x + print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0] end program \ No newline at end of file diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index ab196dfe4..374ea9a93 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -42,12 +42,14 @@ module stdlib_linalg_iterative_solvers interface solve_cg_generic #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) + module subroutine solve_cg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) class(linop_${s}$), intent(in) :: A ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$, intent(in) :: tol + logical(1), intent(in) :: di(:) integer, intent(in) :: maxiter + logical, intent(in) :: restart type(solver_workspace_${s}$), intent(inout) :: workspace end subroutine #:endfor @@ -57,7 +59,7 @@ module stdlib_linalg_iterative_solvers interface solve_cg #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_${matrix}$_${s}$(A,b,x,tol,maxiter,workspace) + module subroutine solve_cg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -66,7 +68,9 @@ module stdlib_linalg_iterative_solvers ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$, intent(in), optional :: tol + logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter + logical, intent(in), optional :: restart type(solver_workspace_${s}$), optional, intent(inout), target :: workspace end subroutine #:endfor diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 66c226eb8..6fce75903 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -15,43 +15,63 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) + module subroutine solve_cg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) class(linop_${s}$), intent(in) :: A ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) + logical(1), intent(in) :: di(:) integer, intent(in) :: maxiter + logical, intent(in) :: restart type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter - ${t}$ :: rtr, rtrold, alpha, beta, norm0_sq + ${t}$ :: norm_sq, norm_sq_old, norm_sq0, residual0, residual + ${t}$ :: alpha, beta, tolsq !------------------------- - associate( p => workspace%tmp(:,1), & - r => workspace%tmp(:,2), & + associate( P => workspace%tmp(:,1), & + R => workspace%tmp(:,2), & Ap => workspace%tmp(:,3)) - x = zero_${s}$ - rtr = A%inner_product(r, r) - norm0_sq = A%inner_product(b, b) - p = b + + norm_sq0 = A%inner_product(B, B) + residual0 = sqrt(norm_sq0) + + if(restart) X = zero_${s}$ + X = merge( B, X, di ) !> copy dirichlet load conditions encoded in B and indicated by di + + call A%matvec(X, R) + R = merge( zero_${s}$, B - R , di ) !> R = B - A*X + norm_sq = A%inner_product(R, R) + + P = R + + tolsq = tol*tol beta = zero_${s}$ - iter = 1 - do while( rtr > tol**2 * norm0_sq .and. iter < maxiter) - p = r + beta * p - call A%matvec(p,Ap) - alpha = rtr / A%inner_product(p, Ap) - x = x + alpha * p - r = r - alpha * Ap - rtrold = rtr - rtr = A%inner_product(r, r) - beta = rtr / rtrold + iter = 0 + do while( norm_sq > tolsq * norm_sq0 .and. iter < maxiter) + call A%matvec(P,Ap) + Ap = merge( zero_${s}$, Ap, di ) + + alpha = norm_sq / A%inner_product(P, Ap) + + X = X + alpha * P + R = R - alpha * Ap + + norm_sq_old = norm_sq + norm_sq = A%inner_product(R, R) + beta = norm_sq / norm_sq_old + + P = R + beta * P + iter = iter + 1 end do + residual = sqrt(norm_sq) end associate end subroutine #:endfor #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_${matrix}$_${s}$(A,b,x,tol,maxiter,workspace) + module subroutine solve_cg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -60,30 +80,57 @@ contains ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$, intent(in), optional :: tol + logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter + logical, intent(in), optional :: restart type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !------------------------- type(linop_${s}$) :: op type(solver_workspace_${s}$), pointer :: workspace_ integer :: n, maxiter_ ${t}$ :: tol_ + logical :: restart_ + logical(1), pointer :: di_(:) !------------------------- n = size(b) + maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter + restart_ = .true.; if(present(restart)) restart_ = restart + tol_ = 1.e-4_${s}$; if(present(tol)) tol_ = tol + + !------------------------- + ! internal memory setup op%matvec => default_matvec op%inner_product => default_dot - - maxiter_ = n - if(present(maxiter)) maxiter_ = maxiter - tol_ = 1.e-4_${s}$ - if(present(tol)) tol_ = tol + if(present(di))then + di_ => di + else + allocate(di_(n),source=.false._1) + end if if(present(workspace)) then if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,3) ) workspace_ => workspace else - allocate( workspace_%tmp(n,3) ) + allocate( workspace_ ) + allocate( workspace_%tmp(n,3), source = zero_${s}$ ) + end if + !------------------------- + ! main call to the solver + call solve_cg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) + + !------------------------- + ! internal memory cleanup + if(present(di))then + di_ => null() + else + deallocate(di_) + end if + if(present(workspace)) then + workspace_ => null() + else + deallocate( workspace_%tmp ) + deallocate( workspace_ ) end if - call solve_cg_generic(op,b,x,tol,maxiter_,workspace_) contains subroutine default_matvec(x,y) diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index e71653973..23e71d436 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -25,7 +25,7 @@ contains type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter - ${t}$ :: norm_sq, norm_sq0, norm_sq_old, residual0 + ${t}$ :: norm_sq, norm_sq0, norm_sq_old, residual0, residual ${t}$ :: zr1, zr2, zv2, alpha, beta, tolsq !------------------------- associate( R => workspace%tmp(:,1), & @@ -37,14 +37,13 @@ contains residual0 = sqrt(norm_sq0) if ( norm_sq0 > zero_${s}$ ) then if(restart) X = zero_${s}$ - where( di ) X = B + X = merge( B, X, di ) !> copy dirichlet load conditions encoded in B and indicated by di call A%matvec(X, R) - R = B - R - where( di ) R = zero_${s}$ + R = merge( zero_${s}$, B - R , di ) !> R = B - A*X - call A%preconditionner(R,P) - where( di ) P = zero_${s}$ + call A%preconditionner(R,P) !> P = M^{-1}*R + P = merge( zero_${s}$, P, di ) tolsq = tol*tol iter = 0 @@ -52,8 +51,8 @@ contains zr2 = one_${s}$ do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) ) - call A%preconditionner(R,S) - where ( di ) S = zero_${s}$ + call A%preconditionner(R,S) !> S = M^{-1}*R + S = merge( zero_${s}$, S, di ) zr2 = A%inner_product( R, S ) if (iter>0) then @@ -62,7 +61,7 @@ contains end if call A%matvec(P, Q) - where( di ) Q = zero_${s}$ + Q = merge( zero_${s}$, Q, di ) zv2 = A%inner_product( P, Q ) alpha = zr2 / zv2 @@ -75,6 +74,7 @@ contains iter = iter + 1 end do end if + residual = sqrt(norm_sq) end associate end subroutine #:endfor @@ -103,18 +103,15 @@ contains logical(1), pointer :: di_(:) !------------------------- n = size(b) + maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter + restart_ = .true.; if(present(restart)) restart_ = restart + tol_ = 1.e-4_${s}$; if(present(tol)) tol_ = tol + !------------------------- + ! internal memory setup op%matvec => default_matvec op%inner_product => default_dot op%preconditionner => default_preconditionner - - maxiter_ = n - if(present(maxiter)) maxiter_ = maxiter - restart_ = .true. - if(present(restart)) restart_ = restart - tol_ = 1.e-4_${s}$ - if(present(tol)) tol_ = tol - if(present(di))then di_ => di else @@ -126,12 +123,25 @@ contains workspace_ => workspace else allocate( workspace_ ) - allocate( workspace_%tmp(n,4) ) + allocate( workspace_%tmp(n,4) , source = zero_${s}$ ) end if - - call solve_pccg_generic(op,b,x,di_,tol,maxiter_,restart_,workspace_) + !------------------------- + ! main call to the solver + call solve_pccg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) - if(.not.present(di)) deallocate(di_) + !------------------------- + ! internal memory cleanup + if(present(di))then + di_ => null() + else + deallocate(di_) + end if + if(present(workspace)) then + workspace_ => null() + else + deallocate( workspace_%tmp ) + deallocate( workspace_ ) + end if contains subroutine default_matvec(x,y) From 19167d5e440695c850e4a4918a09a137bbaa82c1 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Wed, 23 Apr 2025 21:27:17 +0200 Subject: [PATCH 05/26] add other sparse matrices --- src/stdlib_linalg_iterative_solvers.fypp | 2 +- src/stdlib_linalg_iterative_solvers_cg.fypp | 2 +- src/stdlib_linalg_iterative_solvers_pccg.fypp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 374ea9a93..1f2916cc2 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -1,7 +1,7 @@ #: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 MATRIX_TYPES = ["dense", "CSR"] +#:set MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"] #:set RANKS = range(1, 2+1) module stdlib_linalg_iterative_solvers diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 6fce75903..d87473906 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -1,7 +1,7 @@ #: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 MATRIX_TYPES = ["dense", "CSR"] +#:set MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"] #:set RANKS = range(1, 2+1) submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 23e71d436..51ea18961 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -1,7 +1,7 @@ #: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 MATRIX_TYPES = ["dense", "CSR"] +#:set MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"] #:set RANKS = range(1, 2+1) submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg From 93249712adaf84e0667020b9a863a76afbff21e2 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Wed, 23 Apr 2025 21:29:12 +0200 Subject: [PATCH 06/26] add example for custom solver extending the generic interface --- example/linalg/CMakeLists.txt | 1 + example/linalg/example_solve_custom.f90 | 116 ++++++++++++++++++++++++ 2 files changed, 117 insertions(+) create mode 100644 example/linalg/example_solve_custom.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 1621919cc..80dd31fb7 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -43,6 +43,7 @@ ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) ADD_EXAMPLE(solve_cg) ADD_EXAMPLE(solve_pccg) +ADD_EXAMPLE(solve_custom) ADD_EXAMPLE(sparse_from_ijv) ADD_EXAMPLE(sparse_data_accessors) ADD_EXAMPLE(sparse_spmv) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 new file mode 100644 index 000000000..26014ad02 --- /dev/null +++ b/example/linalg/example_solve_custom.f90 @@ -0,0 +1,116 @@ +module custom_solver + use stdlib_kinds, only: dp + use stdlib_sparse + use stdlib_linalg_iterative_solvers, only: linop_dp, solver_workspace_dp, solve_pccg_generic + implicit none +contains + subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) + type(CSR_dp_type), intent(in) :: A + real(dp), intent(in) :: b(:) + real(dp), intent(inout) :: x(:) + real(dp), intent(in), optional :: tol + logical(1), intent(in), optional, target :: di(:) + integer, intent(in), optional :: maxiter + logical, intent(in), optional :: restart + type(solver_workspace_dp), optional, intent(inout), target :: workspace + !------------------------- + type(linop_dp) :: op + type(solver_workspace_dp), pointer :: workspace_ + integer :: n, maxiter_ + real(dp) :: tol_ + logical :: restart_ + logical(1), pointer :: di_(:) + real(dp), allocatable :: diagonal(:) + !------------------------- + n = size(b) + maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter + restart_ = .true.; if(present(restart)) restart_ = restart + tol_ = 1.e-4_dp; if(present(tol)) tol_ = tol + + !------------------------- + ! internal memory setup + op%matvec => my_matvec + op%inner_product => my_dot + op%preconditionner => jacobi_preconditionner + if(present(di))then + di_ => di + else + allocate(di_(n),source=.false._1) + end if + + if(present(workspace)) then + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) ) + workspace_ => workspace + else + allocate( workspace_ ) + allocate( workspace_%tmp(n,4) , source = 0._dp ) + end if + call diag(A,diagonal) + where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal + !------------------------- + ! main call to the solver + call solve_pccg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) + + !------------------------- + ! internal memory cleanup + if(present(di))then + di_ => null() + else + deallocate(di_) + end if + if(present(workspace)) then + workspace_ => null() + else + deallocate( workspace_%tmp ) + deallocate( workspace_ ) + end if + contains + + subroutine my_matvec(x,y) + real(dp), intent(in) :: x(:) + real(dp), intent(inout) :: y(:) + call spmv( A , x, y ) + end subroutine + subroutine jacobi_preconditionner(x,y) + real(dp), intent(in) :: x(:) + real(dp), intent(inout) :: y(:) + y = diagonal * x + end subroutine + pure real(dp) function my_dot(x,y) result(r) + real(dp), intent(in) :: x(:) + real(dp), intent(in) :: y(:) + r = dot_product(x,y) + end function + end subroutine + +end module custom_solver + + +program example_solve_custom + use custom_solver + implicit none + + type(CSR_dp_type) :: laplacian_csr + type(COO_dp_type) :: COO + real(dp) :: laplacian(5,5) + real(dp) :: x(5), load(5) + logical(1) :: dirichlet(5) + + laplacian = reshape( [1, -1, 0, 0, 0,& + -1, 2, -1, 0, 0,& + 0, -1, 2, -1, 0,& + 0, 0, -1, 2, -1,& + 0, 0, 0, -1, 1] , [5,5]) + call dense2coo(laplacian,COO) + call coo2csr(COO,laplacian_csr) + + x = 0._dp + load = dble( [0,0,5,0,0] ) + + dirichlet = .false._1 + dirichlet([1,5]) = .true._1 + + call solve_pccg_custom(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) + print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0] + +end program example_solve_custom \ No newline at end of file From 85a70bacae1de90a692d325a6394c2500a088708 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Fri, 25 Apr 2025 17:17:56 +0200 Subject: [PATCH 07/26] small simplifications for working data --- example/linalg/example_solve_custom.f90 | 18 ++++++++---------- src/stdlib_linalg_iterative_solvers_cg.fypp | 16 ++++++---------- src/stdlib_linalg_iterative_solvers_pccg.fypp | 16 ++++++---------- 3 files changed, 20 insertions(+), 30 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index 26014ad02..74c224bab 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -39,12 +39,13 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) end if if(present(workspace)) then - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) ) workspace_ => workspace else allocate( workspace_ ) - allocate( workspace_%tmp(n,4) , source = 0._dp ) end if + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) , source = 0.d0 ) + !------------------------- + ! Jacobi preconditionner factorization call diag(A,diagonal) where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal !------------------------- @@ -53,17 +54,14 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) !------------------------- ! internal memory cleanup - if(present(di))then - di_ => null() - else - deallocate(di_) - end if - if(present(workspace)) then - workspace_ => null() - else + if(.not.present(di)) deallocate(di_) + di_ => null() + + if(.not.present(workspace)) then deallocate( workspace_%tmp ) deallocate( workspace_ ) end if + workspace_ => null() contains subroutine my_matvec(x,y) diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index d87473906..76c91077e 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -108,29 +108,25 @@ contains end if if(present(workspace)) then - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,3) ) workspace_ => workspace else allocate( workspace_ ) - allocate( workspace_%tmp(n,3), source = zero_${s}$ ) end if + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,3), source = zero_${s}$ ) !------------------------- ! main call to the solver call solve_cg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) !------------------------- ! internal memory cleanup - if(present(di))then - di_ => null() - else - deallocate(di_) - end if - if(present(workspace)) then - workspace_ => null() - else + if(.not.present(di)) deallocate(di_) + di_ => null() + + if(.not.present(workspace)) then deallocate( workspace_%tmp ) deallocate( workspace_ ) end if + workspace_ => null() contains subroutine default_matvec(x,y) diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 51ea18961..fd76191c7 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -119,29 +119,25 @@ contains end if if(present(workspace)) then - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) ) workspace_ => workspace else allocate( workspace_ ) - allocate( workspace_%tmp(n,4) , source = zero_${s}$ ) end if + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) , source = zero_${s}$ ) !------------------------- ! main call to the solver call solve_pccg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) !------------------------- ! internal memory cleanup - if(present(di))then - di_ => null() - else - deallocate(di_) - end if - if(present(workspace)) then - workspace_ => null() - else + if(.not.present(di)) deallocate(di_) + di_ => null() + + if(.not.present(workspace)) then deallocate( workspace_%tmp ) deallocate( workspace_ ) end if + workspace_ => null() contains subroutine default_matvec(x,y) From 84f4bc9fcd08b4c34ed8983dbb310c7e526228a9 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Fri, 25 Apr 2025 20:46:29 +0200 Subject: [PATCH 08/26] make default inner_product point to a default dot_product add enum for workspace sizes --- src/stdlib_linalg_iterative_solvers.fypp | 26 +++++++++++++++++-- src/stdlib_linalg_iterative_solvers_cg.fypp | 10 ++----- src/stdlib_linalg_iterative_solvers_pccg.fypp | 10 ++----- 3 files changed, 28 insertions(+), 18 deletions(-) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 1f2916cc2..ca267442b 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -10,11 +10,17 @@ module stdlib_linalg_iterative_solvers implicit none private + enum, bind(c) + enumerator :: size_wksp_cg = 3 + enumerator :: size_wksp_pccg = 4 + end enum + public :: size_wksp_cg, size_wksp_pccg + #:for k, t, s in R_KINDS_TYPES type, public :: linop_${s}$ procedure(vector_sub_${s}$), nopass, pointer :: matvec => null() procedure(vector_sub_${s}$), nopass, pointer :: preconditionner => null() - procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => null() + procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => default_dot_${s}$ end type #:endfor @@ -114,6 +120,22 @@ module stdlib_linalg_iterative_solvers #:endfor #:endfor end interface - public :: solve_pccg + public :: solve_pccg + + +contains + + !------------------------------------------------------------------ + ! defaults + !------------------------------------------------------------------ + #:for k, t, s in R_KINDS_TYPES + pure ${t}$ function default_dot_${s}$(x,y) result(r) + use stdlib_intrinsics, only: stdlib_dot_product + ${t}$, intent(in) :: x(:) + ${t}$, intent(in) :: y(:) + r = stdlib_dot_product(x,y) + end function + + #:endfor end module stdlib_linalg_iterative_solvers diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 76c91077e..0baa27487 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -8,7 +8,6 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg use stdlib_kinds use stdlib_sparse use stdlib_constants - use stdlib_intrinsics, only: dot_product => stdlib_dot_product use stdlib_linalg_iterative_solvers implicit none @@ -100,7 +99,7 @@ contains !------------------------- ! internal memory setup op%matvec => default_matvec - op%inner_product => default_dot + ! op%inner_product => default_dot if(present(di))then di_ => di else @@ -112,7 +111,7 @@ contains else allocate( workspace_ ) end if - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,3), source = zero_${s}$ ) + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_cg), source = zero_${s}$ ) !------------------------- ! main call to the solver call solve_cg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) @@ -138,11 +137,6 @@ contains call spmv( A , x, y ) #:endif end subroutine - pure ${t}$ function default_dot(x,y) result(r) - ${t}$, intent(in) :: x(:) - ${t}$, intent(in) :: y(:) - r = dot_product(x,y) - end function end subroutine #:endfor diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index fd76191c7..85c861c8c 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -8,7 +8,6 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg use stdlib_kinds use stdlib_sparse use stdlib_constants - use stdlib_intrinsics, only: dot_product => stdlib_dot_product use stdlib_linalg_iterative_solvers implicit none @@ -110,7 +109,7 @@ contains !------------------------- ! internal memory setup op%matvec => default_matvec - op%inner_product => default_dot + ! op%inner_product => default_dot op%preconditionner => default_preconditionner if(present(di))then di_ => di @@ -123,7 +122,7 @@ contains else allocate( workspace_ ) end if - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) , source = zero_${s}$ ) + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = zero_${s}$ ) !------------------------- ! main call to the solver call solve_pccg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) @@ -154,11 +153,6 @@ contains ${t}$, intent(inout) :: y(:) y = x end subroutine - pure ${t}$ function default_dot(x,y) result(r) - ${t}$, intent(in) :: x(:) - ${t}$, intent(in) :: y(:) - r = dot_product(x,y) - end function end subroutine #:endfor From 3ec23a4e932078fcb1688a762755c66516d5b828 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Fri, 25 Apr 2025 21:02:36 +0200 Subject: [PATCH 09/26] make preconditionner a linop --- example/linalg/example_solve_custom.f90 | 5 +++-- src/stdlib_linalg_iterative_solvers.fypp | 7 +++--- src/stdlib_linalg_iterative_solvers_pccg.fypp | 22 +++++++++++++------ 3 files changed, 22 insertions(+), 12 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index 74c224bab..cfbc6fee6 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -15,6 +15,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) type(solver_workspace_dp), optional, intent(inout), target :: workspace !------------------------- type(linop_dp) :: op + type(linop_dp) :: M type(solver_workspace_dp), pointer :: workspace_ integer :: n, maxiter_ real(dp) :: tol_ @@ -31,7 +32,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) ! internal memory setup op%matvec => my_matvec op%inner_product => my_dot - op%preconditionner => jacobi_preconditionner + M%matvec => jacobi_preconditionner if(present(di))then di_ => di else @@ -50,7 +51,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal !------------------------- ! main call to the solver - call solve_pccg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) + call solve_pccg_generic(op,M,b,x,di_,tol_,maxiter_,restart_,workspace_) !------------------------- ! internal memory cleanup diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index ca267442b..3d8765f76 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -19,7 +19,6 @@ module stdlib_linalg_iterative_solvers #:for k, t, s in R_KINDS_TYPES type, public :: linop_${s}$ procedure(vector_sub_${s}$), nopass, pointer :: matvec => null() - procedure(vector_sub_${s}$), nopass, pointer :: preconditionner => null() procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => default_dot_${s}$ end type #:endfor @@ -86,8 +85,9 @@ module stdlib_linalg_iterative_solvers interface solve_pccg_generic #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_pccg_generic_${s}$(A,M,b,x,di,tol,maxiter,restart,workspace) class(linop_${s}$), intent(in) :: A + class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$, intent(in) :: tol @@ -103,7 +103,7 @@ module stdlib_linalg_iterative_solvers interface solve_pccg #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,M,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -115,6 +115,7 @@ module stdlib_linalg_iterative_solvers logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter logical, intent(in), optional :: restart + class(linop_${s}$), optional , intent(in), target :: M !> preconditionner type(solver_workspace_${s}$), optional, intent(inout), target :: workspace end subroutine #:endfor diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 85c861c8c..2c12f1f6b 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -14,8 +14,9 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_pccg_generic_${s}$(A,M,b,x,di,tol,maxiter,restart,workspace) class(linop_${s}$), intent(in) :: A + class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) logical(1), intent(in) :: di(:) @@ -41,7 +42,7 @@ contains call A%matvec(X, R) R = merge( zero_${s}$, B - R , di ) !> R = B - A*X - call A%preconditionner(R,P) !> P = M^{-1}*R + call M%matvec(R,P) !> P = M^{-1}*R P = merge( zero_${s}$, P, di ) tolsq = tol*tol @@ -50,7 +51,7 @@ contains zr2 = one_${s}$ do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) ) - call A%preconditionner(R,S) !> S = M^{-1}*R + call M%matvec(R,S) !> S = M^{-1}*R S = merge( zero_${s}$, S, di ) zr2 = A%inner_product( R, S ) @@ -80,7 +81,7 @@ contains #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,M,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -92,9 +93,11 @@ contains logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter logical, intent(in), optional :: restart + class(linop_${s}$), optional , intent(in), target :: M !> preconditionner type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !------------------------- type(linop_${s}$) :: op + type(linop_${s}$), pointer :: M_ => null() type(solver_workspace_${s}$), pointer :: workspace_ integer :: n, maxiter_ ${t}$ :: tol_ @@ -109,8 +112,12 @@ contains !------------------------- ! internal memory setup op%matvec => default_matvec - ! op%inner_product => default_dot - op%preconditionner => default_preconditionner + if(present(M)) then + M_ => M + else + allocate( M_ ) + M_%matvec => default_preconditionner + end if if(present(di))then di_ => di else @@ -125,7 +132,7 @@ contains if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = zero_${s}$ ) !------------------------- ! main call to the solver - call solve_pccg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) + call solve_pccg_generic(op,M_,b,x,di_,tol_,maxiter_,restart_,workspace_) !------------------------- ! internal memory cleanup @@ -136,6 +143,7 @@ contains deallocate( workspace_%tmp ) deallocate( workspace_ ) end if + M_ => null() workspace_ => null() contains From bfafaa5a77e00ceb2074d58dc4f30d04aa3cbd60 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Fri, 25 Apr 2025 21:07:21 +0200 Subject: [PATCH 10/26] use facility size --- example/linalg/example_solve_custom.f90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index cfbc6fee6..444a20190 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -1,7 +1,10 @@ module custom_solver use stdlib_kinds, only: dp use stdlib_sparse - use stdlib_linalg_iterative_solvers, only: linop_dp, solver_workspace_dp, solve_pccg_generic + use stdlib_linalg_iterative_solvers, only: linop_dp, & + solver_workspace_dp, & + solve_pccg_generic, & + size_wksp_pccg implicit none contains subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) @@ -44,7 +47,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) else allocate( workspace_ ) end if - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,4) , source = 0.d0 ) + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = 0.d0 ) !------------------------- ! Jacobi preconditionner factorization call diag(A,diagonal) From 0b01dbd07d65da99a0221f8d832ff0e329e481b2 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sat, 26 Apr 2025 13:17:02 +0200 Subject: [PATCH 11/26] Add a customizable logger facility, change linop matvec to apply --- example/linalg/example_solve_custom.f90 | 15 ++++++++++--- src/stdlib_linalg_iterative_solvers.fypp | 12 +++++++++- src/stdlib_linalg_iterative_solvers_cg.fypp | 16 ++++++++------ src/stdlib_linalg_iterative_solvers_pccg.fypp | 22 ++++++++++--------- 4 files changed, 44 insertions(+), 21 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index 444a20190..5692f3eeb 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -25,17 +25,18 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) logical :: restart_ logical(1), pointer :: di_(:) real(dp), allocatable :: diagonal(:) + real(dp) :: norm_sq0 !------------------------- n = size(b) maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter restart_ = .true.; if(present(restart)) restart_ = restart tol_ = 1.e-4_dp; if(present(tol)) tol_ = tol - + norm_sq0 = 0.d0 !------------------------- ! internal memory setup - op%matvec => my_matvec + op%apply => my_matvec op%inner_product => my_dot - M%matvec => jacobi_preconditionner + M%apply => jacobi_preconditionner if(present(di))then di_ => di else @@ -48,6 +49,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) allocate( workspace_ ) end if if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = 0.d0 ) + workspace_%callback => my_logger !------------------------- ! Jacobi preconditionner factorization call diag(A,diagonal) @@ -83,6 +85,13 @@ pure real(dp) function my_dot(x,y) result(r) real(dp), intent(in) :: y(:) r = dot_product(x,y) end function + subroutine my_logger(x,norm_sq,iter) + real(dp), intent(in) :: x(:) + real(dp), intent(in) :: norm_sq + integer, intent(in) :: iter + if(iter == 0) norm_sq0 = norm_sq + print *, "Iteration: ", iter, " Residual: ", sqrt(norm_sq), " Relative: ", sqrt(norm_sq)/sqrt(norm_sq0) + end subroutine end subroutine end module custom_solver diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 3d8765f76..792355eea 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -10,15 +10,18 @@ module stdlib_linalg_iterative_solvers implicit none private + !> brief workspace size for the iterative solvers + !> details The size of the workspace is defined by the number of vectors used in the iterative solver. enum, bind(c) enumerator :: size_wksp_cg = 3 enumerator :: size_wksp_pccg = 4 end enum public :: size_wksp_cg, size_wksp_pccg + !> linear operator class for the iterative solvers #:for k, t, s in R_KINDS_TYPES type, public :: linop_${s}$ - procedure(vector_sub_${s}$), nopass, pointer :: matvec => null() + procedure(vector_sub_${s}$), nopass, pointer :: apply => null() procedure(reduction_sub_${s}$), nopass, pointer :: inner_product => default_dot_${s}$ end type #:endfor @@ -26,6 +29,7 @@ module stdlib_linalg_iterative_solvers #:for k, t, s in R_KINDS_TYPES type, public :: solver_workspace_${s}$ ${t}$, allocatable :: tmp(:,:) + procedure(logger_sub_${s}$), pointer, nopass :: callback => null() end type #:endfor @@ -42,6 +46,12 @@ module stdlib_linalg_iterative_solvers ${t}$, intent(in) :: x(:) ${t}$, intent(in) :: y(:) end function + subroutine logger_sub_${s}$(x,norm_sq,iter) + import :: ${k}$ + ${t}$, intent(in) :: x(:) + ${t}$, intent(in) :: norm_sq + integer, intent(in) :: iter + end subroutine #:endfor end interface diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 0baa27487..73245fa6e 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -24,20 +24,21 @@ contains type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter - ${t}$ :: norm_sq, norm_sq_old, norm_sq0, residual0, residual + ${t}$ :: norm_sq, norm_sq_old, norm_sq0 ${t}$ :: alpha, beta, tolsq !------------------------- + iter = 0 associate( P => workspace%tmp(:,1), & R => workspace%tmp(:,2), & Ap => workspace%tmp(:,3)) norm_sq0 = A%inner_product(B, B) - residual0 = sqrt(norm_sq0) + if(associated(workspace%callback)) call workspace%callback(x, norm_sq0, iter) if(restart) X = zero_${s}$ X = merge( B, X, di ) !> copy dirichlet load conditions encoded in B and indicated by di - call A%matvec(X, R) + call A%apply(X, R) R = merge( zero_${s}$, B - R , di ) !> R = B - A*X norm_sq = A%inner_product(R, R) @@ -45,9 +46,9 @@ contains tolsq = tol*tol beta = zero_${s}$ - iter = 0 + if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) do while( norm_sq > tolsq * norm_sq0 .and. iter < maxiter) - call A%matvec(P,Ap) + call A%apply(P,Ap) Ap = merge( zero_${s}$, Ap, di ) alpha = norm_sq / A%inner_product(P, Ap) @@ -62,8 +63,9 @@ contains P = R + beta * P iter = iter + 1 + + if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) end do - residual = sqrt(norm_sq) end associate end subroutine #:endfor @@ -98,7 +100,7 @@ contains !------------------------- ! internal memory setup - op%matvec => default_matvec + op%apply => default_matvec ! op%inner_product => default_dot if(present(di))then di_ => di diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 2c12f1f6b..9ebe636a7 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -25,33 +25,35 @@ contains type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter - ${t}$ :: norm_sq, norm_sq0, norm_sq_old, residual0, residual + ${t}$ :: norm_sq, norm_sq0, norm_sq_old ${t}$ :: zr1, zr2, zv2, alpha, beta, tolsq !------------------------- + iter = 0 associate( R => workspace%tmp(:,1), & S => workspace%tmp(:,2), & P => workspace%tmp(:,3), & Q => workspace%tmp(:,4)) norm_sq = A%inner_product( b, b ) norm_sq0 = norm_sq - residual0 = sqrt(norm_sq0) + if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) + if ( norm_sq0 > zero_${s}$ ) then if(restart) X = zero_${s}$ X = merge( B, X, di ) !> copy dirichlet load conditions encoded in B and indicated by di - call A%matvec(X, R) + call A%apply(X, R) R = merge( zero_${s}$, B - R , di ) !> R = B - A*X - call M%matvec(R,P) !> P = M^{-1}*R + call M%apply(R,P) !> P = M^{-1}*R P = merge( zero_${s}$, P, di ) tolsq = tol*tol - iter = 0 + zr1 = zero_${s}$ zr2 = one_${s}$ do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) ) - call M%matvec(R,S) !> S = M^{-1}*R + call M%apply(R,S) !> S = M^{-1}*R S = merge( zero_${s}$, S, di ) zr2 = A%inner_product( R, S ) @@ -60,7 +62,7 @@ contains P = S + beta * P end if - call A%matvec(P, Q) + call A%apply(P, Q) Q = merge( zero_${s}$, Q, di ) zv2 = A%inner_product( P, Q ) @@ -72,9 +74,9 @@ contains norm_sq_old = norm_sq zr1 = zr2 iter = iter + 1 + if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) end do end if - residual = sqrt(norm_sq) end associate end subroutine #:endfor @@ -111,12 +113,12 @@ contains !------------------------- ! internal memory setup - op%matvec => default_matvec + op%apply => default_matvec if(present(M)) then M_ => M else allocate( M_ ) - M_%matvec => default_preconditionner + M_%apply => default_preconditionner end if if(present(di))then di_ => di From 07b97cef061e2b345f0fbfdfc1022750c3261eab Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sat, 26 Apr 2025 13:34:48 +0200 Subject: [PATCH 12/26] change internal procedure names for custom example --- example/linalg/example_solve_custom.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index 5692f3eeb..ce1b58cb0 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -34,9 +34,9 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) norm_sq0 = 0.d0 !------------------------- ! internal memory setup - op%apply => my_matvec + op%apply => my_apply op%inner_product => my_dot - M%apply => jacobi_preconditionner + M%apply => my_jacobi_preconditionner if(present(di))then di_ => di else @@ -70,12 +70,12 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) workspace_ => null() contains - subroutine my_matvec(x,y) + subroutine my_apply(x,y) real(dp), intent(in) :: x(:) real(dp), intent(inout) :: y(:) call spmv( A , x, y ) end subroutine - subroutine jacobi_preconditionner(x,y) + subroutine my_jacobi_preconditionner(x,y) real(dp), intent(in) :: x(:) real(dp), intent(inout) :: y(:) y = diagonal * x From 379cd819acea3f9516ca61b3811db9947d90e6fd Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 27 Apr 2025 00:38:24 +0200 Subject: [PATCH 13/26] refactor to remove hard dependency on dirichlet BCs --- example/linalg/example_solve_custom.f90 | 15 ++++--- src/stdlib_linalg_iterative_solvers.fypp | 14 +++---- src/stdlib_linalg_iterative_solvers_cg.fypp | 27 ++++++------- src/stdlib_linalg_iterative_solvers_pccg.fypp | 40 ++++++++++--------- 4 files changed, 50 insertions(+), 46 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index ce1b58cb0..49faf750c 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -56,7 +56,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal !------------------------- ! main call to the solver - call solve_pccg_generic(op,M,b,x,di_,tol_,maxiter_,restart_,workspace_) + call solve_pccg_generic(op,M,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup @@ -70,15 +70,20 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) workspace_ => null() contains - subroutine my_apply(x,y) + subroutine my_apply(x,y,alpha,beta) real(dp), intent(in) :: x(:) real(dp), intent(inout) :: y(:) - call spmv( A , x, y ) + real(dp), intent(in) :: alpha + real(dp), intent(in) :: beta + call spmv( A , x, y , alpha, beta ) + y = merge( 0._dp, y, di_ ) end subroutine - subroutine my_jacobi_preconditionner(x,y) + subroutine my_jacobi_preconditionner(x,y,alpha,beta) real(dp), intent(in) :: x(:) real(dp), intent(inout) :: y(:) - y = diagonal * x + real(dp), intent(in) :: alpha + real(dp), intent(in) :: beta + y = merge( 0._dp, diagonal * x , di_ ) end subroutine pure real(dp) function my_dot(x,y) result(r) real(dp), intent(in) :: x(:) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 792355eea..701f2c5b4 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -36,10 +36,12 @@ module stdlib_linalg_iterative_solvers abstract interface #:for k, t, s in R_KINDS_TYPES - subroutine vector_sub_${s}$(x,y) + subroutine vector_sub_${s}$(x,y,alpha,beta) import :: ${k}$ ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta end subroutine pure ${t}$ function reduction_sub_${s}$(x,y) result(r) import :: ${k}$ @@ -57,14 +59,12 @@ module stdlib_linalg_iterative_solvers interface solve_cg_generic #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$, intent(in) :: tol - logical(1), intent(in) :: di(:) integer, intent(in) :: maxiter - logical, intent(in) :: restart type(solver_workspace_${s}$), intent(inout) :: workspace end subroutine #:endfor @@ -95,15 +95,13 @@ module stdlib_linalg_iterative_solvers interface solve_pccg_generic #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_generic_${s}$(A,M,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_pccg_generic_${s}$(A,M,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$, intent(in) :: tol - logical(1), intent(in) :: di(:) integer, intent(in) :: maxiter - logical, intent(in) :: restart type(solver_workspace_${s}$), intent(inout) :: workspace end subroutine #:endfor @@ -124,7 +122,7 @@ module stdlib_linalg_iterative_solvers ${t}$, intent(in), optional :: tol logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter - logical, intent(in), optional :: restart + logical, intent(in), optional :: restart class(linop_${s}$), optional , intent(in), target :: M !> preconditionner type(solver_workspace_${s}$), optional, intent(inout), target :: workspace end subroutine diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 73245fa6e..241403662 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -14,13 +14,11 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_generic_${s}$(A,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) - logical(1), intent(in) :: di(:) integer, intent(in) :: maxiter - logical, intent(in) :: restart type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter @@ -35,11 +33,8 @@ contains norm_sq0 = A%inner_product(B, B) if(associated(workspace%callback)) call workspace%callback(x, norm_sq0, iter) - if(restart) X = zero_${s}$ - X = merge( B, X, di ) !> copy dirichlet load conditions encoded in B and indicated by di - - call A%apply(X, R) - R = merge( zero_${s}$, B - R , di ) !> R = B - A*X + R = B + call A%apply(X, R, alpha= -one_${s}$, beta=one_${s}$) !> R = B - A*X norm_sq = A%inner_product(R, R) P = R @@ -48,8 +43,7 @@ contains beta = zero_${s}$ if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) do while( norm_sq > tolsq * norm_sq0 .and. iter < maxiter) - call A%apply(P,Ap) - Ap = merge( zero_${s}$, Ap, di ) + call A%apply(P,Ap, alpha= one_${s}$, beta=zero_${s}$) !> Ap = A*P alpha = norm_sq / A%inner_product(P, Ap) @@ -116,7 +110,9 @@ contains if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_cg), source = zero_${s}$ ) !------------------------- ! main call to the solver - call solve_cg_generic(op,b,x,di_,tol_,maxiter_,restart_,workspace_) + if(restart_) x = zero_${s}$ + x = merge( b, x, di_ ) !> copy dirichlet load conditions encoded in B and indicated by di + call solve_cg_generic(op,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup @@ -130,14 +126,17 @@ contains workspace_ => null() contains - subroutine default_matvec(x,y) + subroutine default_matvec(x,y,alpha,beta) ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta #:if matrix == "dense" - y = matmul(A,x) + y = alpha * matmul(A,x) + beta * y #:else - call spmv( A , x, y ) + call spmv( A , x, y , alpha, beta ) #:endif + y = merge( zero_${s}$, y, di_ ) end subroutine end subroutine diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 9ebe636a7..ec1b83d6a 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -14,14 +14,12 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_generic_${s}$(A,M,b,x,di,tol,maxiter,restart,workspace) + module subroutine solve_pccg_generic_${s}$(A,M,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) - logical(1), intent(in) :: di(:) integer, intent(in) :: maxiter - logical, intent(in) :: restart type(solver_workspace_${s}$), intent(inout) :: workspace !------------------------- integer :: iter @@ -38,14 +36,11 @@ contains if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) if ( norm_sq0 > zero_${s}$ ) then - if(restart) X = zero_${s}$ - X = merge( B, X, di ) !> copy dirichlet load conditions encoded in B and indicated by di - call A%apply(X, R) - R = merge( zero_${s}$, B - R , di ) !> R = B - A*X + R = B + call A%apply(X, R, alpha= -one_${s}$, beta=one_${s}$) !> R = B - A*X - call M%apply(R,P) !> P = M^{-1}*R - P = merge( zero_${s}$, P, di ) + call M%apply(R,P, alpha= one_${s}$, beta=zero_${s}$) !> P = M^{-1}*R tolsq = tol*tol @@ -53,8 +48,7 @@ contains zr2 = one_${s}$ do while ( (iter < maxiter) .AND. (norm_sq > tolsq * norm_sq0) ) - call M%apply(R,S) !> S = M^{-1}*R - S = merge( zero_${s}$, S, di ) + call M%apply(R,S, alpha= one_${s}$, beta=zero_${s}$) !> S = M^{-1}*R zr2 = A%inner_product( R, S ) if (iter>0) then @@ -62,16 +56,17 @@ contains P = S + beta * P end if - call A%apply(P, Q) - Q = merge( zero_${s}$, Q, di ) + call A%apply(P, Q, alpha= one_${s}$, beta=zero_${s}$) !> Q = A*P zv2 = A%inner_product( P, Q ) alpha = zr2 / zv2 X = X + alpha * P R = R - alpha * Q + norm_sq = A%inner_product( R, R ) norm_sq_old = norm_sq + zr1 = zr2 iter = iter + 1 if(associated(workspace%callback)) call workspace%callback(x, norm_sq, iter) @@ -134,7 +129,9 @@ contains if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = zero_${s}$ ) !------------------------- ! main call to the solver - call solve_pccg_generic(op,M_,b,x,di_,tol_,maxiter_,restart_,workspace_) + if(restart_) x = zero_${s}$ + x = merge( b, x, di_ ) !> copy dirichlet load conditions encoded in B and indicated by di + call solve_pccg_generic(op,M_,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup @@ -149,19 +146,24 @@ contains workspace_ => null() contains - subroutine default_matvec(x,y) + subroutine default_matvec(x,y,alpha,beta) ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta #:if matrix == "dense" - y = matmul(A,x) + y = alpha * matmul(A,x) + beta * y #:else - call spmv( A , x, y ) + call spmv( A , x, y , alpha, beta ) #:endif + y = merge( zero_${s}$, y, di_ ) end subroutine - subroutine default_preconditionner(x,y) + subroutine default_preconditionner(x,y,alpha,beta) ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) - y = x + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + y = merge( zero_${s}$, x, di_ ) end subroutine end subroutine From 5e15a333e50f99ec61db15da92a5795584cdec08 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Mon, 28 Apr 2025 23:20:53 +0200 Subject: [PATCH 14/26] add forward/backward solvers for preconditionning --- src/stdlib_linalg_iterative_aux.fypp | 137 ++++++++++++++++++ src/stdlib_linalg_iterative_solvers.fypp | 37 ++++- src/stdlib_linalg_iterative_solvers_cg.fypp | 2 +- src/stdlib_linalg_iterative_solvers_pccg.fypp | 7 +- 4 files changed, 180 insertions(+), 3 deletions(-) create mode 100644 src/stdlib_linalg_iterative_aux.fypp diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp new file mode 100644 index 000000000..057ce5c48 --- /dev/null +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -0,0 +1,137 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set MATRIX_TYPES = ["dense", "CSR"] +#:set RANKS = range(1, 2+1) + +submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_aux + use stdlib_kinds + use stdlib_sparse + use stdlib_constants + use stdlib_linalg_iterative_solvers + use stdlib_constants + implicit none + + integer, parameter :: ilp = int32 + +contains + + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_forward_triangular_dense_${s}$(L,b,x) + ${t}$, intent(in) :: L(:,:) + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$ :: aux + + integer(ilp) :: i, j, m + + do i = 1, size(L,dim=1) + aux = zero_${s}$ + do j = 1, i-1 + aux = aux + L(i,j)*x(j) + end do + x(i) = b(i) - aux + end do + end subroutine + + module subroutine solve_backward_triangular_dense_${s}$(U,b,x) + ${t}$, intent(in) :: U(:,:) + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$ :: aux + ${t}$ :: baux(size(x)) + + integer(ilp) :: i, j, m + + baux = zero_${s}$ + do i = size(U,dim=1), 1, -1 + x(i) = x(i) - baux(i) + do j = i+1, size(U,dim=2) + baux(j) = baux(j) + U(i,j)*x(i) + end do + end do + end subroutine + + #:endfor + + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_forward_triangular_csr_${s}$(L,b,x) + type(CSR_${s}$_type), intent(in) :: L + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$ :: aux + + integer(ilp) :: i, j, m + select case (L%storage) + case(sparse_full) + do i = 1, L%nrows + aux = zero_${s}$ + do m = L%rowptr(i), L%rowptr(i+1)-1 + j = L%col(m) + if(j>i) cycle !> skip upper part of the matrix + aux = aux + L%data(m)*x(j) + end do + x(i) = b(i) - aux + end do + case(sparse_lower) + do i = 1, L%nrows + aux = zero_${s}$ + do m = L%rowptr(i), L%rowptr(i+1)-2 + j = L%col(m) + aux = aux + L%data(m)*x(j) + end do + x(i) = b(i) - aux + end do + case(sparse_upper) !> treates as lower triangular (thus transpose) + do i = 1, L%nrows + x(i) = b(i) + do m = L%rowptr(i)+1, L%rowptr(i+1)-1 + j = L%col(m) + x(j) = x(j) - L%data(m)*x(i) + end do + end do + end select + end subroutine + + module subroutine solve_backward_triangular_csr_${s}$(U,b,x) + type(CSR_${s}$_type), intent(in) :: U + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + ${t}$ :: aux + ${t}$ :: baux(size(x)) + + integer(ilp) :: i, j, m + + baux = zero_${s}$ + select case (U%storage) + case(sparse_full) + do i = U%nrows, 1, -1 + x(i) = b(i) - baux(i) + do m = U%rowptr(i), U%rowptr(i+1)-1 + j = U%col(m) + if(j skip lower part of the matrix + baux(j) = baux(j) + U%data(m)*x(i) + end do + end do + case(sparse_lower) + do i = U%nrows, 1, -1 + x(i) = b(i) - baux(i) + do m = U%rowptr(i), U%rowptr(i+1)-2 + j = U%col(m) + baux(j) = baux(j) + U%data(m)*x(i) + end do + end do + case(sparse_upper) + do i = U%nrows, 1, -1 + x(i) = b(i) + do m = U%rowptr(i)+1, U%rowptr(i+1)-1 + j = U%col(m) + x(i) = x(i) - U%data(m)*x(j) + end do + end do + end select + end subroutine + + #:endfor + + +end submodule stdlib_linalg_iterative_aux \ No newline at end of file diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 701f2c5b4..eb8300ac4 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -1,7 +1,7 @@ #: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 MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"] +#:set MATRIX_TYPES = ["dense", "CSR"] #:set RANKS = range(1, 2+1) module stdlib_linalg_iterative_solvers @@ -131,6 +131,41 @@ module stdlib_linalg_iterative_solvers end interface public :: solve_pccg + interface solve_forward_triangular + !> Solve the system L x = b, where L is a strictly lower triangular matrix (diagonal assumed = 1). + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_forward_triangular_dense_${s}$(L,b,x) + ${t}$, intent(in) :: L(:,:) + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + end subroutine + module subroutine solve_forward_triangular_csr_${s}$(L,b,x) + type(CSR_${s}$_type), intent(in) :: L + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + end subroutine + #:endfor + end interface + public :: solve_forward_triangular + + + interface solve_backward_triangular + !> Solve the system U x = b, where U is a strictly upper triangular matrix (diagonal assumed = 1). + #:for k, t, s in R_KINDS_TYPES + module subroutine solve_backward_triangular_dense_${s}$(U,b,x) + ${t}$, intent(in) :: U(:,:) + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + end subroutine + module subroutine solve_backward_triangular_csr_${s}$(U,b,x) + type(CSR_${s}$_type), intent(in) :: U + ${t}$, intent(in) :: b(:) + ${t}$, intent(inout) :: x(:) + end subroutine + #:endfor + end interface + public :: solve_backward_triangular + contains diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index 241403662..eca50bf6c 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -1,7 +1,7 @@ #: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 MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"] +#:set MATRIX_TYPES = ["dense", "CSR"] #:set RANKS = range(1, 2+1) submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index ec1b83d6a..69628615e 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -1,7 +1,7 @@ #: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 MATRIX_TYPES = ["dense", "CSR", "CSC", "ELL", "SELLC"] +#:set MATRIX_TYPES = ["dense", "CSR"] #:set RANKS = range(1, 2+1) submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg @@ -115,6 +115,11 @@ contains allocate( M_ ) M_%apply => default_preconditionner end if + !> TODO: factorize the preconditionner + !> Jacobi: M = diag(A)^{-1} + !> SOR: M = (w/(2-w)) (D/w+L) D^{-1} (D/w+U) ; A = D + L + U + !> ILU: M = ( LU )^{-1} ; A = LU + !> ICHOL: M = (L*D*U)^{-1} ; A = L*D*U if(present(di))then di_ => di else From 8c2aa908823a9e44a6dcc06cccb67cd2d9abed72 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Tue, 29 Apr 2025 18:22:12 +0200 Subject: [PATCH 15/26] fix solve forward/backward --- src/stdlib_linalg_iterative_aux.fypp | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp index 057ce5c48..502c15094 100644 --- a/src/stdlib_linalg_iterative_aux.fypp +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -20,16 +20,15 @@ contains ${t}$, intent(in) :: L(:,:) ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) - ${t}$ :: aux integer(ilp) :: i, j, m - do i = 1, size(L,dim=1) - aux = zero_${s}$ - do j = 1, i-1 - aux = aux + L(i,j)*x(j) + x = zero_${s}$ + do j = 1, size(L,dim=1) + x(j) = x(j) + b(j) + do i = j+1, size(L,dim=1) + x(i) = x(i) - L(i,j)*x(j) end do - x(i) = b(i) - aux end do end subroutine @@ -37,16 +36,13 @@ contains ${t}$, intent(in) :: U(:,:) ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) - ${t}$ :: aux - ${t}$ :: baux(size(x)) integer(ilp) :: i, j, m - baux = zero_${s}$ - do i = size(U,dim=1), 1, -1 - x(i) = x(i) - baux(i) - do j = i+1, size(U,dim=2) - baux(j) = baux(j) + U(i,j)*x(i) + do j = size(U,dim=1), 1, -1 + x(j) = b(j) + do i = 1, j-1 + x(i) = x(i) - U(i,j)*x(j) end do end do end subroutine From 517291ae8dc30907320d192b99661af8caf5eea1 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Wed, 30 Apr 2025 08:30:13 +0200 Subject: [PATCH 16/26] fix colum index --- src/stdlib_linalg_iterative_aux.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp index 502c15094..82c83336b 100644 --- a/src/stdlib_linalg_iterative_aux.fypp +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -24,7 +24,7 @@ contains integer(ilp) :: i, j, m x = zero_${s}$ - do j = 1, size(L,dim=1) + do j = 1, size(L,dim=2) x(j) = x(j) + b(j) do i = j+1, size(L,dim=1) x(i) = x(i) - L(i,j)*x(j) @@ -39,7 +39,7 @@ contains integer(ilp) :: i, j, m - do j = size(U,dim=1), 1, -1 + do j = size(U,dim=2), 1, -1 x(j) = b(j) do i = 1, j-1 x(i) = x(i) - U(i,j)*x(j) From 367987a6f876751eaacf6c3855d42860d76a03c9 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Wed, 30 Apr 2025 21:49:42 +0200 Subject: [PATCH 17/26] add preconditionners --- src/stdlib_linalg_iterative_solvers.fypp | 3 +- src/stdlib_linalg_iterative_solvers_pccg.fypp | 76 ++++++++++++++++--- 2 files changed, 67 insertions(+), 12 deletions(-) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index eb8300ac4..b1451b8b0 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -111,7 +111,7 @@ module stdlib_linalg_iterative_solvers interface solve_pccg #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,M,workspace) + module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -123,6 +123,7 @@ module stdlib_linalg_iterative_solvers logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter logical, intent(in), optional :: restart + integer, intent(in), optional :: precond !> preconditionner method flag class(linop_${s}$), optional , intent(in), target :: M !> preconditionner type(solver_workspace_${s}$), optional, intent(inout), target :: workspace end subroutine diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 69628615e..126685e4b 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -7,10 +7,18 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg use stdlib_kinds use stdlib_sparse + use stdlib_linalg, only: diag use stdlib_constants use stdlib_linalg_iterative_solvers implicit none + enum, bind(c) + enumerator :: none = 0 + enumerator :: jacobi + enumerator :: ssor + enumerator :: ichol + end enum + contains #:for k, t, s in R_KINDS_TYPES @@ -78,7 +86,7 @@ contains #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,M,workspace) + module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -90,6 +98,7 @@ contains logical(1), intent(in), optional, target :: di(:) integer, intent(in), optional :: maxiter logical, intent(in), optional :: restart + integer, intent(in), optional :: precond !> preconditionner method flag class(linop_${s}$), optional , intent(in), target :: M !> preconditionner type(solver_workspace_${s}$), optional, intent(inout), target :: workspace !------------------------- @@ -101,25 +110,53 @@ contains logical :: restart_ logical(1), pointer :: di_(:) !------------------------- + ! working data for preconditionner + integer :: precond_ + ${t}$, allocatable :: diagonal(:) + #:if matrix == "dense" + ${t}$, allocatable:: L(:,:) !> lower triangular + #:else + type(${matrix}$_${s}$_type) :: L !> lower triangular + #:endif + !------------------------- n = size(b) maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter restart_ = .true.; if(present(restart)) restart_ = restart tol_ = 1.e-4_${s}$; if(present(tol)) tol_ = tol - + precond_ = none ; if(present(precond)) precond_ = precond !------------------------- ! internal memory setup - op%apply => default_matvec + op%apply => matvec if(present(M)) then M_ => M else allocate( M_ ) - M_%apply => default_preconditionner + select case(precond_) + case(jacobi) + #:if matrix == "dense" + diagonal = diag(A) + #:else + call diag(A,diagonal) + #:endif + where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal + M_%apply => precond_jacobi + case(ssor) + #:if matrix == "dense" + diagonal = diag(A) + #:else + call diag(A,diagonal) + #:endif + where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal + ! TODO: Extract the lower triangular part of A in L + M_%apply => precond_ldl + case(ichol) + ! TODO: Add LDL^t factorization + M_%apply => precond_ldl + case default + M_%apply => precond_none + end select end if - !> TODO: factorize the preconditionner - !> Jacobi: M = diag(A)^{-1} - !> SOR: M = (w/(2-w)) (D/w+L) D^{-1} (D/w+U) ; A = D + L + U - !> ILU: M = ( LU )^{-1} ; A = LU - !> ICHOL: M = (L*D*U)^{-1} ; A = L*D*U + if(present(di))then di_ => di else @@ -151,7 +188,7 @@ contains workspace_ => null() contains - subroutine default_matvec(x,y,alpha,beta) + subroutine matvec(x,y,alpha,beta) ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) ${t}$, intent(in) :: alpha @@ -163,13 +200,30 @@ contains #:endif y = merge( zero_${s}$, y, di_ ) end subroutine - subroutine default_preconditionner(x,y,alpha,beta) + + subroutine precond_none(x,y,alpha,beta) ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) ${t}$, intent(in) :: alpha ${t}$, intent(in) :: beta y = merge( zero_${s}$, x, di_ ) end subroutine + subroutine precond_jacobi(x,y,alpha,beta) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + y = merge( zero_${s}$, diagonal * x, di_ ) !> inverted diagonal + end subroutine + subroutine precond_ldl(x,y,alpha,beta) + ${t}$, intent(in) :: x(:) + ${t}$, intent(inout) :: y(:) + ${t}$, intent(in) :: alpha + ${t}$, intent(in) :: beta + call solve_forward_triangular( L , x , y ) + y = merge( zero_${s}$, diagonal * y, di_ ) !> inverted diagonal + call solve_backward_triangular( L , y , y ) + end subroutine end subroutine #:endfor From 05c076b4224829a312d7e0406fa9b89d82a66388 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Thu, 1 May 2025 14:23:28 +0200 Subject: [PATCH 18/26] fix cmake --- src/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 55224c833..a9b9f79c2 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -92,6 +92,7 @@ set(cppFiles stdlib_linalg_blas.fypp stdlib_linalg_lapack.fypp + stdlib_linalg_iterative_aux.fypp stdlib_linalg_iterative_solvers.fypp stdlib_linalg_iterative_solvers_cg.fypp stdlib_linalg_iterative_solvers_pccg.fypp From f027017e265c54fc94a39d8536d9d7707b81950c Mon Sep 17 00:00:00 2001 From: jalvesz Date: Thu, 1 May 2025 15:55:10 +0200 Subject: [PATCH 19/26] add factorizations --- src/stdlib_linalg_iterative_aux.fypp | 179 ++++++++++++++++++ src/stdlib_linalg_iterative_solvers.fypp | 34 ++++ src/stdlib_linalg_iterative_solvers_pccg.fypp | 23 ++- 3 files changed, 228 insertions(+), 8 deletions(-) diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp index 82c83336b..4c70b4614 100644 --- a/src/stdlib_linalg_iterative_aux.fypp +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -130,4 +130,183 @@ contains #:endfor + !============================================================== + ! Factorization for preconditioners + !============================================================== + #:for k, t, s in R_KINDS_TYPES + module subroutine factorize_ssor_dense_${s}$(A, w, L, D) + ${t}$, intent(in) :: A(:,:) + ${t}$, intent(in) :: w + ${t}$, intent(out) :: L(:,:) + ${t}$, intent(out) :: D(:) + + integer(ilp) :: i, j, n + ${t}$ :: inv_w + + n = size(A,dim=1) + inv_w = 1._${s}$ / w + do j = 1, n + D(j) = A(j,j) * (2._${s}$-w)*inv_w + do i = 1, j-1 + L(i,j) = 0._${s}$ + end do + L(j,j) = L(j,j) * inv_w + do i = j+1, n + L(i,j) = A(i,j) + end do + end do + + end subroutine + + #:endfor + + #:for k, t, s in R_KINDS_TYPES + module subroutine factorize_ssor_csr_${s}$(A, w, L, D) + type(CSR_${s}$_type), intent(in) :: A + ${t}$, intent(in) :: w + type(CSR_${s}$_type), intent(out) :: L + ${t}$, intent(out) :: D(:) + + integer(ilp) :: i, j, n + ${t}$ :: inv_w + + inv_w = 1._${s}$ / w + L%data = A%data + select case(A%storage) + case(sparse_lower) + do i = 1, A%nrows + D(i) = A%data( A%rowptr(i+1)-1 ) * (2._${s}$-w)*inv_w + L%data( A%rowptr(i+1)-1 ) = L%data( A%rowptr(i+1)-1 ) * inv_w + end do + case(sparse_upper) + do i = 1, A%nrows + D(i) = A%data( A%rowptr(i) ) * (2._${s}$-w)*inv_w + L%data( A%rowptr(i) ) = L%data( A%rowptr(i) ) * inv_w + end do + case(sparse_full) + do i = 1, A%nrows + do j = A%rowptr(i), A%rowptr(i+1)-1 + if( A%col(j) == i ) then + D(i) = A%data(j) * (2._${s}$-w)*inv_w + L%data(j) = L%data(j) * inv_w + exit + end if + end do + end do + end select + + end subroutine + + #:endfor + + !> Bunch-Kaufman factorization of a symmetric positive definite matrix A. + !> The matrix A is assumed to be symmetric and positive definite. + #:for k, t, s in R_KINDS_TYPES + module subroutine factorize_ldlt_dense_${s}$(A, L, D) + ${t}$, intent(in) :: A(:,:) + ${t}$, intent(out) :: L(:,:) + ${t}$, intent(out) :: D(:) + + ${t}$:: auxsum + integer(ilp) :: i, j, k, n + + ! Initialize L to zero + n = size(A,dim=1) + L = zero_${s}$ + + do i = 1, n + ! Compute D(i) + auxsum = A(i, i) + do k = 1, i - 1 + auxsum = auxsum - (L(i, k) ** 2) * D(k) + end do + D(i) = auxsum + + ! Check for positive definiteness + if (D(i) <= zero_${s}$ ) D(i) = one_${s}$ + + ! Set unit diagonal + L(i, i) = one_${s}$ + + ! Compute L(j, i) for j = i+1 to n + do j = i + 1, n + auxsum = A(j, i) + do k = 1, i - 1 + auxsum = auxsum - L(j, k) * L(i, k) * D(k) + end do + L(j, i) = auxsum / D(i) + end do + end do + end subroutine + + #:endfor + + #:for k, t, s in R_KINDS_TYPES + module subroutine factorize_ldlt_csr_${s}$(A, L, D) + type(CSR_${s}$_type), intent(in) :: A + type(CSR_${s}$_type), intent(out) :: L + ${t}$, intent(out) :: D(:) + + ${t}$:: auxsum + integer(ilp) :: i, j, k, n, p, q + + integer, allocatable :: col_to_pos(:) + ${t}$, allocatable :: temp(:) + + allocate(col_to_pos(n), source=-1) + allocate(temp(n),source=zero_${s}$) + + n = A%nrows + do i = 1, n + ! Zero the temp workspace + temp = zero_${s}$ + + ! Load row i into temp (lower triangle only) + do p = A%rowptr(i), A%rowptr(i+1) - 1 + j = A%col(p) + if (j > i) cycle + temp(j) = A%data(p) + col_to_pos(j) = p ! map column to position in L%data + end do + + ! Compute L(i, k) * D(k) * L(i, k) contributions + do p = A%rowptr(i), A%rowptr(i+1) - 1 + k = A%col(p) + if (k >= i) exit + + auxsum = temp(k) + + ! Compute L(i, 1:k-1) ⋅ D(1:k-1) ⋅ L(k, 1:k-1) + do q = A%rowptr(k), A%rowptr(k+1) - 1 + j = A%col(q) + if (j >= k) exit + auxsum = auxsum - L%data(col_to_pos(j)) * D(j) * L%data(q) + end do + + L%data(p) = auxsum / D(k) + temp(k) = zero_${s}$ ! Clear temp slot + end do + + ! Compute D(i) + auxsum = temp(i) + do p = A%rowptr(i), A%rowptr(i+1) - 1 + k = A%col(p) + if (k >= i) exit + auxsum = auxsum - L%data(p)**2 * D(k) + end do + + D(i) = merge( auxsum, one_${s}$ , auxsum > zero_${s}$) + + ! Set L(i, i) = 1.0 + L%data(col_to_pos(i)) = one_${s}$ + + ! Reset col_to_pos map + do p = A%rowptr(i), A%rowptr(i+1) - 1 + col_to_pos(A%col(p)) = -1 + end do + end do + end subroutine + + #:endfor + end submodule stdlib_linalg_iterative_aux \ No newline at end of file diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index b1451b8b0..9c5d20ad1 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -167,6 +167,40 @@ module stdlib_linalg_iterative_solvers end interface public :: solve_backward_triangular + interface factorize_ldlt + #:for k, t, s in R_KINDS_TYPES + module subroutine factorize_ldlt_dense_${s}$(A, L, D) + ${t}$, intent(in) :: A(:,:) + ${t}$, intent(out) :: L(:,:) + ${t}$, intent(out) :: D(:) + end subroutine + module subroutine factorize_ldlt_csr_${s}$(A, L, D) + type(CSR_${s}$_type), intent(in) :: A + type(CSR_${s}$_type), intent(out) :: L + ${t}$, intent(out) :: D(:) + end subroutine + #:endfor + end interface + public :: factorize_ldlt + + interface factorize_ssor + #:for k, t, s in R_KINDS_TYPES + module subroutine factorize_ssor_dense_${s}$(A, w, L, D) + ${t}$, intent(in) :: A(:,:) + ${t}$, intent(in) :: w + ${t}$, intent(out) :: L(:,:) + ${t}$, intent(out) :: D(:) + end subroutine + module subroutine factorize_ssor_csr_${s}$(A, w, L, D) + type(CSR_${s}$_type), intent(in) :: A + ${t}$, intent(in) :: w + type(CSR_${s}$_type), intent(out) :: L + ${t}$, intent(out) :: D(:) + end subroutine + #:endfor + end interface + public :: factorize_ssor + contains diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 126685e4b..4b4a68b75 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -7,7 +7,6 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg use stdlib_kinds use stdlib_sparse - use stdlib_linalg, only: diag use stdlib_constants use stdlib_linalg_iterative_solvers implicit none @@ -16,7 +15,7 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg enumerator :: none = 0 enumerator :: jacobi enumerator :: ssor - enumerator :: ichol + enumerator :: ldlt end enum contains @@ -88,6 +87,7 @@ contains #:for k, t, s in R_KINDS_TYPES module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) #:if matrix == "dense" + use stdlib_linalg, only: diag ${t}$, intent(in) :: A(:,:) #:else type(${matrix}$_${s}$_type), intent(in) :: A @@ -126,11 +126,13 @@ contains precond_ = none ; if(present(precond)) precond_ = precond !------------------------- ! internal memory setup - op%apply => matvec + ! Preconditionner if(present(M)) then M_ => M else allocate( M_ ) + allocate(diagonal(n),source=zero_${s}$) + select case(precond_) case(jacobi) #:if matrix == "dense" @@ -138,7 +140,6 @@ contains #:else call diag(A,diagonal) #:endif - where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal M_%apply => precond_jacobi case(ssor) #:if matrix == "dense" @@ -146,23 +147,29 @@ contains #:else call diag(A,diagonal) #:endif - where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal - ! TODO: Extract the lower triangular part of A in L + L = A !> copy A structure to L + call factorize_ssor( A , one_${s}$ , L, diagonal ) M_%apply => precond_ldl - case(ichol) - ! TODO: Add LDL^t factorization + case(ldlt) + L = A !> copy A structure to L + call factorize_ldlt( A , L, diagonal ) M_%apply => precond_ldl case default M_%apply => precond_none end select + where(abs(diagonal)>epsilon(zero_${s}$)) diagonal = one_${s}$/diagonal end if + ! matvec for the operator + op%apply => matvec + ! direchlet boundary conditions mask if(present(di))then di_ => di else allocate(di_(n),source=.false._1) end if + ! workspace for the solver if(present(workspace)) then workspace_ => workspace else From 7fd2586c60714cc99c76268500cc42e6125946c7 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Thu, 1 May 2025 21:59:14 +0200 Subject: [PATCH 20/26] backward solve to use Lt marix --- src/stdlib_linalg_iterative_aux.fypp | 38 ++++++++++++------------ src/stdlib_linalg_iterative_solvers.fypp | 8 ++--- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp index 4c70b4614..7361bb5a9 100644 --- a/src/stdlib_linalg_iterative_aux.fypp +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -32,17 +32,17 @@ contains end do end subroutine - module subroutine solve_backward_triangular_dense_${s}$(U,b,x) - ${t}$, intent(in) :: U(:,:) + module subroutine solve_backward_triangular_dense_${s}$(Lt,b,x) + ${t}$, intent(in) :: Lt(:,:) ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) integer(ilp) :: i, j, m - do j = size(U,dim=2), 1, -1 + do j = size(Lt,dim=2), 1, -1 x(j) = b(j) do i = 1, j-1 - x(i) = x(i) - U(i,j)*x(j) + x(i) = x(i) - Lt(j,i)*x(j) end do end do end subroutine @@ -88,8 +88,8 @@ contains end select end subroutine - module subroutine solve_backward_triangular_csr_${s}$(U,b,x) - type(CSR_${s}$_type), intent(in) :: U + module subroutine solve_backward_triangular_csr_${s}$(Lt,b,x) + type(CSR_${s}$_type), intent(in) :: Lt ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) ${t}$ :: aux @@ -98,30 +98,30 @@ contains integer(ilp) :: i, j, m baux = zero_${s}$ - select case (U%storage) + select case (Lt%storage) case(sparse_full) - do i = U%nrows, 1, -1 + do i = Lt%nrows, 1, -1 x(i) = b(i) - baux(i) - do m = U%rowptr(i), U%rowptr(i+1)-1 - j = U%col(m) + do m = Lt%rowptr(i), Lt%rowptr(i+1)-1 + j = Lt%col(m) if(j skip lower part of the matrix - baux(j) = baux(j) + U%data(m)*x(i) + baux(j) = baux(j) + Lt%data(m)*x(i) end do end do case(sparse_lower) - do i = U%nrows, 1, -1 + do i = Lt%nrows, 1, -1 x(i) = b(i) - baux(i) - do m = U%rowptr(i), U%rowptr(i+1)-2 - j = U%col(m) - baux(j) = baux(j) + U%data(m)*x(i) + do m = Lt%rowptr(i), Lt%rowptr(i+1)-2 + j = Lt%col(m) + baux(j) = baux(j) + Lt%data(m)*x(i) end do end do case(sparse_upper) - do i = U%nrows, 1, -1 + do i = Lt%nrows, 1, -1 x(i) = b(i) - do m = U%rowptr(i)+1, U%rowptr(i+1)-1 - j = U%col(m) - x(i) = x(i) - U%data(m)*x(j) + do m = Lt%rowptr(i)+1, Lt%rowptr(i+1)-1 + j = Lt%col(m) + x(i) = x(i) - Lt%data(m)*x(j) end do end do end select diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 9c5d20ad1..2fe9e85f0 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -153,13 +153,13 @@ module stdlib_linalg_iterative_solvers interface solve_backward_triangular !> Solve the system U x = b, where U is a strictly upper triangular matrix (diagonal assumed = 1). #:for k, t, s in R_KINDS_TYPES - module subroutine solve_backward_triangular_dense_${s}$(U,b,x) - ${t}$, intent(in) :: U(:,:) + module subroutine solve_backward_triangular_dense_${s}$(Lt,b,x) + ${t}$, intent(in) :: Lt(:,:) ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) end subroutine - module subroutine solve_backward_triangular_csr_${s}$(U,b,x) - type(CSR_${s}$_type), intent(in) :: U + module subroutine solve_backward_triangular_csr_${s}$(Lt,b,x) + type(CSR_${s}$_type), intent(in) :: Lt ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) end subroutine From c596ac085352a5330df43f7eaf994821eeb4894b Mon Sep 17 00:00:00 2001 From: jalvesz Date: Thu, 1 May 2025 21:59:48 +0200 Subject: [PATCH 21/26] start unit testing --- test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_solve_iterative.fypp | 129 +++++++++++++++++++ 2 files changed, 131 insertions(+) create mode 100644 test/linalg/test_linalg_solve_iterative.fypp diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index c496bd2b3..c1ab37f85 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -13,6 +13,7 @@ set( "test_linalg_determinant.fypp" "test_linalg_qr.fypp" "test_linalg_schur.fypp" + "test_linalg_solve_iterative.fypp" "test_linalg_svd.fypp" "test_linalg_matrix_property_checks.fypp" "test_linalg_sparse.fypp" @@ -32,6 +33,7 @@ ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) ADDTEST(linalg_qr) ADDTEST(linalg_schur) +ADDTEST(linalg_solve_iterative) ADDTEST(linalg_svd) ADDTEST(blas_lapack) ADDTEST(linalg_sparse) diff --git a/test/linalg/test_linalg_solve_iterative.fypp b/test/linalg/test_linalg_solve_iterative.fypp new file mode 100644 index 000000000..2859e158c --- /dev/null +++ b/test/linalg/test_linalg_solve_iterative.fypp @@ -0,0 +1,129 @@ +#:include "common.fypp" +#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX)) +#:set MATRIX_TYPES = ["dense", "CSR"] +! Test linear system iterative solvers +module test_linalg_solve_iterative + use stdlib_kinds + use stdlib_sparse + use stdlib_linalg_iterative_solvers + use testdrive, only: error_type, check, new_unittest, unittest_type + + implicit none + private + + public :: test_linear_systems + + contains + + subroutine test_linear_systems(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + tests = [ new_unittest("factorize_ldlt",test_factorize_ldlt), & + new_unittest("solve_ldlt",test_solve_ldlt) ] + + end subroutine test_linear_systems + + + !> Simple linear system + subroutine test_factorize_ldlt(error) + type(error_type), allocatable, intent(out) :: error + + #:for k, t, s in R_KINDS_TYPES + block + ${t}$, parameter :: tol = 1.e-4_${k}$ + ${t}$ :: A(5,5) = reshape([${t}$ :: 5, -2, 0, -2, -2, & + -2, 5, -2, 0, 0, & + 0, -2, 5, -2, 0, & + -2, 0, -2, 5, -2, & + -2, 0, 0, -2, 5], [5,5]) + ${t}$ :: Lref(5,5) = transpose(reshape([${t}$ :: 1. , 0. , 0. , 0. , 0.,& + -0.4, 1. , 0. , 0. , 0.,& + 0. , -0.47619048, 1. , 0. , 0.,& + -0.4, -0.19047619, -0.58823529, 1. , 0., & + -0.4, -0.19047619, -0.09411765, -1.2, 1.], [5,5])) + + ${t}$ :: Dref(5) = [${t}$ :: 5, 4.2, 4.04761904761904744987, 2.64705882352941124225, 0.2] + + ${t}$ :: L(5,5), D(5) + + call factorize_ldlt(A, L, D) + + call check(error, all(abs(L-Lref) 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_solve_iterative From acefaafc77eb6607629d7244834d46f0e4828f15 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Mon, 5 May 2025 22:44:45 +0200 Subject: [PATCH 22/26] review csr factorization --- src/stdlib_linalg_iterative_aux.fypp | 102 +++++++++-------------- src/stdlib_linalg_iterative_solvers.fypp | 16 ++-- 2 files changed, 47 insertions(+), 71 deletions(-) diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp index 7361bb5a9..ac49bee96 100644 --- a/src/stdlib_linalg_iterative_aux.fypp +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -137,8 +137,8 @@ contains module subroutine factorize_ssor_dense_${s}$(A, w, L, D) ${t}$, intent(in) :: A(:,:) ${t}$, intent(in) :: w - ${t}$, intent(out) :: L(:,:) - ${t}$, intent(out) :: D(:) + ${t}$, intent(inout) :: L(:,:) + ${t}$, intent(inout) :: D(:) integer(ilp) :: i, j, n ${t}$ :: inv_w @@ -164,8 +164,8 @@ contains module subroutine factorize_ssor_csr_${s}$(A, w, L, D) type(CSR_${s}$_type), intent(in) :: A ${t}$, intent(in) :: w - type(CSR_${s}$_type), intent(out) :: L - ${t}$, intent(out) :: D(:) + type(CSR_${s}$_type), intent(inout) :: L + ${t}$, intent(inout) :: D(:) integer(ilp) :: i, j, n ${t}$ :: inv_w @@ -204,8 +204,8 @@ contains #:for k, t, s in R_KINDS_TYPES module subroutine factorize_ldlt_dense_${s}$(A, L, D) ${t}$, intent(in) :: A(:,:) - ${t}$, intent(out) :: L(:,:) - ${t}$, intent(out) :: D(:) + ${t}$, intent(inout) :: L(:,:) + ${t}$, intent(inout) :: D(:) ${t}$:: auxsum integer(ilp) :: i, j, k, n @@ -244,67 +244,43 @@ contains #:for k, t, s in R_KINDS_TYPES module subroutine factorize_ldlt_csr_${s}$(A, L, D) type(CSR_${s}$_type), intent(in) :: A - type(CSR_${s}$_type), intent(out) :: L - ${t}$, intent(out) :: D(:) + type(CSR_${s}$_type), intent(inout) :: L + ${t}$, intent(inout) :: D(:) - ${t}$:: auxsum - integer(ilp) :: i, j, k, n, p, q - - integer, allocatable :: col_to_pos(:) - ${t}$, allocatable :: temp(:) - - allocate(col_to_pos(n), source=-1) - allocate(temp(n),source=zero_${s}$) + ${t}$:: aux_diag, nondiag + integer(ilp) :: i, j, k, m, n, p, q, r, ad1, ad2 - n = A%nrows - do i = 1, n - ! Zero the temp workspace - temp = zero_${s}$ - - ! Load row i into temp (lower triangle only) - do p = A%rowptr(i), A%rowptr(i+1) - 1 - j = A%col(p) - if (j > i) cycle - temp(j) = A%data(p) - col_to_pos(j) = p ! map column to position in L%data - end do + D = zero_${s}$ + L%data = zero_${s}$ - ! Compute L(i, k) * D(k) * L(i, k) contributions - do p = A%rowptr(i), A%rowptr(i+1) - 1 - k = A%col(p) - if (k >= i) exit - - auxsum = temp(k) - - ! Compute L(i, 1:k-1) ⋅ D(1:k-1) ⋅ L(k, 1:k-1) - do q = A%rowptr(k), A%rowptr(k+1) - 1 - j = A%col(q) - if (j >= k) exit - auxsum = auxsum - L%data(col_to_pos(j)) * D(j) * L%data(q) - end do - - L%data(p) = auxsum / D(k) - temp(k) = zero_${s}$ ! Clear temp slot - end do - - ! Compute D(i) - auxsum = temp(i) - do p = A%rowptr(i), A%rowptr(i+1) - 1 - k = A%col(p) - if (k >= i) exit - auxsum = auxsum - L%data(p)**2 * D(k) - end do - - D(i) = merge( auxsum, one_${s}$ , auxsum > zero_${s}$) - - ! Set L(i, i) = 1.0 - L%data(col_to_pos(i)) = one_${s}$ - - ! Reset col_to_pos map - do p = A%rowptr(i), A%rowptr(i+1) - 1 - col_to_pos(A%col(p)) = -1 + select case(A%storage) + case(sparse_lower) + do i = 1, A%nrows + ad1 = A%rowptr(i) + ad2 = A%rowptr(i+1)-1 + L%data(ad2) = one_${s}$ + aux_diag = zero_${s}$ + do m = ad1, ad2-1 + j = A%col(m) + nondiag = A%data(m) + q = ad1 + do p = A%rowptr(j), A%rowptr(j+1)-2 + k = A%col(p) + do r = q, m-1 + n = A%col(r) + if( k == n ) then + nondiag = nondiag - L%data(r)*L%data(p)*D(k) ! Aij - Lik Ljk Dk + q = r + 1 + end if + end do + end do + L%data(m) = nondiag / D(j) ! Lij = (Aij - Lik Ljk Dk)/Dj + aux_diag = aux_diag + (L%data(m) * nondiag) + end do + D(i) = A%data(ad2) - aux_diag ! Dj = Ajj - Ljk^2 Dk + if (D(i) <= zero_${s}$) D(i)= one_${s}$ end do - end do + end select end subroutine #:endfor diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 2fe9e85f0..929e1e051 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -171,13 +171,13 @@ module stdlib_linalg_iterative_solvers #:for k, t, s in R_KINDS_TYPES module subroutine factorize_ldlt_dense_${s}$(A, L, D) ${t}$, intent(in) :: A(:,:) - ${t}$, intent(out) :: L(:,:) - ${t}$, intent(out) :: D(:) + ${t}$, intent(inout) :: L(:,:) + ${t}$, intent(inout) :: D(:) end subroutine module subroutine factorize_ldlt_csr_${s}$(A, L, D) type(CSR_${s}$_type), intent(in) :: A - type(CSR_${s}$_type), intent(out) :: L - ${t}$, intent(out) :: D(:) + type(CSR_${s}$_type), intent(inout) :: L + ${t}$, intent(inout) :: D(:) end subroutine #:endfor end interface @@ -188,14 +188,14 @@ module stdlib_linalg_iterative_solvers module subroutine factorize_ssor_dense_${s}$(A, w, L, D) ${t}$, intent(in) :: A(:,:) ${t}$, intent(in) :: w - ${t}$, intent(out) :: L(:,:) - ${t}$, intent(out) :: D(:) + ${t}$, intent(inout) :: L(:,:) + ${t}$, intent(inout) :: D(:) end subroutine module subroutine factorize_ssor_csr_${s}$(A, w, L, D) type(CSR_${s}$_type), intent(in) :: A ${t}$, intent(in) :: w - type(CSR_${s}$_type), intent(out) :: L - ${t}$, intent(out) :: D(:) + type(CSR_${s}$_type), intent(inout) :: L + ${t}$, intent(inout) :: D(:) end subroutine #:endfor end interface From f413cbf64b3e01c0df7cb5b32d5de54a7256f902 Mon Sep 17 00:00:00 2001 From: Jose Alves Date: Mon, 5 May 2025 22:47:06 +0200 Subject: [PATCH 23/26] change name generic for kernel --- example/linalg/example_solve_custom.f90 | 4 ++-- src/stdlib_linalg_iterative_solvers.fypp | 12 ++++++------ src/stdlib_linalg_iterative_solvers_cg.fypp | 4 ++-- src/stdlib_linalg_iterative_solvers_pccg.fypp | 4 ++-- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index 49faf750c..e4be3c99e 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -3,7 +3,7 @@ module custom_solver use stdlib_sparse use stdlib_linalg_iterative_solvers, only: linop_dp, & solver_workspace_dp, & - solve_pccg_generic, & + solve_pccg_kernel, & size_wksp_pccg implicit none contains @@ -56,7 +56,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal !------------------------- ! main call to the solver - call solve_pccg_generic(op,M,b,x,tol_,maxiter_,workspace_) + call solve_pccg_kernel(op,M,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index 929e1e051..c9d6041db 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -57,9 +57,9 @@ module stdlib_linalg_iterative_solvers #:endfor end interface - interface solve_cg_generic + interface solve_cg_kernel #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) + module subroutine solve_cg_kernel_${s}$(A,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A ${t}$, intent(in) :: b(:) ${t}$, intent(inout) :: x(:) @@ -69,7 +69,7 @@ module stdlib_linalg_iterative_solvers end subroutine #:endfor end interface - public :: solve_cg_generic + public :: solve_cg_kernel interface solve_cg #:for matrix in MATRIX_TYPES @@ -93,9 +93,9 @@ module stdlib_linalg_iterative_solvers end interface public :: solve_cg - interface solve_pccg_generic + interface solve_pccg_kernel #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_generic_${s}$(A,M,b,x,tol,maxiter,workspace) + module subroutine solve_pccg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:) @@ -106,7 +106,7 @@ module stdlib_linalg_iterative_solvers end subroutine #:endfor end interface - public :: solve_pccg_generic + public :: solve_pccg_kernel interface solve_pccg #:for matrix in MATRIX_TYPES diff --git a/src/stdlib_linalg_iterative_solvers_cg.fypp b/src/stdlib_linalg_iterative_solvers_cg.fypp index eca50bf6c..77d45f6fc 100644 --- a/src/stdlib_linalg_iterative_solvers_cg.fypp +++ b/src/stdlib_linalg_iterative_solvers_cg.fypp @@ -14,7 +14,7 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_cg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_cg_generic_${s}$(A,b,x,tol,maxiter,workspace) + module subroutine solve_cg_kernel_${s}$(A,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A ${t}$, intent(in) :: b(:), tol ${t}$, intent(inout) :: x(:) @@ -112,7 +112,7 @@ contains ! main call to the solver if(restart_) x = zero_${s}$ x = merge( b, x, di_ ) !> copy dirichlet load conditions encoded in B and indicated by di - call solve_cg_generic(op,b,x,tol_,maxiter_,workspace_) + call solve_cg_kernel(op,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 4b4a68b75..ebc76d7c7 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -21,7 +21,7 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_generic_${s}$(A,M,b,x,tol,maxiter,workspace) + module subroutine solve_pccg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:), tol @@ -180,7 +180,7 @@ contains ! main call to the solver if(restart_) x = zero_${s}$ x = merge( b, x, di_ ) !> copy dirichlet load conditions encoded in B and indicated by di - call solve_pccg_generic(op,M_,b,x,tol_,maxiter_,workspace_) + call solve_pccg_kernel(op,M_,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup From ea7d35e5cfc9dcb62fe69eb9bb19ecffb75c0975 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sat, 10 May 2025 19:16:29 +0200 Subject: [PATCH 24/26] shorten factorization procedures names --- src/stdlib_linalg_iterative_aux.fypp | 8 ++++---- src/stdlib_linalg_iterative_solvers.fypp | 16 +++++++-------- src/stdlib_linalg_iterative_solvers_pccg.fypp | 20 +++++++++---------- test/linalg/test_linalg_solve_iterative.fypp | 5 ++--- 4 files changed, 24 insertions(+), 25 deletions(-) diff --git a/src/stdlib_linalg_iterative_aux.fypp b/src/stdlib_linalg_iterative_aux.fypp index ac49bee96..03da4e4b3 100644 --- a/src/stdlib_linalg_iterative_aux.fypp +++ b/src/stdlib_linalg_iterative_aux.fypp @@ -134,7 +134,7 @@ contains ! Factorization for preconditioners !============================================================== #:for k, t, s in R_KINDS_TYPES - module subroutine factorize_ssor_dense_${s}$(A, w, L, D) + module subroutine ssor_dense_${s}$(A, w, L, D) ${t}$, intent(in) :: A(:,:) ${t}$, intent(in) :: w ${t}$, intent(inout) :: L(:,:) @@ -161,7 +161,7 @@ contains #:endfor #:for k, t, s in R_KINDS_TYPES - module subroutine factorize_ssor_csr_${s}$(A, w, L, D) + module subroutine ssor_csr_${s}$(A, w, L, D) type(CSR_${s}$_type), intent(in) :: A ${t}$, intent(in) :: w type(CSR_${s}$_type), intent(inout) :: L @@ -202,7 +202,7 @@ contains !> Bunch-Kaufman factorization of a symmetric positive definite matrix A. !> The matrix A is assumed to be symmetric and positive definite. #:for k, t, s in R_KINDS_TYPES - module subroutine factorize_ldlt_dense_${s}$(A, L, D) + module subroutine ldlt_dense_${s}$(A, L, D) ${t}$, intent(in) :: A(:,:) ${t}$, intent(inout) :: L(:,:) ${t}$, intent(inout) :: D(:) @@ -242,7 +242,7 @@ contains #:endfor #:for k, t, s in R_KINDS_TYPES - module subroutine factorize_ldlt_csr_${s}$(A, L, D) + module subroutine ldlt_csr_${s}$(A, L, D) type(CSR_${s}$_type), intent(in) :: A type(CSR_${s}$_type), intent(inout) :: L ${t}$, intent(inout) :: D(:) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index c9d6041db..efa8ca307 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -167,31 +167,31 @@ module stdlib_linalg_iterative_solvers end interface public :: solve_backward_triangular - interface factorize_ldlt + interface ldlt #:for k, t, s in R_KINDS_TYPES - module subroutine factorize_ldlt_dense_${s}$(A, L, D) + module subroutine ldlt_dense_${s}$(A, L, D) ${t}$, intent(in) :: A(:,:) ${t}$, intent(inout) :: L(:,:) ${t}$, intent(inout) :: D(:) end subroutine - module subroutine factorize_ldlt_csr_${s}$(A, L, D) + module subroutine ldlt_csr_${s}$(A, L, D) type(CSR_${s}$_type), intent(in) :: A type(CSR_${s}$_type), intent(inout) :: L ${t}$, intent(inout) :: D(:) end subroutine #:endfor end interface - public :: factorize_ldlt + public :: ldlt - interface factorize_ssor + interface ssor #:for k, t, s in R_KINDS_TYPES - module subroutine factorize_ssor_dense_${s}$(A, w, L, D) + module subroutine ssor_dense_${s}$(A, w, L, D) ${t}$, intent(in) :: A(:,:) ${t}$, intent(in) :: w ${t}$, intent(inout) :: L(:,:) ${t}$, intent(inout) :: D(:) end subroutine - module subroutine factorize_ssor_csr_${s}$(A, w, L, D) + module subroutine ssor_csr_${s}$(A, w, L, D) type(CSR_${s}$_type), intent(in) :: A ${t}$, intent(in) :: w type(CSR_${s}$_type), intent(inout) :: L @@ -199,7 +199,7 @@ module stdlib_linalg_iterative_solvers end subroutine #:endfor end interface - public :: factorize_ssor + public :: ssor contains diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index ebc76d7c7..78a83fccf 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -12,10 +12,10 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg implicit none enum, bind(c) - enumerator :: none = 0 - enumerator :: jacobi - enumerator :: ssor - enumerator :: ldlt + enumerator :: pc_none = 0 + enumerator :: pc_jacobi + enumerator :: pc_ssor + enumerator :: pc_ldlt end enum contains @@ -123,7 +123,7 @@ contains maxiter_ = n; if(present(maxiter)) maxiter_ = maxiter restart_ = .true.; if(present(restart)) restart_ = restart tol_ = 1.e-4_${s}$; if(present(tol)) tol_ = tol - precond_ = none ; if(present(precond)) precond_ = precond + precond_ = pc_none; if(present(precond)) precond_ = precond !------------------------- ! internal memory setup ! Preconditionner @@ -134,25 +134,25 @@ contains allocate(diagonal(n),source=zero_${s}$) select case(precond_) - case(jacobi) + case(pc_jacobi) #:if matrix == "dense" diagonal = diag(A) #:else call diag(A,diagonal) #:endif M_%apply => precond_jacobi - case(ssor) + case(pc_ssor) #:if matrix == "dense" diagonal = diag(A) #:else call diag(A,diagonal) #:endif L = A !> copy A structure to L - call factorize_ssor( A , one_${s}$ , L, diagonal ) + call ssor( A , one_${s}$ , L, diagonal ) M_%apply => precond_ldl - case(ldlt) + case(pc_ldlt) L = A !> copy A structure to L - call factorize_ldlt( A , L, diagonal ) + call ldlt( A , L, diagonal ) M_%apply => precond_ldl case default M_%apply => precond_none diff --git a/test/linalg/test_linalg_solve_iterative.fypp b/test/linalg/test_linalg_solve_iterative.fypp index 2859e158c..007c4a3b7 100644 --- a/test/linalg/test_linalg_solve_iterative.fypp +++ b/test/linalg/test_linalg_solve_iterative.fypp @@ -49,13 +49,12 @@ module test_linalg_solve_iterative ${t}$ :: L(5,5), D(5) - call factorize_ldlt(A, L, D) + call ldlt(A, L, D) call check(error, all(abs(L-Lref) Date: Sat, 17 May 2025 18:20:16 +0200 Subject: [PATCH 25/26] change precond ldl name --- src/stdlib_linalg_iterative_solvers_pccg.fypp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pccg.fypp index 78a83fccf..cb21b0efb 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pccg.fypp @@ -149,11 +149,11 @@ contains #:endif L = A !> copy A structure to L call ssor( A , one_${s}$ , L, diagonal ) - M_%apply => precond_ldl + M_%apply => precond_ldlt case(pc_ldlt) L = A !> copy A structure to L call ldlt( A , L, diagonal ) - M_%apply => precond_ldl + M_%apply => precond_ldlt case default M_%apply => precond_none end select @@ -222,7 +222,7 @@ contains ${t}$, intent(in) :: beta y = merge( zero_${s}$, diagonal * x, di_ ) !> inverted diagonal end subroutine - subroutine precond_ldl(x,y,alpha,beta) + subroutine precond_ldlt(x,y,alpha,beta) ${t}$, intent(in) :: x(:) ${t}$, intent(inout) :: y(:) ${t}$, intent(in) :: alpha From b239028d92928653803c4b03882b0e5af92333c8 Mon Sep 17 00:00:00 2001 From: jalvesz Date: Sun, 18 May 2025 16:33:27 +0200 Subject: [PATCH 26/26] rename to pcg --- example/linalg/CMakeLists.txt | 2 +- example/linalg/example_solve_cg.f90 | 2 +- example/linalg/example_solve_custom.f90 | 12 ++++++------ ...mple_solve_pccg.f90 => example_solve_pcg.f90} | 8 ++++---- src/CMakeLists.txt | 2 +- src/stdlib_linalg_iterative_solvers.fypp | 16 ++++++++-------- ... => stdlib_linalg_iterative_solvers_pcg.fypp} | 12 ++++++------ test/linalg/test_linalg_solve_iterative.fypp | 2 +- 8 files changed, 28 insertions(+), 28 deletions(-) rename example/linalg/{example_solve_pccg.f90 => example_solve_pcg.f90} (77%) rename src/{stdlib_linalg_iterative_solvers_pccg.fypp => stdlib_linalg_iterative_solvers_pcg.fypp} (96%) diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 80dd31fb7..2aa9ca2e8 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -42,7 +42,7 @@ ADD_EXAMPLE(solve1) ADD_EXAMPLE(solve2) ADD_EXAMPLE(solve3) ADD_EXAMPLE(solve_cg) -ADD_EXAMPLE(solve_pccg) +ADD_EXAMPLE(solve_pcg) ADD_EXAMPLE(solve_custom) ADD_EXAMPLE(sparse_from_ijv) ADD_EXAMPLE(sparse_data_accessors) diff --git a/example/linalg/example_solve_cg.f90 b/example/linalg/example_solve_cg.f90 index 626cf1184..4ab7261b4 100644 --- a/example/linalg/example_solve_cg.f90 +++ b/example/linalg/example_solve_cg.f90 @@ -1,4 +1,4 @@ -program example_solve_pccg +program example_solve_pcg use stdlib_kinds, only: dp use stdlib_linalg_iterative_solvers, only: solve_cg diff --git a/example/linalg/example_solve_custom.f90 b/example/linalg/example_solve_custom.f90 index e4be3c99e..c7d6c82df 100644 --- a/example/linalg/example_solve_custom.f90 +++ b/example/linalg/example_solve_custom.f90 @@ -3,11 +3,11 @@ module custom_solver use stdlib_sparse use stdlib_linalg_iterative_solvers, only: linop_dp, & solver_workspace_dp, & - solve_pccg_kernel, & - size_wksp_pccg + solve_pcg_kernel, & + size_wksp_pcg implicit none contains - subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) + subroutine solve_pcg_custom(A,b,x,di,tol,maxiter,restart,workspace) type(CSR_dp_type), intent(in) :: A real(dp), intent(in) :: b(:) real(dp), intent(inout) :: x(:) @@ -48,7 +48,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) else allocate( workspace_ ) end if - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = 0.d0 ) + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pcg) , source = 0.d0 ) workspace_%callback => my_logger !------------------------- ! Jacobi preconditionner factorization @@ -56,7 +56,7 @@ subroutine solve_pccg_custom(A,b,x,di,tol,maxiter,restart,workspace) where(abs(diagonal)>epsilon(0.d0)) diagonal = 1._dp/diagonal !------------------------- ! main call to the solver - call solve_pccg_kernel(op,M,b,x,tol_,maxiter_,workspace_) + call solve_pcg_kernel(op,M,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup @@ -126,7 +126,7 @@ program example_solve_custom dirichlet = .false._1 dirichlet([1,5]) = .true._1 - call solve_pccg_custom(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) + call solve_pcg_custom(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0] end program example_solve_custom \ No newline at end of file diff --git a/example/linalg/example_solve_pccg.f90 b/example/linalg/example_solve_pcg.f90 similarity index 77% rename from example/linalg/example_solve_pccg.f90 rename to example/linalg/example_solve_pcg.f90 index 6c5b47b40..289fc0f05 100644 --- a/example/linalg/example_solve_pccg.f90 +++ b/example/linalg/example_solve_pcg.f90 @@ -1,7 +1,7 @@ -program example_solve_pccg +program example_solve_pcg use stdlib_kinds, only: dp use stdlib_sparse - use stdlib_linalg_iterative_solvers, only: solve_pccg + use stdlib_linalg_iterative_solvers, only: solve_pcg type(CSR_dp_type) :: laplacian_csr type(COO_dp_type) :: COO @@ -23,10 +23,10 @@ program example_solve_pccg dirichlet = .false._1 dirichlet([1,5]) = .true._1 - call solve_pccg(laplacian, load, x, tol=1.d-6, di=dirichlet) + call solve_pcg(laplacian, load, x, tol=1.d-6, di=dirichlet) print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0] x = 0._dp - call solve_pccg(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) + call solve_pcg(laplacian_csr, load, x, tol=1.d-6, di=dirichlet) print *, x !> solution: [0.0, 2.5, 5.0, 2.5, 0.0] end program \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4c767c8b9..9c409bb38 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -95,7 +95,7 @@ set(cppFiles stdlib_linalg_iterative_aux.fypp stdlib_linalg_iterative_solvers.fypp stdlib_linalg_iterative_solvers_cg.fypp - stdlib_linalg_iterative_solvers_pccg.fypp + stdlib_linalg_iterative_solvers_pcg.fypp ) add_subdirectory(blas) diff --git a/src/stdlib_linalg_iterative_solvers.fypp b/src/stdlib_linalg_iterative_solvers.fypp index efa8ca307..8fcdf86e0 100644 --- a/src/stdlib_linalg_iterative_solvers.fypp +++ b/src/stdlib_linalg_iterative_solvers.fypp @@ -14,9 +14,9 @@ module stdlib_linalg_iterative_solvers !> details The size of the workspace is defined by the number of vectors used in the iterative solver. enum, bind(c) enumerator :: size_wksp_cg = 3 - enumerator :: size_wksp_pccg = 4 + enumerator :: size_wksp_pcg = 4 end enum - public :: size_wksp_cg, size_wksp_pccg + public :: size_wksp_cg, size_wksp_pcg !> linear operator class for the iterative solvers #:for k, t, s in R_KINDS_TYPES @@ -93,9 +93,9 @@ module stdlib_linalg_iterative_solvers end interface public :: solve_cg - interface solve_pccg_kernel + interface solve_pcg_kernel #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace) + module subroutine solve_pcg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:) @@ -106,12 +106,12 @@ module stdlib_linalg_iterative_solvers end subroutine #:endfor end interface - public :: solve_pccg_kernel + public :: solve_pcg_kernel - interface solve_pccg + interface solve_pcg #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) + module subroutine solve_pcg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) #:if matrix == "dense" ${t}$, intent(in) :: A(:,:) #:else @@ -130,7 +130,7 @@ module stdlib_linalg_iterative_solvers #:endfor #:endfor end interface - public :: solve_pccg + public :: solve_pcg interface solve_forward_triangular !> Solve the system L x = b, where L is a strictly lower triangular matrix (diagonal assumed = 1). diff --git a/src/stdlib_linalg_iterative_solvers_pccg.fypp b/src/stdlib_linalg_iterative_solvers_pcg.fypp similarity index 96% rename from src/stdlib_linalg_iterative_solvers_pccg.fypp rename to src/stdlib_linalg_iterative_solvers_pcg.fypp index 78a83fccf..fb7b48a0a 100644 --- a/src/stdlib_linalg_iterative_solvers_pccg.fypp +++ b/src/stdlib_linalg_iterative_solvers_pcg.fypp @@ -4,7 +4,7 @@ #:set MATRIX_TYPES = ["dense", "CSR"] #:set RANKS = range(1, 2+1) -submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg +submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pcg use stdlib_kinds use stdlib_sparse use stdlib_constants @@ -21,7 +21,7 @@ submodule(stdlib_linalg_iterative_solvers) stdlib_linalg_iterative_pccg contains #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace) + module subroutine solve_pcg_kernel_${s}$(A,M,b,x,tol,maxiter,workspace) class(linop_${s}$), intent(in) :: A class(linop_${s}$), intent(in) :: M !> preconditionner ${t}$, intent(in) :: b(:), tol @@ -85,7 +85,7 @@ contains #:for matrix in MATRIX_TYPES #:for k, t, s in R_KINDS_TYPES - module subroutine solve_pccg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) + module subroutine solve_pcg_${matrix}$_${s}$(A,b,x,di,tol,maxiter,restart,precond,M,workspace) #:if matrix == "dense" use stdlib_linalg, only: diag ${t}$, intent(in) :: A(:,:) @@ -175,12 +175,12 @@ contains else allocate( workspace_ ) end if - if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pccg) , source = zero_${s}$ ) + if(.not.allocated(workspace_%tmp)) allocate( workspace_%tmp(n,size_wksp_pcg) , source = zero_${s}$ ) !------------------------- ! main call to the solver if(restart_) x = zero_${s}$ x = merge( b, x, di_ ) !> copy dirichlet load conditions encoded in B and indicated by di - call solve_pccg_kernel(op,M_,b,x,tol_,maxiter_,workspace_) + call solve_pcg_kernel(op,M_,b,x,tol_,maxiter_,workspace_) !------------------------- ! internal memory cleanup @@ -236,4 +236,4 @@ contains #:endfor #:endfor -end submodule stdlib_linalg_iterative_pccg \ No newline at end of file +end submodule stdlib_linalg_iterative_pcg \ No newline at end of file diff --git a/test/linalg/test_linalg_solve_iterative.fypp b/test/linalg/test_linalg_solve_iterative.fypp index 007c4a3b7..7b5db0a62 100644 --- a/test/linalg/test_linalg_solve_iterative.fypp +++ b/test/linalg/test_linalg_solve_iterative.fypp @@ -97,7 +97,7 @@ module test_linalg_solve_iterative #:endfor end subroutine test_solve_ldlt - ! TODO: add tests for solve_cg, solve_pccg, etc. + ! TODO: add tests for solve_cg, solve_pcg, etc. end module test_linalg_solve_iterative