program qfepde
!
!***********************************************************************
!
!! QFEPDE is a simple finite element program 
!
!
!  Discussion:
!
!    QFEPDE solves a single partial differential equation on a rectangular 
!    domain using the finite element method with quadratic elements.
!
!  Modified:
!
!    08 April 2002
!
!  Parameters:
!
!  a       =  matrix containing coefficients defined by the finite element
!           method.  The solution of a*u = f is the solution to the
!           original partial differential equation.
!
!  area    =  for each element, the area of the triangle.
!
!  f       =  vector which at first stores the right hand side of the
!           set of linear equations to be solved, and then the solution.
!
!  indx    =  for a given node, the number of the associated variable,
!           if any, or zero (for boundary nodes).
!
!  maxban  =  maximum bandwidth allowed in matrix.
!
!  maxvar  =  maximum number of variables (unknowns) allowed, which is equal
!           to the number of nodes, minus the boundary nodes.
!
!  maxnod  =  maximum number of nodes allowed.
!
!  maxtri  =  maximum number of triangles allowed.
!
!  nnt     =  number of nodes per triangle, fixed at 6 for this program.
!  node    =  an array which contains the 6 nodes that define a given triangle.
!  np      =  number of nodes for this problem.
!  nt      =  number of triangles for this problem.
!  nu      =  actual number of variables for this problem.
!  nx      =  number of elements in x direction.
!  ny      =  number of elements in y direction.
!  we      =  weights associated with approximate integration of the error.
!  wk      =  work space needed by the linear equation solver.
!  wq      =  weights associated with approximate integration of the
!           finite element coefficients.
!  xe      =  x-coordinates of points used for approximate integration
!           of the error.
!  xl      =  x-coordinate of left of region.
!  xn      =  a vector which contains the x coordinate of each node.
!  xq      =  x-coordinates of points associated with approximate integration
!           of the finite element coefficients.
!  xr      =  x-coordinate of right of region.
!  yb      =  y-coordinate of bottom of region.
!  ye      =  y-coordinates of points associated with the
!             approximate integration of the error.
!  yn      =  a vector which contains the y coordinate of each node.
!  yq      =  y-coordinates of points associated with the approximate
!             integration of the finite element coefficients.
!  yt      =  y-coordinate of top of region.
!
  implicit none
!
  integer, parameter :: maxvar = 81
  integer, parameter :: maxban = 25
  integer, parameter :: maxtri = 250
  integer, parameter :: maxnod = 121
  integer, parameter :: nnt = 6
  integer, parameter :: nq = 3
  integer, parameter :: ne = 13
