Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
class(sparse_t), | intent(in) | :: | A | |||
real(kind=wp), | intent(in), | dimension(:) | :: | b | ||
real(kind=wp), | intent(in), | optional | dimension(:) | :: | x0 | |
real(kind=wp), | intent(in), | optional | :: | tol | ||
integer, | intent(in), | optional | :: | maxIts |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module.
function minimumResidual(A,b,x0,tol,maxIts) result(x)
class(sparse_t),intent(in)::A
real(wp),dimension(:),intent(in)::b
real(wp),dimension(:),intent(in),optional::x0
real(wp),intent(in),optional::tol
integer,intent(in),optional::maxIts
real(wp),dimension(:),allocatable::x
real(wp),dimension(:),allocatable::D,r,Ar
real(wp)::alpha
integer::lmaxIts
real(wp)::ltol
real(wp)::rss0,rss
integer::k
lmaxIts = 1000000
ltol = 1.0E-6_wp
if(present(maxIts)) lmaxIts = maxIts
if(present(tol )) ltol = tol
D = A%getDiagonal()
x = b/D
if(present(x0)) x = x0
r = b-matmul(A,x)
rss0 = sqrt(sum(r**2))
call startReport('Minimum Residual Solver')
do k=1,lmaxIts
Ar = matmul(A,r)
alpha = sum(Ar*r)/sum(Ar*Ar)
x = x+alpha*r
r = b-matmul(A,x)
rss = sqrt(sum(r**2))
call solverProgress(k,rss/rss0,ltol,'Minimum Residual Solver')
if(rss/rss0<ltol) exit
end do
end function minimumResidual