Add a customizable logger facility, change linop matvec to apply

Author: jalvesz

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

src/stdlib_linalg_iterative_solvers.fypp

@@ -10,22 +10,26 @@ 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
 
     #: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
 

src/stdlib_linalg_iterative_solvers_cg.fypp

@@ -24,30 +24,31 @@ 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)
 
             P = R
 
             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

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