program hw2_runbcg
c
c Compile e.g. with 
c
c   gfortran -Wall -o hw2_runbcg -Wall   hw2_runbcg.f  bcg2.f
c
c This depends on the reference implementation of bicg-stab (bcg2.f) from the 
c web-page of H.v.d.Vorst (one of the authors of BiCGstab)
c   http://www.staff.science.uu.nl/~vorst102/software.html
c Use "simple version" of BiCGstab(ell)
c 
      parameter (n=1000)
      parameter (l=1)
      double precision x(n)
      double precision b(n)

      double precision z(n)
      double precision tol

      integer mxmv
      double precision work(n,(3+3+2*(l+1)))
      integer info, i
      external bicgstab2
      external matvec_prec
      external hw2_2c
      external hw2_2b
      external precond
      external tridiag_solve
      double precision resnorm

      tol=1E-10
      mxmv=300
      
      
      do i = 1, n
c right-hand side
         b(i) = 3.0
c initial guess
         x(i) = -0.4
      end do


c Run BiCGStab
      call bicgstab2 (1, l, n, x, b, hw2_2b, 1, tol,
     +      'rel', mxmv, work, n*10000, info)

c print the output
      write (*,*) 'info=',info, ' total matvecs=', mxmv

c Compute and print the residual ||Ax-b||
      call hw2_2b(n,x,z)
      resnorm=0.0
      do  i = 1, n
         resnorm=resnorm+ (b(i)-z(i))**2.0
      end do
      resnorm=sqrt(resnorm)
      write (*,*) '||Ax-b||', resnorm

      stop

      end


c     Delivers the matrix vector product y=A*x for A given in HW2.2c
      subroutine hw2_2c (n,x,y)
      integer n
      integer i
      parameter (pi=3.141592653589793)
      double precision x(n), y(n) 
      double precision h

      h=1.0/n

c     compute the matrix vector product corresponding to tridiagonal matrix
      do i = 1, n

         y(i)=(-5.0/2.0)*x(i)/(h*h)
     +          +10.0*sin(pi*(1.0*i)/(1.0*n))*x(n-i+1)
         if (i .GT. 1) then
            y(i)=y(i)+(4.0/3.0)*x(i-1)/(h*h)
         end if
         if (i .GT. 2) then
            y(i)=y(i)-(1.0/12.0)*x(i-2)/(h*h)
         end if

         if (i .LT. n) then
            y(i)=y(i)+(4.0/3.0)*x(i+1)/(h*h)
         end if
         if (i .LT. n-1) then
            y(i)=y(i)-(1.0/12.0)*x(i+2)/(h*h)
         end if
      end do

      end

c   Delivers the matrix vector product for homework 2.2b
      subroutine hw2_2b(n,x,y)
      integer n
      double precision x(n), y(n) 
      integer i

c     compute the matrix vector product corresponding to tridiagonal matrix
      do i = 1, n
         y(i)=(-1.0)*x(i)
         if (i .GT. 1) then
            y(i)=y(i)+0.1*x(i-1)
         end if
         if (i .LT. n) then
            y(i)=y(i)-0.1*x(i+1)
         end if
      end do

      end





      subroutine tridiag_solve(n,alpha,beta,b,x)
c
c     Solve a linear system for a symmetric tridiagonal with constant diagonal
c     entries.
c
c        T= alpha beta  0     0      0 
c           beta  alpha beta   
c           0       .     .   .
c                      beta   alpha beta
c           0       0     0   beta  alpha
c    The input b contains the right-hand side 
c    The output x contains the soluction Ax=b
c    
c
      integer n
      double precision alpha, beta
      double precision b(n), x(n)
      double precision vp(n)
      double precision bp(n)
      double precision m
      integer i

      bp(1)=alpha
      vp(1)=b(1)



c    
      do i = 2,n
         m = beta/bp(i-1)
         bp(i) = alpha - m*beta
         vp(i) = b(i) - m*vp(i-1)
      end do 

      x(n) = vp(n)/bp(n)

c   Backward subs
      do i = n-1, 1, -1
         x(i) = (vp(i) - beta*x(i+1))/bp(i)
      end do 
      
      end




c     Preconditioned version of matrix vector product
c     
      subroutine matvec_prec (n,x,y)
      integer n
      double precision x(n), y(n) 
      external precond
      external matvec

      call hw2_2c(n,x,y)
      call precond(n,y)

      end
      

      subroutine precond (n,z)
      integer n
      double precision z(n), b(n)
      double precision alpha, beta
      double precision h
      integer i
      external tridiag_solve 
      do i = 1, n
         b(i)=z(i)
      end do
      h=1.0/n
      alpha=10.0
      beta=10.0
      call tridiag_solve(n, alpha, beta, b, z)

      end