module var
        real hi(2)
        real hx
        real hy
        real Nx
        real Ny
        real x0
        real xn
        real y0
        real yn
      end module var

      program main
      use var

        include "link_f90_static.h"
        integer lda,i,j,p,n,shift
        parameter (nxpoint=10,nypoint=10,n=(nypoint-2)*nxpoint,ncoda=nxpoint,lda=n+ncoda)

        integer ijob

        real*8 a(lda,ncoda+2),uu(n),errr,er
        real k1,k2,f,u,g1,g2,g3,g4,x,y
        external lslpb


        open(1,file='OUTPUT2.dat')

        write (1,*)
        write (1,*) '---------------------------------------------------&
     &-----------------------------------------'
      Nx=Nxpoint
      Ny=Nypoint
      x0=1.0
      xn=10.0
      y0=1.0
      yn=10.0
      hx=(xn-x0)/(Nx-1)
      hy=(yn-y0)/(Ny-1)
      hi(1)=1.0
      hi(2)=1.0
      shift=NCODA


        write (1,1) 'n = '
        write (1,*) n
        write (1,*)

      !=======================================================================
!         calculates a()
!=======================================================================
        ! j=2
        ! i=1
        a(ncoda+1,1)=k1(x(1+0.5),y(2.0))*hy/hx+hi(1)*hy+           (k2(x(1.0),y(2.0+0.5))+k2(x(1.0),y(2.0-0.5)))*hx/hy/2
        a(ncoda+1+1,2)=-k1(x(1.0+0.5),y(2.0))*hy/hx
        a(ncoda+1+nx,ncoda+1)=-k2(x(1.0),y(2.0+0.5))*hx/hy/2
        a(ncoda+1,ncoda+2)=hx*hy/2*f(x(1.0),y(2.0))+g1(y(1.0))*hy+    &   k2(x(1.0),y(1.5))*g3(x(1.0))*hx/hy/2
        ! j=2
        ! i=2,...,n-1
        do i=2,nx-1
          a(ncoda+i,1)=(k1(x(i+0.5),y(2.0))+k1(x(i-0.5),y(2.0)))*     &     hy/hx+(k2(x(real(i)),y(2.5))+k2(x(real(i)),y(1.5)))*hx/hy
          a(ncoda+i+1,2)=-k1(x(i+0.5),y(2.0))*hy/hx
          a(ncoda+i+nx,ncoda+1)=-k2(x(real(i)),y(2.0+0.5))*hx/hy
          a(ncoda+i,ncoda+2)=hx*hy*f(x(real(i)),y(2.0))-k2(x(real(i)),y(1.5))*g3(x(1.0))*hx/hy
        end do
        ! j=2
        ! i=n
        a(ncoda+nx,1)=hi(2)*hy+k1(x(nx-0.5),y(2.0))*hy/hx+     &   (k2(x(real(nx)),y(2+0.5))+k2(x(real(nx)),y(2.0-0.5)))*hx/hy/2 !
          !
        a(ncoda+nx+nx,ncoda+1)=-k2(x(real(nx)),y(2.0+0.5))*hx/hy/2                                                                    !
        a(ncoda+nx,ncoda+2)=hx*hy/2*f(x(real(nx)),y(2.0))+k2(x(real(nx)),y(1.5))*g3(x(real(nx)))*hx/hy/2+g2(y(real(nx)))*hy         !
        ! j=3,...,n-2
        ! i=1
        do j=3,ny-2
          a(ncoda+1+(j-2)*nx,1)=k1(x(1+0.5),y(real(j)))*hy/hx+hi(1)*hy+       &     (k2(x(1.0),y(j+0.5))+k2(x(1.0),y(j-0.5)))*hx/hy/2    !
          a(ncoda+1+(j-2)*nx+1,2)=-k1(x(1+0.5),y(real(j)))*hy/hx                                                    !
          a(ncoda+1+(j-2)*nx+nx,ncoda+1)=-k2(x(1.0),y(j+0.5))*hx/hy/2                                      !
          a(ncoda+1+(j-2)*nx,ncoda+2)=g1(y(real(j)))*hy+hx*hy/2*f(x(1.0),y(real(j)))                  !
        end do
        ! j=3,...,n-2
        ! i=2,...,n-1
        do i=2,nx-1
          do j=3,ny-2
            a(ncoda+i+(j-2)*nx,1)=(k1(x(i+0.5),y(real(j)))+k1(x(i-0.5),y(real(j))))*          hy/hx+(k2(x(real(i)),y(j+0.5))+k2(x(real(i)),y(j-0.5)))*hx/hy !
            a(ncoda+i+(j-2)*nx+1,2)=-k1(x(i+0.5),y(real(j)))*hy/hx                                                      !
            a(ncoda+i+(j-2)*nx+nx,ncoda+1)=-k2(x(real(i)),y(j+0.5))*hx/hy                                               !
            a(ncoda+i+(j-2)*nx,ncoda+2)=f(x(real(i)),y(real(j)))*hx*hy                                !
          end do
        end do
        ! j=3,...,n-2
        ! i=n
        do j=3,ny-2
          a(ncoda+nx+(j-2)*nx,1)=hi(2)*hy+k1(x(nx-0.5),y(real(j)))*hy/hx+     &   (k2(x(real(nx)),y(j+0.5))+k2(x(real(nx)),y(j-0.5)))*hx/hy/2                !
          a(ncoda+nx+(j-2)*nx+nx,ncoda+1)=-k2(x(real(nx)),y(j+0.5))*hx/hy/2                                                         !
          a(ncoda+nx+(j-2)*nx,ncoda+2)=g2(y(real(j)))*hy+hx*hy/2*f(x(real(nx)),y(real(j)))                                                !
        end do
        ! j=n-1
        ! i=1
        a(ncoda+1+(nx-1-2)*nx,1)=k1(x(1+0.5),y(ny-1.0))*hy/hx+hi(1)*hy+       &     (k2(x(1.0),y(ny-1.0+0.5))+k2(x(1.0),y(ny-1.0-0.5)))*hx/hy/2         !
        a(ncoda+1+(nx-1-2)*nx+1,2)=-k1(x(1+0.5),y(n-1.0))*hy/hx                                                                                          !
        a(ncoda+1+(nx-1-2)*nx,ncoda+2)=hx*hy/2*f(x(1.0),y(ny-1.0))+g1(y(ny-1.0))*hy+    &   k2(x(1.0),y(ny-0.5))*g4(x(1.0))*hx/hy/2
        ! j=n-1
        ! i=2,...,n-1
        do i=2,nx-1
          a(ncoda+i+(nx-1-2)*nx,1)=(k1(x(i+0.5),y(ny-1.0))+k1(x(i-0.5),y(ny-1.0)))*     &     hy/hx+(k2(x(real(i)),y(ny-1.0+0.5))+k2(x(real(i)),y(ny-1.0-0.5)))*hx/hy !
          a(ncoda+i+(nx-1-2)*nx+1,2)=-k1(x(i+0.5),y(n-1.0))*hy/hx                                                                                     !
          a(ncoda+i+(nx-1-2)*nx,ncoda+2)=k2(x(real(i)),y(ny-0.5))*g4(x(real(i)))*hx/hy/2+hx*hy/2*f(x(real(i)),y(ny-1.0))
        end do
        ! j=n-1
        ! i=n
        a(ncoda+nx+(nx-1-2)*nx,1)=hi(2)*hy+k1(x(nx-0.5),y(n-1.0))*hy/hx+     &   (k2(x(real(nx)),y(n-1.0+0.5))+k2(x(real(nx)),y(n-1.0-0.5)))*hx/hy/2 !
        a(ncoda+nx+(nx-1-2)*nx,ncoda+2)=k2(x(real(nx)),y(ny-0.5))*g4(x(1.0))*hx/hy/2+g2(y(ny-1.0))*hy+hx*hy/2*f(x(real(nx)),y(ny-1.0))
