! complex_real.h: compile real or complex versions from the same source: ************************************************************************* ** TEST PROGRAM FOR various types of sparse solvers ************************************************************************* module gm_data ! real*8 :: V_morgan = 1.0 ! real*8 :: V_morgan = 41.0 real*8 :: V_morgan = 0.0 end module program validation * use gm_data implicit none integer lda, ldstrt, lwork parameter (lda = 2000, ldstrt = 60) parameter (lwork = ldstrt**2 + ldstrt*(lda+5) + 5*lda + 1) * integer i, j, n, nx, 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 a(lda,lda), work(lwork), xtrue(lda) 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: nx = 42 ! including 2 boundary points dx = 1.0/(nx-1) n = nx**2 ! n = nx write(*,*) "Matrix size nx=", nx, "N = nx^2 = ", n if (n.gt.lda) then write(*,*) 'You are asking for a too large matrix' goto 100 endif ! initial guess: do j = 1,n work(j) = ZERO xtrue(j) = 0.0 enddo !---------------------------------------- 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 work(n+j) = ZERO enddo ! one or two delta functions: ! work(n+n/2+1) = -1/dx**2 work(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_rhs_2d(nx, work(n+1)) endif !! !----------------------------------- ***************************************** ** 1st or 2cd order Richardson iterative solution of A x = b: ***************************************** * tol = 1.e-9 call dr_gmres(n, Ax_sub_2d, work(1), work(n+1), 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,*) work(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(work(j)), real(xtrue(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. use gm_data implicit none integer n, nx complex*16 x(nx,nx), Ax(nx,nx) real*8 Delta integer i, j ! 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) = -(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: Ax(:,1) = 0.0 Ax(:,nx) = 0.0 Ax(1,:) = 0.0 Ax(nx,:) = 0.0 return end !----------------------------------------------------------------------- subroutine set_rhs_2d(nx, b) ! set rhs for the 2D Morgan convection-diffusion problem implicit none integer nx complex*16 b(nx,nx) integer i, j b = 0.0 do i = 2, nx-1 do j = 2, nx-1 b(i,j) = -41**2 enddo enddo return end !----------------------------------------------------------------------- subroutine riter1(n,A,lda,x,b, info, rinfo) ! solve A x = b using 1st order Richardson iteration ! on input, x contains the initial guess. implicit none integer n, lda complex*16 A(lda,lda) 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, err integer iter, niter, maxiter dt=0.4 maxiter = 100 write(*,*) "b=", b do iter = 1, maxiter ! perform the matrix vector product Ax = A * x call zgemv('N',n,n,ONE,A,lda,x,1,ZERO,Ax,1) res = Ax - b ! write(*,*) "x =", x ! write(*,*) "res =", res xn = x + dt*res x = xn err=maxval(abs(res)) write(*,*) "err =", err if(err .le. 1.e-10) goto 20 enddo 20 niter = iter info(2)=niter rinfo(2)=err return end subroutine riter1 !----------------------------------------------------------------------- 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, lda, 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