# 1 "dr_gmres_test.F" ************************************************************************* ** TEST PROGRAM FOR various types of sparse solvers ************************************************************************* ! ! Solving A x = b ! Will also denote x with Phi. ! # 1 "complex_real.h" 1 ! complex_real.h: compile real or complex versions from the same source: # 10 "dr_gmres_test.F" 2 module gm_data ! real*8 :: V_morgan = 1.0 ! real*8 :: V_morgan = 41.0 real*8 :: V_morgan = 0.0 ! nx = # of grid points in 1 direction, including 2 bdy points. integer, parameter :: nx = 42 integer, parameter :: n = nx**2 complex*16 Phi(nx, nx) complex*16 Phi1d(n) ! equivalent 1-D version of Phi equivalence (Phi,Phi1d) complex*16 b(nx, nx) complex*16 b1d(n) ! equivalent 1-D version of b equivalence (b,b1d) logical :: bc_grad = .false. ! = .false. b.c. on Phi (Dirichlet b.c.'s) ! = .true. b.c. on grad(Phi) (Neumann b.c.'s) end module program validation * use gm_data implicit none * integer i, j, m integer revcom, colx, coly, colz, nbscal integer irc(5), icntl(8), info(3) * integer matvec, precondLeft, precondRight, dotProd parameter (matvec=1, precondLeft=2, precondRight=3, dotProd=4) * integer nout * complex*16 xtrue(n) real*8 cntl(5), rinfo(2) real*8 dx, x, pi real*8 tol ! relative error tolerance * complex*16 ZERO, ONE parameter (ZERO = (0.0d0, 0.0d0), ONE = (1.0d0, 0.0d0)) ! external Ax_sub_1d ! subroutine to evaluate A.x for 1d Poisson Eq. external Ax_sub_2d ! subroutine to evaluate A.x * logical :: problem_delta = .false. pi = abs(acos(-1.0)) * *************************************************************** ** Generate the test matrix A and set the right-hand side ** in positions (n+1) to 2n of the array work. ** The right-hand side is chosen such that the exact solution ** is the vector of all ones. *************************************************************** * write(*,*) '***********************************************' write(*,*) "This code is an example of an algorithm I'll call" write(*,*) 'optimized dynamic relaxation. Another possible' write(*,*) 'name is recycled GMRES(2), which is like standard' write(*,*) 'GMRES(2) but with a modified restarting method that' write(*,*) 'recycles previous search directions rather than just' write(*,*) & 'throwing them away. This is equivalent to LGMRES(1,1)' write(*,*) '(Loose GMRES), and can also be though of as an' write(*,*) 'optimized 2cd Order Richardson iteration, or' write(*,*) 'a modified form of 2cd order Chebyshev iteration.' write(*,*) 'Results are written in output files' write(*,*) 'fort.* and sol_zTestgmres.' write(*,*) '***********************************************' write(*,*) ! do 2-D Poisson equation eventually, start with 1-d: dx = 1.0/(nx-1) write(*,*) "Matrix size nx=", nx, "N = nx^2 = ", n !---------------------------------------- if(problem_delta) then print *, "can only do Morgan problem at present" stop ! ! or give an explicit b and ignore xtrue: do j = 1,n b1d(n+j) = ZERO enddo ! one or two delta functions: ! b1d(n+n/2+1) = -1/dx**2 b1d(n+n*0.75+1) = -1/dx**2 else !---------------------------------------- ! RHS for the Morgan problem: print *, "doing the Morgan problem with V_morgan = ", V_morgan call set_x0_rhs_2d endif !! !----------------------------------- ***************************************** ** 1st or 2cd order Richardson iterative solution of A x = b: ***************************************** * tol = 1.e-9 call dr_gmres(n, Ax_sub_2d, Phi, b, tol, info, & rinfo, nx) * ******************************* * dump the solution on a file ******************************* * nout = 11 open(nout,FILE='sol_zTestgmres',STATUS='unknown') ! write(nout,*) 'Restart : ', m write(nout,*) 'info(1) = ',info(1),' info(2) = ',info(2) write(nout,*) 'rinfo(1) = ',rinfo(1),' rinfo(2) = ',rinfo(2) write(nout,*) 'Optimal workspace = ', info(3) write(nout,*) 'Solution : ' do j=1,n write(nout,*) Phi1d(j) enddo write(nout,*) ' ' * write(12,*) "# x Phi(x), Phi_true(x)" write(12,*) "0 0.0 0.0" do j=1,n write(12,*) j, real(Phi1d(j)) enddo write(12,*) n+1," 0.0 0.0" write(*,*) "# iterations =", info(2), "relative residual err=", & rinfo(2) write(*,*) "average convergence rate /iteration", & rinfo(2)**(1./info(2)) 100 continue * stop end !----------------------------------------------------------------------- subroutine Ax_sub_1d(n, x, Ax, nx) ! evaluate Matrix-Vector product A x and return the result in Ax implicit none integer n, nx complex*16 x(n), Ax(n) real*8 Delta integer i ! for a 1-D Poisson equation, on the domain [0,1] with a step size 1/(n-1), ! with x=0 boundary conditions. Delta = 1.0/(nx-1) ! skip boundary points at i=1, n do i=2,n-1 Ax(i) = (x(i+1)-2*x(i)+x(i-1))/Delta**2 enddo return end !----------------------------------------------------------------------- subroutine Ax_sub_2d(n, x, Ax, nx) ! evaluate Matrix-Vector product A x and return the result in Ax ! for the 2-D Poisson equation or the isotropic Diffusion problem. use gm_data, only: v_morgan, b, bc_grad implicit none integer n, nx complex*16 x(nx,nx), Ax(nx,nx) real*8 Delta, unit, dt, D integer i, j ! For 2-D Poisson problem: dt = 1 unit = 0 D = 1 ! For isotropic diffusion problem: dt = 200 unit = 1 D = 0.01 ! for a 2-D Poisson equation, on the domain [0,1] x [0,1] ! with a step size 1/(n-1), with x=0 boundary conditions. Delta = 1.0/(nx-1) ! V = 0.0 ! V = 1.0 ! V = 41.0 ! V = 41.0**2 ! skip boundary points at i=1, n do i = 2 , nx-1 do j = 2, nx-1 Ax(i,j) = unit*x(i,j) & + dt*D*( -(x(i+1,j )-2*x(i,j)+x(i-1,j ))/Delta**2 & -(x(i ,j+1)-2*x(i,j)+x(i ,j-1))/Delta**2 ) & +V_morgan*(x(i+1,j)-x(i-1,j))/(2*Delta) enddo enddo ! set results at boundaries if(.not. bc_grad) then !Dirichlet B.C.'s, (Phi=b) Ax(:,1) = x(:,1) Ax(:,nx) = x(:,nx) Ax(1,:) = x(1,:) Ax(nx,:) = x(nx,:) else ! Neumann B.C.'s (grad(Phi)=0 on bdy's) Ax(:,1) = Ax(:,2) Ax(:,nx) = Ax(:,nx-1) Ax(1,:) = Ax(2,:) Ax(nx,:) = Ax(nx-1,:) endif return end subroutine Ax_sub_2d !----------------------------------------------------------------------- subroutine set_x0_rhs_2d ! set initial x and rhs for the 2D Morgan convection-diffusion problem: use gm_data, only: nx, b, Phi, bc_grad implicit none integer i, j ! initialize the boundary values of b first: if(bc_grad) then b = 0 else b = 50.0 ! boundary values for x ! ! b = 0.0 ! boundary values for x endif do i = 2, nx-1 do j = 2, nx-1 ! Morgan problem: ! b(i,j) = -41**2 ! Put in spatially varying temperature for diffusion problem: b(i,j) = -41**2*sin(i*2*3.14159/(1.0*nx)) enddo enddo ! initial guess at solution x0 (or Phi0): Phi = 0 ! boundary conditions: ! initial x0 must satisfy the b.c.'s too: if(.not. bc_grad) then ! Dirichlet B.C.'s: Phi(:,1) = b(:,1) Phi(:,nx) = b(:,nx) Phi(1,:) = b(1,:) Phi(nx,:) = b(nx,:) else ! Neumann B.C.'s (grad(Phi)=0 on bdy's) Phi(:,1) = Phi(:,2) Phi(:,nx) = Phi(:,nx-1) Phi(1,:) = Phi(2,:) Phi(nx,:) = Phi(nx-1,:) endif return end !----------------------------------------------------------------------- subroutine riter2(n, Ax_sub, x, b, info, rinfo, nx) ! solve A x = b using 2cd order Richardson iteration ! on input, x contains the initial guess. implicit none integer n, nx complex*16 x(n), b(n) integer info(2) real*8 rinfo(2) complex*16 Ax(n), xp(n), xn(n), res(n) complex*16, parameter :: ONE = (1.0d0, 0.0d0) complex*16, parameter :: ZERO = (0.0d0, 0.0d0) real*8 dt, nu, err, pi, kmin, alpha, beta1, nudt_half integer iter, niter, maxiter external Ax_sub dt=0.5*2.0 nu=2*2*3.14159/n * 1.1 maxiter = 1000 pi = 2*abs(acos(0.0)) ! Calculated optimal iteration parameters: kmin=2*pi/(n+1) alpha = 1-2*kmin beta1 = 1/4.+alpha/4+sqrt(alpha)/2 nudt_half=(1-alpha)/(1+alpha) dt = sqrt(beta1*(1+nudt_half)) ! dt = (1+2*kmin)/(1+kmin)/2*(1-kmin+sqrt(1-2*kmin)) ! dt = sqrt(dt) nu=nudt_half*2/dt*0.33 ! must be an error in the above calculation of the optimal parameters, because ! I empirically find a better fit is obtained if I multiply this nu by 0.33... write(*,*) "alpha, beta1 =", alpha, beta1 write(*,*) "dt, nu =", dt, nu write(*,*) "predicted error reduction / iteration =", 1/(1+nu*dt) write(*,*) "b=", b write(15,*) "# iteration #, err" do iter = 1, maxiter ! perform the matrix vector product Ax = A * x call Ax_sub(n, x, Ax, nx) res = Ax - b ! write(*,*) "x =", x ! write(*,*) "res =", res xn = x + (x-xp)*(1-nu*dt/2)/(1+nu*dt/2) + dt**2*res/(1+nu*dt/2) xp = x x = xn err=maxval(abs(res)) write(15,*) iter, err if(err .le. 1.e-10) goto 20 enddo 20 niter = iter info(2)=niter rinfo(2)=err return end subroutine riter2