Initial commit from Polytechnique Montreal

1 files changed, 181 insertions, 0 deletions

diff --git a/Trivac/src/FLDMRA.f90 b/Trivac/src/FLDMRA.f90
new file mode 100755
index 0000000..ba42c07
--- /dev/null
+++ b/Trivac/src/FLDMRA.f90
@@ -0,0 +1,181 @@
+!
+!-----------------------------------------------------------------------
+!
+!Purpose:
+! GMRES(m) linear equation solver.
+!
+!Copyright:
+! Copyright (c) 2019 Ecole Polytechnique de Montreal
+! this library is free software; you can redistribute it and/or
+! modify it under the terms of the gnu lesser general public
+! license as published by the free software foundation; either
+! version 2.1 of the license, or (at your option) any later version
+!
+!Author(s): A. Hebert
+!
+!Reference:
+! based on a Matlab script by C. T. Kelley, July 10, 1994.
+!
+!Parameters: input
+! B       fixed source
+! atv     function pointer for the matrix-vector product returning
+!         X+M*(B-A*X) where X is the unknown and B is the source.
+!         The format for atv is "Y=atv(X,B,n,...)"
+! n       order of matrix A
+! ertol   iteration convergence criterion
+! nstart  restarts the GMRES method every nstart iterations
+! maxit   maximum number of GMRES iterations.
+! impx    print parameter: =0: no print; =1: minimum printing.
+! iptrk   L_TRACK pointer to the tracking information
+! ipsys   L_SYSTEM pointer to system matrices
+! ipflux  L_FLUX pointer to the solution
+!
+!Parameters: input/output
+! X       initial estimate / solution of the linear system.
+!
+!Parameters: output
+! iter    actual number of iterations
+!
+!----------------------------------------------------------------------------
+!
+subroutine FLDMRA(B,atv,n,ertol,nstart,maxit,impx,iptrk,ipsys,ipflux,X,iter)
+  use GANLIB
+  implicit real(kind=8) (a-h,o-z)
+  !----
+  !  subroutine arguments
+  !----
+  real(kind=8), dimension(n), intent(in) :: B
+  integer, intent(in) :: nstart,maxit,impx
+  real(kind=8), intent(in) :: ertol
+  interface
+    function atv(X,B,n,iptrk,ipsys,ipflux) result(Y)
+      use GANLIB
+      integer, intent(in) :: n
+      real(kind=8), dimension(n), intent(in) :: X, B
+      real(kind=8), dimension(n) :: Y
+      type(c_ptr) iptrk,ipsys,ipflux
+    end function atv
+  end interface
+  real(kind=8), dimension(n), intent(inout) :: X
+  integer, intent(out) :: iter
+  type(c_ptr) iptrk,ipsys,ipflux
+  !----
+  !  local variables
+  !----
+  integer, parameter :: iunout=6
+  !----
+  !  allocatable arays
+  !----
+  real(kind=8), allocatable, dimension(:) :: r,qq,g,c,s
+  real(kind=8), allocatable, dimension(:,:) :: v,h
+  !----
+  !  scratch storage allocation
+  !----
+  allocate(v(n,nstart+1),g(nstart+1),h(nstart+1,nstart+1), &
+  c(nstart+1),s(nstart+1))
+  !----
+  !  global GMRES(m) iteration.
+  !----
+  allocate(r(n),qq(n))
+  eps1=ertol*sqrt(dot_product(B(:n),B(:n)))
+  rho=1.0d10
+  iter=0
+  do while((rho > eps1).and.(iter < maxit))
+    r(:)=atv(X,B,n,iptrk,ipsys,ipflux)-X(:)
+    rho=sqrt(dot_product(r(:n),r(:n)))
+    !----
+    !  test for termination on entry
+    !----
+    if(rho < eps1) then
+       deallocate(qq,r)
+       go to 100
+    endif
+    !
+    g(:nstart+1)=0.0d0
+    h(:nstart,:nstart)=0.0d0
+    v(:n,:nstart+1)=0.0d0
+    c(:nstart+1)=0.0d0
+    s(:nstart+1)=0.0d0
+    g(1)=rho
+    v(:n,1)=r(:n)/rho
+    !----
+    !  gmres(1) iteration
+    !----
+    k=0
+    do while((rho > eps1).and.(k < nstart).and.(iter < maxit))
+      k=k+1
+      iter=iter+1
+      if(impx > 2) write(iunout,200) iter,rho,eps1
+      qq(:n)=0.0d0
+      r(:)=atv(v(:,k),qq,n,iptrk,ipsys,ipflux)
+      v(:n,k+1)=v(:n,k)-r(:n)
+      !----
+      !  modified Gram-Schmidt
+      !----
+      do j=1,k
+        hr=dot_product(v(:n,j),v(:n,k+1))
+        h(j,k)=hr
+        v(:n,k+1)=v(:n,k+1)-hr*v(:n,j)
+      enddo
+      h(k+1,k)=sqrt(dot_product(v(:n,k+1),v(:n,k+1)))
+      !----
+      !  reorthogonalize
+      !----
+      do j=1,k
+        hr=dot_product(v(:n,j),v(:n,k+1))
+        h(j,k)=h(j,k)+hr
+        v(:n,k+1)=v(:n,k+1)-hr*v(:n,j)
+      enddo
+      h(k+1,k)=sqrt(dot_product(v(:n,k+1),v(:n,k+1)))
+      !----
+      !  watch out for happy breakdown 
+      !----
+      if(h(k+1,k) /= 0.0) then
+        v(:n,k+1)=v(:n,k+1)/h(k+1,k)
+      endif
+      !----
+      !  form and store the information for the new Givens rotation
+      !----
+      do i=1,k-1
+        w1=c(i)*h(i,k)-s(i)*h(i+1,k)
+        w2=s(i)*h(i,k)+c(i)*h(i+1,k)
+        h(i,k)=w1
+        h(i+1,k)=w2
+      enddo
+      znu=sqrt(h(k,k)**2+h(k+1,k)**2)
+      if(znu /= 0.0) then
+        c(k)=h(k,k)/znu
+        s(k)=-h(k+1,k)/znu
+        h(k,k)=c(k)*h(k,k)-s(k)*h(k+1,k)
+        h(k+1,k)=0.0d0
+        w1=c(k)*g(k)-s(k)*g(k+1)
+        w2=s(k)*g(k)+c(k)*g(k+1)
+        g(k)=w1
+        g(k+1)=w2
+      endif
+      !----
+      !  update the residual norm
+      !----
+      rho=abs(g(k+1))
+    enddo
+    !----
+    !  at this point either k > nstart or rho < eps1.
+    !  it's time to compute x and cycle.
+    !----
+    h(:k,k+1)=g(:k)
+    call ALSBD(k,1,h,ier,nstart+1)
+    if(ier /= 0) call XABORT('FLDMRA: singular matrix.')
+    do i=1,n
+      X(i)=X(i)+dot_product(v(i,:k),h(:k,k+1))
+    enddo
+  enddo
+  deallocate(qq,r)
+  !----
+  !  scratch storage deallocation
+  !----
+  100 deallocate(s,c,h,g,v)
+  return
+  !
+  200 format(24h FLDMRA: outer iteration,i4,10h  L2 norm=,1p,e11.4, &
+  6h eps1=,e11.4)
+end subroutine FLDMRA