!=======================================================================



      !do j=2,ny-1
      !      do i=1,nx
      !            a(ncoda+I+(j-2)*(nx),1) =-cp(real(i),real(j))
      !            a(ncoda+i+(j-2)*nx,ncoda+2)=gp(real(i),real(j))
      !            a(ncoda+i+(j-2)*nx,ncoda)=-dp(real(i),real(j))
      !            p=I+(j-2)*nx
      !
      !      enddo
      !enddo
      !do j=2,ny-2
      !      do i=1,nx
      !            a(ncoda+i+(j-1)*nx,ncoda+1)=-ep(real(i),real(j))
      !      enddo
      !enddo




        ijob=1
        call dlslpb(n,a,lda,ncoda,ijob,uu)


        write(1,*)
        do j=2,ny-1
          do i=1,nx
            write(1,4) a(ncoda+i+(j-2)*nx,ncoda+2)
          end do
          write(1,*)
        end do

        write(1,*)


        write(1,*)
        write (1,*) '--------------------------------------------------------------------------------------------'

        errr=0.d0
        do j=2,ny-1
          do i=1,nx
              er=abs(u(x(real(i)),y(real(j)))-a(ncoda+i+(j-2)*nx,ncoda+2))
              if (er>errr) errr=er
          end do
        end do

        write(1,1) '||U-Unum|| = '
        write(1,3) errr
        write(*,*) '||U-Unum|| = '
        write(*,*) errr
        read(*,*)