!
  real a(maxvar,maxban)
  real area(maxtri)
  real eh1
  real el2
  real f(maxvar)
  integer ib
  integer ierr
  integer indx(maxvar)
  integer node(maxtri,nnt)
  integer np
  integer nt
  integer nu
  integer nx
  integer ny
  real we(ne)
  real wk(maxvar)
  real wq(nq)
  real xe(maxtri,ne)
  real, parameter :: xl = 0.0
  real xn(maxnod)
  real xq(maxtri,nq)
  real, parameter :: xr = 1.0
  real, parameter :: yb = 0.0
  real ye(maxtri,ne)
  real yn(maxnod)
  real yq(maxtri,nq)
  real, parameter :: yt = 1.0
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QFEPDE'
  write ( *, '(a)' ) '  A finite element program.'
  write ( *, '(a)' ) '  The region is a 2 dimensional rectangle.'
  write ( *, '(a)' ) '  Quadratic basis elements are used.'

  nx = 4
  ny = 4

  call geometry ( xn, yn, area, node, indx, nt, nu, np, nnt, &
    xl, xr, yb, yt, nx, ny, maxtri, maxvar, maxnod )

  call setint ( nq, wq, xq, yq, maxtri, node, nt, xn, yn, maxnod, nnt )

  call setint ( ne, we, xe, ye, maxtri, node, nt, xn, yn, maxnod, nnt )

  call assemble ( a, f, ib, xn, yn, xq, yq, wq, nq, area, node, indx, &
    nnt, nu, nt, maxtri, maxvar, maxban, maxnod )

  call bepp ( maxvar, nu, ib, ib, a, 1, f, wk, ierr, maxban )

  if ( ierr /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QFEPDE - Fatal error!'
    write ( *, '(a)' ) '  The matrix is singular!'
    stop
  end if

  call prisol ( f, indx, maxtri, nnt, node, np, nt, xn, yn, maxvar, maxnod )

  call eror ( el2, eh1, we, xe, ye, ne, area, node, indx, xn, yn, &
    f, nt, nnt, maxtri, maxvar, maxnod )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QFEPDE:'
  write ( *, '(a,g14.6)' ) '  L2 solution error = ', el2
  write ( *, '(a,g14.6)' ) '  H1 solution error = ', eh1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'QFEPDE'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine setint ( n, w, x, y, maxtri, node, nt, xn, yn, maxnod, nnt )
!
!***********************************************************************
!
!! SETINT sets up a Gauss quadrature rule.
!
!
!  Discussion:
!
!    SETINT sets the weights and the coordinates of
!    the points involved in the gaussian quadrature.
!
!    The gaussian quadrature is performed twice in the
!    code; the first time with a three-points formula
!    in order to calculate the integrals that appear
!    in the stiffness matrix, the second time with a
!    thirteen-points formula in order to check the
!    accuracy of the results in the error analysis.
!
!    The parameter N, passed to setint as an argument,
!    tells which of the two quadratures is being performed.
!
  implicit none
!
  integer maxnod
  integer maxtri
  integer n
  integer nnt
!
  integer i
  integer ii
  integer iii
  integer ip1 
  integer ip2
  integer ip3
  integer it
  integer node(maxtri,nnt)
  integer nt
  real w(n)
  real x(maxtri,n)
  real x1
  real x2
  real x3
  real xn(maxnod)
  real y(maxtri,n)
  real y1
  real y2
  real y3
  real yn(maxnod)
  real z1
  real z2
  real z3
  real z4
  real z5
  real z6
  real z7
!
  if ( n == 3 ) then

    w(1:3) = 1.0 / 3.0

    do it = 1, nt
      ip1 = node(it,1)
      ip2 = node(it,2)
      ip3 = node(it,3)
      x1 = xn(ip1)
      x2 = xn(ip2)
      x3 = xn(ip3)
      y1 = yn(ip1)
      y2 = yn(ip2)
      y3 = yn(ip3)
      x(it,1) = 0.5*(x1+x2)
      x(it,2) = 0.5*(x2+x3)
      x(it,3) = 0.5*(x3+x1)
      y(it,1) = 0.5*(y1+y2)
      y(it,2) = 0.5*(y2+y3)
      y(it,3) = 0.5*(y3+y1)
    end do

  else

    do i = 1, 3
      w(i) = 0.175615257433204
      ii = i+3
      w(ii) = 0.053347235608839
      ii = i+6
      w(ii) = 0.077113760890257
      iii = ii+3
      w(iii) = w(ii)
    end do

    w(13) = 1.0 - sum ( w(1:12) )

    z1 = 0.479308067841923
    z2 = 0.260345966079038
    z3 = 0.869739794195568
    z4 = 0.065130102902216
    z5 = 0.638444188569809
    z6 = 0.312865496004875
    z7 = 0.048690315425316

    do it = 1, nt
      ip1 = node(it,1)
      ip2 = node(it,2)
      ip3 = node(it,3)
      x1 = xn(ip1)
      x2 = xn(ip2)
      x3 = xn(ip3)
      y1 = yn(ip1)
      y2 = yn(ip2)
      y3 = yn(ip3)
      x(it, 1) = z1*x1+z2*x2+z2*x3
      y(it, 1) = z1*y1+z2*y2+z2*y3
      x(it, 2) = z2*x1+z1*x2+z2*x3
      y(it, 2) = z2*y1+z1*y2+z2*y3
      x(it, 3) = z2*x1+z2*x2+z1*x3
      y(it, 3) = z2*y1+z2*y2+z1*y3
      x(it, 4) = z3*x1+z4*x2+z4*x3
      y(it, 4) = z3*y1+z4*y2+z4*y3
      x(it, 5) = z4*x1+z3*x2+z4*x3
      y(it, 5) = z4*y1+z3*y2+z4*y3
      x(it, 6) = z4*x1+z4*x2+z3*x3
      y(it, 6) = z4*y1+z4*y2+z3*y3
      x(it, 7) = z5*x1+z6*x2+z7*x3
      y(it, 7) = z5*y1+z6*y2+z7*y3
      x(it, 8) = z5*x1+z7*x2+z6*x3
      y(it, 8) = z5*y1+z7*y2+z6*y3
      x(it, 9) = z6*x1+z5*x2+z7*x3
      y(it, 9) = z6*y1+z5*y2+z7*y3
      x(it,10) = z6*x1+z7*x2+z5*x3
      y(it,10) = z6*y1+z7*y2+z5*y3
      x(it,11) = z7*x1+z5*x2+z6*x3
      y(it,11) = z7*y1+z5*y2+z6*y3
      x(it,12) = z7*x1+z6*x2+z5*x3
      y(it,12) = z7*y1+z6*y2+z5*y3
      x(it,13) = (x1+x2+x3)/3.0
      y(it,13) = (y1+y2+y3)/3.0
    end do

  end if

  return
end
subroutine geometry ( xn, yn, area, node, indx, nt, nu, np, nnt, &
  xl, xr, yb, yt, nx, ny, maxtri, maxvar, maxnod )
!
!******************************************************************************
!
!! GEOMETRY carries out a variety of computations related to the geometry.
!
!
!  Discussion:
!
!    The routine
!
!  (1)  set the origin of the coordinate at the lower lefthand corner of
!       the rectangular domain. the x-coordinate matches the horizontal edge
!       of the domain. the y-coordinate matches the vertical edge of the domain.
!  (2)  divide the rectangular domain into (nx-1)*(ny-1) small sub-rectangular
!       domains.
!  (3)  each sub-rectangular domain is divided into 2 triangular elements by a
!       diagonal line from the lower lefthand corner to the upper righthand
!       corner.
!  (4)  the triangular element's no. counts by row-wise (horizontal) order
!       from the left bottom of the domain in the sub-rectangular base.
!       in each sub-rectangular domain, it counts the left upper triangular
!       element first,then the right lower triangular element second.
!       there are 2*(nx-1)*(ny-1) triangular elements totally.
!  (5)  each triangular element contains six nodes,three at the tips of the
!       triangle and the other three at the middle points of each side of the
!       the triangule seperately.
!       there are (nx+(nx-1))*(ny+(ny-1)) nodes totally.
!  (6)  the node's no. counts by row-wise order from the left bottom of the
!       domain.
!  (7)  the node's order in each triangular is followings:
!       (i)  the first 3 nodes are the tip's nodes, count from the smallest
!            no. of the nodes (lower left corner), then follow the cw (clock
!            wise) direction for the upper triangular elements and the ccw
!            (counter clock wise) direction for the lower triangular elements.
!       (ii) the nodes of 4,5 and 6 are the three middle nodes, count from
!            the smallest no. of the nodes, then follow the direction as (i)
!            mentioned.
!  (8)  calculate the coordinates of each node.
!  (9)  count the no. of unknown variables and all node's index no.
!  (10) calculate the area of each triangular element.
!
!  Parameters:
!
!   grdx   :external function to calculate the real node's x-coordinate
!   grdy   :external function to calculate the real node's y-coordinate
!   nxm    :no. of triangular elements in x-direction
!   nym    :no. of triangular elements in y-direction
!   hx     :normalized side length of triangle in x-direction
!   hy     :normalized side length of triangle in y-direction
!   xx     :normalized node's x-coordinate
!   yy     :normalized node's y-coordonate
!   nx2    :no. of nodes in x-direction
!   ny2    :no. of nodes in y-direction
!   ip     :node's no.
!   it1    :upper triangular element's no.
!   it2    :lower triangular element's no.
!
  implicit none
!
  integer maxnod
  integer maxtri
  integer maxvar
  integer nnt
!
  real area(maxtri)
  real grdx
  real grdy
  real hx
  real hy
  integer i
  integer ic
  integer indx(maxvar)
  integer ip
  integer ip1
  integer ip2
  integer ip3
  integer it
  integer it1
  integer it2
  integer ix
  integer iy
  integer j
  integer jc
  integer node(maxtri,nnt)
  integer np
  integer nt
  integer nu
  integer nx
  integer nx2
  integer nxm
  integer ny
  integer ny2
  integer nym
  real tol
  real x
  real x1
  real x2
  real x3
  real xl
  real xn(maxnod)
  real xr
  real xx
  real y
  real y1
  real y2
  real y3
  real yb
  real yn(maxnod)
  real yt
  real yy
!
!  GRDX and GRDY are "statement functions", which can be used to
!  adjust the arrangement of nodes in the X and Y directions.
!
  grdx(x) = x
  grdy(y) = y

  nxm = nx-1
  hx = (xr-xl) / real(nxm)
  nym = ny-1
  hy = (yt-yb) / real(nym)

  yy = -hy / 2.0

  nx2 = nx+nxm
  ny2 = ny+nym

  nu = 0
  ip = 0
  it1 = -1

  do jc = 1, ny2

    yy = yy + hy/2.0

    xx = -hx / 2.0

    iy = mod ( jc, 2 )

    do ic = 1, nx2

      xx = xx + hx/2.0
      ix = mod ( ic, 2 )

      ip = ip+1

      xn(ip) = grdx(xx)
      yn(ip) = grdy(yy)

      if ( ix == 1 .and. iy == 1 ) then

        if ( ic /= nx2 .and. jc /= ny2 ) then

          it1 = it1 + 2
          it2 = it1 + 1

          node(it1,1) = ip
          node(it1,2) = ip+nx2+nx2
          node(it1,3) = ip+nx2+nx2+2
          node(it1,4) = ip+nx2
          node(it1,5) = ip+nx2+nx2+1
          node(it1,6) = ip+nx2+1
          node(it2,1) = ip
          node(it2,2) = ip+2
          node(it2,3) = ip+nx2+nx2+2
          node(it2,4) = ip+1
          node(it2,5) = ip+nx2+2
          node(it2,6) = ip+nx2+1

        end if

      end if

      tol = 0.0002

      if ( ic /= 1 .and. ic /= nx2 .and. jc /= 1 .and. jc /= ny2 ) then

        if ( nu >= maxvar ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'GEOMETRY - Fatal error!'
          write ( *, '(a)' ) '  Too many variables!'
          stop
        end if

        nu = nu+1
        indx(ip) = nu

      else

        indx(ip) = 0

      end if

    end do

  end do

  np = ip
  nt = it2

  do it = 1, nt
    ip1 = node(it,1)
    ip2 = node(it,2)
    ip3 = node(it,3)
    x1 = xn(ip1)
    x2 = xn(ip2)
    x3 = xn(ip3)
    y1 = yn(ip1)
    y2 = yn(ip2)
    y3 = yn(ip3)
    area(it) = 0.5 * abs ( y1*(x2-x3) + y2*(x3-x1) + y3*(x1-x2) )
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '       Triangle #           Node #'
  write ( *, '(a)' ) ' '

  do i = 1, nt
    write ( *, '(5x,i6,8x,6i6)' ) i, node(i,1:6)
  end do

  return
end
subroutine assemble ( a, f, ib, xn, yn, xq, yq, wq, nq, area, node, indx, &
  nnt, nu, nt, maxtri, maxvar, maxban, maxnod )
!
!***********************************************************************
!
!! ASSEMBLE assembles the stiffness matrix and load vector.
!
!
!  Discussion:
!
!    The stiffness matrix is stored in a matrix with as many rows as
!    there are degrees of freedom and as many columns as the max
!    bandwidth of the system of equations.
!
  implicit none
!
  integer maxban
  integer maxnod
  integer maxtri
  integer maxvar
  integer nnt
  integer nq
!
  real a(maxvar,maxban)
  real aa
  real aij
  real ar
  real area(maxtri)
  real bb
  real bbb
  real bbx
  real bby
  real bx
  real by
  real f(maxvar)
  integer i
  integer ib
  integer ib1
  integer ibt
  integer ii
  integer ij
  integer in
  integer indx(maxvar)
  integer inn
  integer ip
  integer ipp
  integer iq
  integer it
  integer j
  integer node(maxtri,nnt)
  integer nt
  integer nu
  real rhs
  real ubc
  real wq(nq)
  real x
  real xn(maxnod)
  real xq(maxtri,nq)
  real xx
  real y
  real yn(maxnod)
  real yq(maxtri,nq)
  real yy
!
  ib = 0
  do it = 1, nt

    do in = 1, nnt

      ip = node(it,in)

      i = indx(ip)

      if ( i /= 0 ) then

        do inn = 1, nnt
          ipp = node(it,inn)
          j = indx(ipp)
          ij = j-i
          ib = max ( ib, ij )
        end do

      end if

    end do
  end do

  ib1 = ib+1
  ibt = ib+ib1

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'ASSEMBLE:'
  write ( *, '(a,i6)' ) '  The matrix bandwidth is ', ibt

  if ( ibt > maxban ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'ASSEMBLE:'
    write ( *, '(a)' ) '  Fatal error!'
    write ( *, '(a)' ) '  Exceeding the maximum bandwidth of ', maxban
    stop
  end if

  f(1:nu) = 0.0

  a(1:nu,1:ibt) = 0.0

  do it = 1, nt

    do iq = 1, nq

      x = xq(it,iq)
      y = yq(it,iq)
      ar = area(it) * wq(iq)

      do in = 1, nnt

        ip = node(it,in)
        call qbf ( x, y, it, in, bb, bx, by, xn, yn, node, maxtri )

        i = indx(ip)

        if ( i /= 0 ) then
          ii = ib1-i

          f(i) = f(i) + rhs(x,y) * bb * ar

          do inn = 1, nnt

            ipp = node(it,inn)
            call qbf(x,y,it,inn,bbb,bbx,bby,xn,yn,node,maxtri)
            j = indx(ipp)

            aij = (bx*bbx+by*bby) * aa(x,y)

            if ( j /= 0 ) then
              ij = ii+j
              a(i,ij) = a(i,ij) + aij*ar
            else
              xx = xn(ipp)
              yy = yn(ipp)
              f(i) = f(i) - ar * aij*ubc(xx,yy)
            end if

          end do

        end if
      end do

    end do

  end do

  return
end
subroutine qbf ( x, y, it, in, b, dbdx, dbdy, xn, yn, node, maxtri )
!
!***********************************************************************
!
!! QBF evaluates the shape functions and their derivatives.
!
!
!  Discussion:
!
!    The coordinate transformation between global (X,Y) coordinates and 
!    local (T,S) coordinates is accomplished by an affine mapping.
!
  implicit none
!
  integer maxtri
!
  real b
  real c
  real d
  real dbdx
  real dbdy
  integer i1
  integer i2
  integer i3
  integer in
  integer in1
  integer in2
  integer in3
  integer inn
  integer it
  integer j1
  integer j2
  integer j3
  integer node(maxtri,*)
  real s
  real t
  real x
  real xn(*)
  real y
  real yn(*)
!
  if ( in <= 3 ) then

    in1 = in
    in2 = mod(in,3)+1
    in3 = mod(in+1,3)+1
    i1 = node(it,in1)
    i2 = node(it,in2)
    i3 = node(it,in3)
    d = (xn(i2)-xn(i1))*(yn(i3)-yn(i1))-(xn(i3)-xn(i1))*(yn(i2)-yn(i1))
    t = 1.0+((yn(i2)-yn(i3))*(x-xn(i1))+(xn(i3)-xn(i2))*(y-yn(i1)))/d
    b = t*(2.0*t-1.0)
    dbdx = (yn(i2)-yn(i3))*(4.0*t-1.0)/d
    dbdy = (xn(i3)-xn(i2))*(4.0*t-1.0)/d

  else

    inn = in-3
    in1 = inn
    in2 = mod(inn,3)+1
    in3 = mod(inn+1,3)+1
    i1 = node(it,in1)
    i2 = node(it,in2)
    i3 = node(it,in3)
    j1 = i2
    j2 = i3
    j3 = i1
    d = (xn(i2)-xn(i1))*(yn(i3)-yn(i1))-(xn(i3)-xn(i1))*(yn(i2)-yn(i1))
    c = (xn(j2)-xn(j1))*(yn(j3)-yn(j1))-(xn(j3)-xn(j1))*(yn(j2)-yn(j1))
    t = 1.0 + ((yn(i2)-yn(i3))*(x-xn(i1))+(xn(i3)-xn(i2))*(y-yn(i1)))/d
    s = 1.0 + ((yn(j2)-yn(j3))*(x-xn(j1))+(xn(j3)-xn(j2))*(y-yn(j1)))/c

    b = 4.0 * s * t
    dbdx = 4.0 * (t*(yn(j2)-yn(j3))/c + s*(yn(i2)-yn(i3))/d)
    dbdy = 4.0 * (t*(xn(j3)-xn(j2))/c + s*(xn(i3)-xn(i2))/d)

  end if

  return
end
function aa ( x, y )
!
!***********************************************************************
!
!! AA evaluates the coefficient function A(X,Y) in the differential equation.
!
!
!  Parameters:
!
!    Input, real X, Y, the point at which the function is to be evaluated.
!
!    Output, real AA, the value of the function at (X,Y).
!
  implicit none
!
  real aa
  real x
  real y
!
  aa = 1.0

  return
end
function uex ( x, y )
!
!***********************************************************************
!
!! UEX returns the exact solution of the differential equation.
!
!
!  Parameters:
!
!    Input, real X, Y, the point at which the function is to be evaluated.
!
!    Output, real UEX, the value of the exact solution at (X,Y).
!
  implicit none
!
  real uex
  real x
  real y
!
  uex = sin ( x ) * cos ( y )

  return
end
function uexx ( x, y )
!
!***********************************************************************
!
!! UEXX returns the X derivative of the exact solution.
!
!
!  Parameters:
!
!    Input, real X, Y, the point at which the function is to be evaluated.
!
!    Output, real UEXX, the value of the X derivative of the
!    exact solution at (X,Y).
!
  implicit none
!
  real uexx
  real x
  real y
!
  uexx = cos ( x ) * cos ( y )

  return
end
function uexy ( x, y )
!
!***********************************************************************
!
!! UEXY returns the Y derivative of the exact solution.
!
!
!  Parameters:
!
!    Input, real X, Y, the point at which the function is to be evaluated.
!
!    Output, real UEXY, the value of the Y derivative of the
!    exact solution at (X,Y).
!
  implicit none
!
  real uexy
  real x
  real y
!
  uexy = - sin ( x ) * sin ( y )

  return
end
function rhs ( x, y )
!
!***********************************************************************
!
!! RHS returns the value of the right hand side of the PDE.
!
!
!  Parameters:
!
!    Input, real X, Y, the point at which the function is to be evaluated.
!
!    Output, real RHS, the value of the right hand side of the
!    differential equation at (X,Y).
!
  implicit none
!
  real rhs
  real tol
  real x
  real y
!
  rhs = - 2.0 * sin ( x ) * cos ( y )

  return
end
function ubc ( x, y )
!
!***********************************************************************
!
!! UBC returns the Dirichlet boundary condition at any point on the boundary.
!
  implicit none
!
  real ubc
  real uex
  real x
  real y
!
  ubc = uex ( x, y )

  return
end
subroutine bepp ( maxvar, n, m1, m2, a, nrhs, b, wk, ierr, maxban )
!
!***********************************************************************
!
!! BEPP performs Gauss elimination of a banded linear system.
!
  implicit none
!
  integer maxban
  integer maxvar
  integer nrhs
!
  real a(maxvar,maxban)
  real b(maxvar,nrhs)
  real fff
  integer i
  integer ierr
  integer ii
  integer j
  integer je
  integer js
  integer k
  integer l
  integer m
  integer m1
  integer m2
  integer n
  real p
  real q
  real sicri
  real wk(maxvar)
  real x
!
  sicri = 1.0e-16
  ierr = 0
  m = m1+m2+1

  do i = 1, n
    p = abs(a(i,1))
    do j = 2, m
      q = abs(a(i,j))
      if ( q > p)p = q
    end do
    fff = abs(p)
    if ( fff < sicri ) then
      ierr = 1
      return
    end if
    wk(i) = p
  end do

  do i = 1, m1
    je = m2+i
    js = m1+1-i
    do j = 1, je
      a(i,j) = a(i,js+j)
    end do
    je = je+1
    do j = je, m
      a(i,j) = 0.0
    end do
  end do

  l = m1

  do k = 1, n-1

    if ( l < n ) then
      l = l+1
    end if

    i = k
    p = abs(a(k,1)) / wk(k)

    do j = k+1, l
      q = abs(a(j,1)) / wk(j)
      if ( q > p ) then
        p = q
        i = j
      end if
    end do

    fff = abs ( p )

    if ( fff < sicri ) then
      ierr = 1
      return
    end if

    if ( i /= k ) then
      do j = 1, m
        x = a(k,j)
        a(k,j) = a(i,j)
        a(i,j) = x
      end do
      do j = 1, nrhs
        x = b(k,j)
        b(k,j) = b(i,j)
        b(i,j) = x
      end do
    end if

    do i = k+1, l
      x = a(i,1)/a(k,1)
      do j = 2, m
        a(i,j-1) = a(i,j) - x * a(k,j)
      end do
      a(i,m) = 0.0
      do j = 1, nrhs
        b(i,j) = b(i,j) - x*b(k,j)
      end do
    end do

  end do

  fff = abs(a(n,1))

  if ( fff < sicri ) then
    ierr = 1
    return
  end if

  do j = 1, nrhs

    l = 1

    do i = n, 1, -1

      fff = abs ( a(i,1) )

      if ( fff < sicri ) then
        ierr = 1
        return
      end if

      x = b(i,j)

      if ( i /= n ) then
        do k = 2, l
          x = x - a(i,k) * b(i-1+k,j)
        end do
      end if
 
      b(i,j) = x / a(i,1)

      if ( l < m ) then
        l = l + 1
      end if

    end do

  end do

  return
end
subroutine eror ( el2, eh1, we, xe, ye, ne, area, node, indx, xn, yn, &
  f, nt, nnt, maxtri, maxvar, maxnod )

!**********************************************************************
!
!! EROR calculates L2 and H1 errors in the solution. 
!
!
!  l2 norm of the difference (u-utrue) and
!     root mean square deviation of the gradient
!     (ux-utruex)**2+(uy-utruey)**2.
!     they are named as el2 and eh1, respectively.
!
  implicit none
!
  integer maxnod
  integer maxtri
  integer maxvar
  integer ne
  integer nnt
!
  real ar
  real area(maxtri)
  real bb
  real bx
  real by
  real eh1
  real el2
  real f(maxvar)
  integer i
  integer in
  integer indx(maxvar)
  integer ip
  integer iq
  integer it
  integer node(maxtri,nnt)
  integer nt
  real uex
  real uexy
  real uexx
  real uh
  real uhx
  real uhy
  real we(ne)
  real x
  real x1
  real xe(maxtri,ne)
  real xn(maxnod)
  real y
  real y1
  real ye(maxtri,ne)
  real yn(maxnod)
!
  el2 = 0.0
  eh1 = 0.0

  do it = 1, nt

    do iq = 1, ne

      ar = area(it) * we(iq)

      x = xe(it,iq)
      y = ye(it,iq)
      uh = 0.0
      uhx = 0.0
      uhy = 0.0
      do  in = 1, nnt
        ip = node(it,in)

        call qbf(x,y,it,in,bb,bx,by,xn,yn,node,maxtri)

        x1 = xn(ip)
        y1 = yn(ip)
        i = indx(ip)

        if ( i > 0 ) then
          uh = uh+bb*f(i)
          uhx = uhx+bx*f(i)
          uhy = uhy+by*f(i)
        else
          uh = uh+bb*uex(x1,y1)
          uhx = uhx+bx*uex(x1,y1)
          uhy = uhy+by*uex(x1,y1)
        end if
      end do

      el2 = el2+ar*(uh-uex(x,y))**2

      eh1 = eh1+ar*((uhx-uexx(x,y))**2+(uhy-uexy(x,y))**2)

    end do

  end do

  el2 = sqrt ( el2 )
  eh1 = sqrt ( eh1 )

  return
end
subroutine prisol ( f, indx, maxtri, nnt, node, np, nt, xn, yn, maxvar, &
  maxnod )
!
!***********************************************************************
!
!! PRISOL prints the solution as a table of 5 columns.
!
!
!  x,y,a(x,y),u(calcul.),utrue.
!     where:
!     x,y:coordinates of node.
!     a(x,y):variable thermal conductivity at node.
!     u: calculated solution at node.
!     utrue:analytical sol. for comparison.
!
  implicit none
!
  integer maxnod
  integer maxtri
  integer maxvar
  integer nnt
!
  real a1
  real aa
  real f(maxvar)
  integer i
  integer indx(maxvar)
  integer ip
  integer node(maxtri,nnt)
  integer np
  integer nt
  real uex
  real uh
  real utrue
  real x
  real xn(maxnod)
  real y
  real yn(maxnod)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) &
    '       x             y            a(x,y)             u           utrue'

  write ( *, '(a)' ) ' '

  do ip = 1, np

    i = indx(ip)
    x = xn(ip)
    y = yn(ip)
    a1 = aa(x,y)

    if ( i > 0 ) then
      uh = f(i)
    else 
      uh = uex ( x, y )
    end if

    utrue = uex ( x, y )

    write ( *, '(5g14.6)' ) x, y, a1, uh, utrue

  end do

  return
end