c-------------------------------------------------------------:
      program biharmonic
c-------------------------------------------------------------:
c----
c---- u    + 2*u    + u    + alpha*(u  +u  ) + beta*u = f(x,y)
c----  xxxx     xxyy   yyyy         xx  yy
c----      
c-------------------------------------------------------------:

      implicit none
      
      external approxit,error,f1,f2,f3,u,rand,gettime
      integer m,n,wn,idf,nmax,mmax
      parameter (mmax=9999)
      parameter (nmax=9999)
      parameter (idf=mmax+2)
      parameter (wn=3*mmax+4*nmax+4*nmax+2*mmax+0.5*(nmax+1)**2+19+100)
      integer iflag,itcg,lw,i,j,problem,k,nstop
      double precision a,b,c,d,alpha,beta,tol
      double precision bda(nmax)
      double precision bdb(nmax)
      double precision bdc(mmax)
      double precision bdd(mmax)
      double precision f(mmax+2,nmax+2)
      double precision w(wn)
      double precision x,y,c0,t,emax,el2,erel
      double precision approxit,f1,f2,f3,u,rand,gettime
 100  format(I5,' ',E16.8,' ',E16.8,' ',E16.8,' ',E16.8)
c-------------------------------------------------------------:
c---- Initialize problem
c-------------------------------------------------------------:
c---- Problem & constants
c-------------------------------------------------------------:
      problem = 2
      c0 = 30
      nstop = 20
      write(*,*)'Solving problem',problem
c-------------------------------------------------------------:
c---- Problem size
c-------------------------------------------------------------:
      m = 199
      n = 199
c-------------------------------------------------------------:
c---- Domain
c-------------------------------------------------------------:
      a = -1
      b = 1
      c = -1
      d = 1
c-------------------------------------------------------------:
c---- Constant of the equation
c-------------------------------------------------------------:
      alpha = 1
      beta  = 1
      lw = wn
c-------------------------------------------------------------:
c---- Write outut file
c-------------------------------------------------------------:
      write(*,*)'Writing ouptufile.'
      open(unit=10,file='bih_output')
c-------------------------------------------------------------:
c---- Loop
c-------------------------------------------------------------:
      do k = 3000,4000
         m = k*2+1
         n = k*2+1
C         c0 = k * 5
C         if (k .eq. 0) c0 =1
c-------------------------------------------------------------:
c---- Option iflag
c-------------------------------------------------------------:
      iflag = 4
c-------------------------------------------------------------:
c---- Exterior derivatives values, f3
c-------------------------------------------------------------:
         do j = 1,n
            x = a
            y = c+j*(d-c)/(n+1)
            bda(j) = f3(problem,x,y,4,c0)
            x = b
            y = c+j*(d-c)/(n+1)
            bdb(j) = f3(problem,x,y,2,c0)
         enddo
         do i = 1,m
            x = a+i*(b-a)/(m+1)
            y = c
            bdc(i) = f3(problem,x,y,1,c0)
            x = a+i*(b-a)/(m+1)
            y = d
            bdd(i) = f3(problem,x,y,3,c0)
         enddo
c-------------------------------------------------------------:
c---- Function values, f1
c-------------------------------------------------------------:
         do i = 2,m+1
            do j = 2,n+1
               x =a+(i-1)*(b-a)/(m+1)
               y =c+(j-1)*(d-c)/(n+1)
               f(i,j) = f1(problem,x,y,alpha,beta,c0)
            enddo
         enddo
c-------------------------------------------------------------:
c---- Boundary values, f2
c-------------------------------------------------------------:
         do i = 1,m+2 
            x = a+(i-1)*(b-a)/(m+1)
            y = c
            f(i,1)   = f2(problem,x,y,c0)
            x = a+(i-1)*(b-a)/(m+1)
            y = d
            f(i,n+2) = f2(problem,x,y,c0)
         end do

         do j = 1,n+2
            x = a
            y = c+(j-1)*(d-c)/(n+1)
            f(1,j)   = f2(problem,x,y,c0)
            y = b
            x = c+(j-1)*(d-c)/(n+1)
            f(m+2,j) = f2(problem,x,y,c0)
         end do     
c-------------------------------------------------------------:
c---- Solve problem
c-------------------------------------------------------------:
         t = gettime()
         call dbihar(a,b,m,bda,bdb,bdc,bdd,c,d,n,f,idf,
     +        alpha,beta,iflag,tol,itcg,w,lw)
         t = gettime()-t
         call error(problem,f,mmax,nmax,m,n,a,b,c,d,c0,emax,el2,erel)
         write(*,100)int(n),emax,el2,erel,t
         write(10,100)int(c0),emax,el2,erel,t
c          do i = 1,m+2
c             do j = 1, n+2
c                x = a + (i-1)*(b-a)/(m+1)
c                y = c + (j-1)*(d-c)/(n+1)
c                write(10,*)i,j,f(i,j)
c                write(10,*)i,j,u(problem,x,y,c0)
c             end do
c          end do
c-------------------------------------------------------------:
c---- End loop
c-------------------------------------------------------------:
          IF(erel .LT. 0.001) stop
      end do
c-------------------------------------------------------------:
c---- Info output of bihar
c-------------------------------------------------------------:
      write(*,*)'Output info of dbihar:'
      write(*,*)'iflag =',iflag
      write(*,*)'tol   =',tol
      write(*,*)'itcg  =',itcg
      write(*,*)'lw    =',lw
      write(*,*)'Ndof  =',m*n
      write(*,*)'Time  =',t
c-------------------------------------------------------------:
c---- Norm-2 Error
c-------------------------------------------------------------:
      write(*,*)'L2    ='
c-------------------------------------------------------------:
c---- Approximation
c-------------------------------------------------------------:
      write(*,*),'Approximation of solution in arbitrary point'
c      write(*,*),'Write x and y'
c      read(*,*),x,y
      x = rand()
      y = rand()
      write(*,*)x,y
      write(*,*)'Approximation=',approxit(f,mmax,nmax,m,n,a,b,c,d,x,y)
      write(*,*)'Exact        =',u(problem,x,y,c0)
      write(*,*)'Difference   =',approxit(f,mmax,nmax,m,n,a,b,c,d,x,y)-u(problem,x,y,c0)
c-------------------------------------------------------------:
c---- Done
c-------------------------------------------------------------:
      close(10)
      write(*,*)'Done.'
      stop
      end