1       format(a13$)
2       format(2i5$)
3       format(e13.3)
4       format(f10.4$)
5       format(i10$)
55      format(i5$)
555     format(i7)
6       format(e10.2$)
7       format(a32$)
8       format(a5$)
9       format(a11$)
10      format(f5.3$)
12      format(e5.3$)
11      format(e16.4$)
      end



      real function x(i)
        use var
        real i
        x=x0+(i-1)*hx
        return
      end

      real function y(j)
        use var
        real j
        y=y0+(j-1)*hy
        return
      end

      real function k1(x,y)
        use var
        real x,y
        k1=1.0
        return
      end

      real function k2(x,y)
        use var
        real x,y
        k2=1.0
        return
      end

      real function f(x,y)
        use var
        real x,y
        !f=-4*(x*y*y+x*x*y) !0.0
        f=0.0
        return
      end

      real function g1(y)
        use var
        real y
        g1=1
        return
      end

      real function g2(y)
        use var
        real y
        g2=1.0
        return
      end

      real function g3(x)
        use var
        real x
        g3=1.0
        return
      end

      real function g4(x)
        use var
        real x
        !g4=100*x*x
        g4=1.0
        return
      end


      real function u(i,j)
        use var
        real i,j
        !u=r(i)*r(i)*z(j)*z(j)*z(j) !r(i)+z(j)
        u=1
        return
      end

      real function Cp(i, j)
        use var
        external k1, k2
        real i, j
        if (i==1.0 ) then
            Cp=-k1(x(i+0.5),y(j))*hy/hx-hi(1)*hy-k2(x(i),y(j+0.5))*hx/hy/2.0-k2(x(i),y(j-0.5))*hx/hy/2.0
            !write(*,*) -r(1.5)*k1(r(1.5),z(j))*hz/hr
            !write(*,*) -hi(1)*r(1.0)*hz
            !write(*,*) -r(1.0)*k2(r(1.0),z(j+0.5))*hr/hz/2.0
            !write(*,*) -r(1.0)*k2(r(1.0),z(j-0.5))*hr/hz/2.0
        else if (i==N) then
            Cp=-hi(2)*hy-k1(x(i-0.5),y(j))*hy/hx-k2(x(i),y(j-0.5))*hx/hy/2.0-k2(x(i),y(j+0.5))*hx/hy/2.0
        else
            Cp=-k1(x(i+0.5),y(j))*hy/hx-k1(x(i-0.5),y(j))*hy/hx-k2(x(i),y(j+0.5))*hx/hy-k2(x(i),y(j-0.5))*hx/hy
        end if
        !write(*,"('Cp='f7.3)") Cp
        return
      end

      real function Dp(i, j)
        use var
        external k1, k2
        real i, j
        Dp=0.0
        if (i>=2.0 .and. i<=Nx-1.0 ) then
            Dp=k1(x(i+0.5),y(j))*hy/hx
            !write(*,*) 11111111
        end if
        if (i==nx) then
            dp=0
        endif
        !if (j>=2.0 .and. j<=Nz-2.0 .and. i>=2.0 .and. i<=Nr-1.0) then
        !    Dp=r(i)*k2(r(i),z(j+0.5))*hr/hz
        !    !write(*,*) 22222222
        !end if
        !write(*,"('Dp='f7.3)") Dp
        return
      end

      real function Ep(i, j)
        use var
        external k1, k2
        real i, j
        Ep=0.0
        if (i==1.0 .or. i==Nx .and. j<=ny-2.0) then
            Ep=k2(x(i),y(j+0.5))*hx/hy/2.0
        end if
        if (i>=2.0 .and. i<=Nx-1.0 .and. j<=ny-2.0) then
            Ep=k2(x(i),y(j+0.5))*hx/hy
            !write(*,*) i, j, Ep
        end if
        !write(*,"('Ep='f7.3)") Ep
        return
      end

      real function Gp(i, j)
        use var
        external k1, k2
        real i, j
        integer flag
        if (j==2.0 .and. i==1.0) then
            Gp=f(x(i),y(j))*hx*hy/2.0+k2(x(i),y(j-0.5))*hx/hy/2.0*g3(x(i))+g1(y(j))*hy !1

        end if
        if (j>=3.0 .and. j<=Ny-2.0 .and. i==1.0) then
            Gp=f(x(i),y(j))*hx*hy/2.0+g1(y(j))*hy !4
        end if
        if (j==Ny-1.0 .and. i==1.0) then
            Gp=f(x(i),y(j))*hx*hy/2.0+k2(x(i),y(j+0.5))*hx/hy/2.0*g4(x(i))+g1(y(j))*hy !7
        end if
        if (j==2.0 .and. i>=2.0 .and. i<=Nx-1.0) then
            Gp=f(x(i),y(j))*hx*hy+k2(x(i),y(j-0.5))*hx/hy*g3(x(i)) !2
        end if
        if (j==2.0 .and. i==Nx) then
            Gp=f(x(i),y(j))*hx*hy/2.0+k2(x(i),y(j-0.5))*hx/hy/2.0*g3(x(i))+g2(y(j))*hy !3
        end if
        if (j==Ny-1.0 .and. i>=2.0 .and. i<=Nx-1.0) then
            Gp=f(x(i),y(j))*hx*hy+k2(x(i),y(j+0.5))*hx/hy*g4(x(i)) !8
        end if
        if (j==Ny-1.0 .and. i==Nx) then
            Gp=f(x(i),y(j))*hx*hy/2.0+k2(x(i),y(j+0.5))*hx/hy/2.0*g4(x(i))+g2(y(j))*hy !9
        end if
        if (j>=3.0 .and. j<=Ny-2.0 .and. i==Nx) then
            Gp= g2(y(j))*hy+f(x(i),y(j))*hx*hy/2.0 !6
        end if
        if (j>=3.0 .and. j<=Ny-2.0 .and. i>=2.0 .and. i<=Nx-1.0) then
            Gp=f(x(i),y(j))*hx*hy !5
        end if
        return
      end