anchor.diff:    end subroutine V3PROD
   
     !======================================================================!
  -  !> @brief VLCPY
  -  !======================================================================!
  -  subroutine VLCPY(A,B,N)
  -    implicit none
  -    !------
  -    integer(kind=kint), intent(in):: N
  -    real(kind=kreal), intent(in):: B(N,N)
  -    real(kind=kreal), intent(out):: A(N,N)
  -    !------
  -    integer(kind=kint):: i
  -
  -    do i = 1, N
  -      A(1:N,i) = B(i,1:N)
  -    enddo
  -  end subroutine VLCPY
  -
  -  !======================================================================!
  -  !> @brief VCOPY
  -  !======================================================================!
  -  subroutine VCOPY(A,C,N)
  -    implicit none
  -    !------
  -    integer(kind=kint), intent(in):: N
  -    real(kind=kreal), intent(in):: A(N)
  -    real(kind=kreal), intent(out):: C(N)
  -    !------
  -    C = A
  -  end subroutine VCOPY
  -
  -  !======================================================================!
  -  !> @brief VPROD
  -  !======================================================================!
  -  subroutine VPROD(A,B,C,N)
  -    implicit none
  -    !------
  -    integer(kind=kint), intent(in):: N
  -    real(kind=kreal), intent(in):: A(N)
  -    real(kind=kreal), intent(in):: B(N)
  -    real(kind=kreal), intent(out):: C(N)
  -    !------
  -    C(1:N) = A(1:N)*B(1:N)
  -  end subroutine VPROD
  -
  -  !======================================================================!
     !> @brief VXPROD
     !======================================================================!
     subroutine VXPROD(Ndeg,Ndegl,Zln,Diag,Zz,N)
  
positive.code:   end subroutine MATINI
  
    !======================================================================!
    !> @brief NUFORM
    !======================================================================!
    subroutine NUFORM(hecMESH,hecMAT,NEQns,NDEg,LEN_colno,NSTop,LEN_dsln,INVp,XLNzr,COLno,DIAg,ZLN,DSLn,RALoc,STAge,Ir)
      implicit none
      !------
      type (HECMWST_LOCAL_MESH), intent(in)::hecMESH
      type (HECMWST_MATRIX), intent(in)::hecMAT
      integer(kind=kint), intent(in):: NEQns
      integer(kind=kint), intent(in):: NDEg
      integer(kind=kint), intent(in):: LEN_colno
      integer(kind=kint), intent(in):: NSTop
      integer(kind=kint), intent(in):: LEN_dsln
      integer(kind=kint), intent(in):: INVp(:)
      integer(kind=kint), intent(in):: XLNzr(:)
      integer(kind=kint), intent(in):: COLno(:)
      real(kind=kreal), pointer :: DIAg(:)
      real(kind=kreal), pointer :: ZLN(:)
      real(kind=kreal), pointer :: DSLn(:)
      integer(kind=kint), intent(inout):: RALoc
      integer(kind=kint), intent(inout):: STAge
      integer(kind=kint), intent(out):: Ir
      !------
      integer(kind=kint):: i
      integer(kind=kint):: idbg
      integer(kind=kint):: ierr
      integer(kind=kint):: j
      integer(kind=kint):: k
      integer(kind=kint):: kk
      integer(kind=kint):: ntotal
      integer(kind=kint):: numnp
      integer(kind=kint):: ndof
      integer(kind=kint):: ndof2
      real(kind=kreal), allocatable :: val(:)
  
      idbg = 0
      numnp = hecMAT%NP
      ndof = hecMESH%N_DOF
      ntotal = numnp*ndof
  
      !*Allocations
      allocate (val(NDEg*NDEg),stat=ierr)
      if ( ierr/=0 ) stop "Allocation error:val"
      write (6,*) "nuform:stage = ", STAge
      kk = 0
      ndof2 = ndof*ndof
  
      do j = 1, numnp
        !*Diagonal
        kk = kk + 1
        call VLCPY(val,hecMAT%D(ndof2*(j-1)+1:ndof2*j),ndof)
        call STAIJ1(0,j,j,val,NEQns,NDEg,LEN_colno,NSTop,LEN_dsln,INVp,XLNzr,COLno,DIAg,ZLN,DSLn,RALoc,STAge,Ir)
  
        do i = 1, ndof
          if ( val((i-1)*ndof+i)<=0 ) write (idbg,*) 'j,j,val:', j, i, val((i-1)*ndof+i)
        enddo
  
        !*Lower
        do k = hecMAT%INDEXL(j-1) + 1, hecMAT%INDEXL(j)
          i = hecMAT%ITEML(k)
          kk = kk + 1
          call VLCPY(val,hecMAT%AL(ndof2*(k-1)+1:ndof2*k),ndof)
          call STAIJ1(0,j,i,val,NEQns,NDEg,LEN_colno,NSTop,LEN_dsln,INVp,XLNzr,COLno,DIAg,ZLN,DSLn,RALoc,STAge,Ir)
        enddo
      enddo
  
      deallocate (val)
    end subroutine NUFORM
  
    !======================================================================!
    !> @brief NUFCT0 performs Cholesky factorization
    !          if(iv(22).eq.0)    normal type
    !          if(iv(22).gt.0)    code generation type
    !======================================================================!
    subroutine NUFCT0(NEQns,NDEg,NSTop,PARent,NCH,XLNzr,COLno,DIAg,ZLN,DSLn,STAge,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: NEQns
      integer(kind=kint), intent(in):: NDEg
      integer(kind=kint), intent(in):: NSTop
      integer(kind=kint), intent(in):: XLNzr(:)
      integer(kind=kint), intent(in):: COLno(:)
      integer(kind=kint), intent(in):: PARent(:)
      integer(kind=kint), intent(inout):: NCH(:)
      real(kind=kreal), intent(inout):: DIAg(:)
      real(kind=kreal), intent(inout):: ZLN(:)
      real(kind=kreal), intent(inout):: DSLn(:)
      integer(kind=kint), intent(inout):: STAge
      integer(kind=kint), intent(out):: Ir
      !------
      integer(kind=kint):: irr
      integer(kind=kint):: ndeg2
      integer(kind=kint):: ndegl
      integer(kind=kint), allocatable :: INDx(:)
      real(kind=kreal), allocatable :: TEMp(:)
  
      if ( STAge/=20 ) then
        print *, '*********Setting Stage 40!*********'
        Ir = 40
        return
      else
        Ir = 0
      endif
  
      allocate (TEMp(NDEg*NDEg*NEQns),stat=IRR)
      if ( IRR/=0 ) then
        write (*,*) '##Error : Not enough memory'
        call HECMW_ABORT(HECMW_COMM_GET_COMM())
        !stop
      endif
      allocate (INDx(NEQns),stat=IRR)
      if ( IRR/=0 ) then
        write (*,*) '##Error : Not enough memory'
        call HECMW_ABORT(HECMW_COMM_GET_COMM())
        !stop
      endif
      !rmiv
      ndegl = NDEg*(NDEg+1)
      ndegl = ndegl/2
      ndeg2 = NDEg*NDEg
      !rmiv
      select case (NDEg)
        case (1)
          call NUFCT(XLNzr,COLno,DSLn,ZLN,DIAg,INDx,TEMp,NEQns,PARent,NCH,NSTop,Ir)
        case (2)
          call NUFCT2(XLNzr,COLno,DSLn,ZLN,DIAg,INDx,TEMp,NEQns,PARent,NCH,NSTop,Ir)
        case (3)
          call NUFCT3(XLNzr,COLno,DSLn,ZLN,DIAg,INDx,TEMp,NEQns,PARent,NCH,NSTop,Ir)
        case (6)
          call NUFCT6(XLNzr,COLno,DSLn,ZLN,DIAg,INDx,TEMp,NEQns,PARent,NCH,NSTop,Ir)
        case default
          call NUFCTX(XLNzr,COLno,DSLn,ZLN,DIAg,INDx,TEMp,NEQns,PARent,NCH,NSTop,NDEg,ndegl,Ir)
      endselect
  
      STAge = 30
      deallocate (TEMp)
      deallocate (INDx)
    end subroutine NUFCT0
  
    !======================================================================!
    !> @brief NUFCT performs cholesky factorization in row order
    !     (i) xlnzr,colno,zln,diag
    !         symbolicaly factorized
    !     (o) zln,diag,dsln
    !======================================================================!
    subroutine NUFCT(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Parent(*)
      integer(kind=kint), intent(out):: Ir
      integer(kind=kint), intent(out):: Indx(*)
      integer(kind=kint), intent(inout):: Nch(*)
      real(kind=kreal), intent(inout):: Zln(*)
      real(kind=kreal), intent(inout):: Diag(*)
      real(kind=kreal), intent(out):: Temp(*)
      real(kind=kreal), intent(inout):: Dsln(*)
      !------
      integer(kind=kint):: ic
      integer(kind=kint):: l
      integer(kind=kint):: ISEm
      real(kind=kreal):: t1
      real(kind=kreal):: t2
      real(kind=kreal):: t3
      real(kind=kreal):: t4
      real(kind=kreal):: t5
      real(kind=kreal):: tt
  
      ISEm = 1
      !
      ! phase I
      !
      call PTIME(t1)
      Diag(1) = 1.0D0/Diag(1)
      l = Parent(1)
      Nch(l) = Nch(l) - 1
      Nch(1) = -1
      do ic = 2, Nstop - 1
        call sum(ic,Xlnzr,Colno,Zln,Diag,Nch,Parent,Temp,Indx,ISEm)
      enddo
      !
      ! phase II
      !
      call PTIME(t2)
      do ic = Nstop, Neqns
        call SUM1(ic,Xlnzr,Colno,Zln,Temp,Indx)
      enddo
      !
      ! phase III
      !
      call PTIME(t3)
      call SUM2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
      !
      ! phase IV
      !
      call PTIME(t4)
      call SUM3(Neqns-Nstop+1,Dsln,Diag(Nstop),Indx,Temp)
      call PTIME(t5)
      tt = t5 - t1
      t1 = t2 - t1
      t2 = t3 - t2
      t3 = t4 - t3
      t4 = t5 - t4
      return
      Ir = 30
    end subroutine NUFCT
  
    !======================================================================!
    !> @brief NUFCT2  performs cholesky factorization in row order
    !     (i) xlnzr,colno,zln,diag
    !         symbolicaly factorized
    !     (o) zln,diag,dsln
    !======================================================================!
    subroutine NUFCT2(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Parent(*)
      integer(kind=kint), intent(out):: Ir
      integer(kind=kint), intent(out):: Indx(*)
      integer(kind=kint), intent(inout):: Nch(*)
      real(kind=kreal), intent(inout):: Zln(4,*)
      real(kind=kreal), intent(inout):: Diag(3,*)
      real(kind=kreal), intent(out):: Temp(4,*)
      real(kind=kreal), intent(inout):: Dsln(4,*)
      !------
      integer(kind=kint):: ic
      integer(kind=kint):: l
      real(kind=kreal):: t1
      real(kind=kreal):: t2
      real(kind=kreal):: t3
      real(kind=kreal):: t4
      real(kind=kreal):: t5
      real(kind=kreal):: tt
  
      !
      ! phase I
      !
      Ir = 0
      call PTIME(t1)
      if ( Nstop>1 ) call INV2(Diag(1,1),Ir)
      l = Parent(1)
      Nch(l) = Nch(l) - 1
      Nch(1) = -1
      do ic = 2, Nstop - 1
        call S2UM(ic,Xlnzr,Colno,Zln,Diag,Nch,Parent,Temp,Indx)
      enddo
      !
      ! phase II
      !
      call PTIME(t2)
      do ic = Nstop, Neqns
        call S2UM1(ic,Xlnzr,Colno,Zln,Temp,Indx)
      enddo
      !
      ! phase III
      !
      call PTIME(t3)
      call S2UM2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
      !
      ! phase IV
      !
      call PTIME(t4)
      call S2UM3(Neqns-Nstop+1,Dsln,Diag(1,Nstop),Indx,Temp)
      call PTIME(t5)
      tt = t5 - t1
      t1 = t2 - t1
      t2 = t3 - t2
      t3 = t4 - t3
      t4 = t5 - t4
    end subroutine NUFCT2
  
    !======================================================================!
    !> @brief NUFCT3 performs cholesky factorization in row order
    !     (i) xlnzr,colno,zln,diag
    !         symbolicaly factorized
    !     (o) zln,diag,dsln
    !======================================================================!
    subroutine NUFCT3(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Parent(*)
      integer(kind=kint), intent(out):: Ir
      integer(kind=kint), intent(out):: Indx(*)
      integer(kind=kint), intent(inout):: Nch(*)
      real(kind=kreal), intent(inout):: Zln(9,*)
      real(kind=kreal), intent(inout):: Diag(6,*)
      real(kind=kreal), intent(out):: Temp(*)
      real(kind=kreal), intent(inout):: Dsln(9,*)
      !------
      integer(kind=kint):: ic
      integer(kind=kint):: l
      real(kind=kreal):: t1
      real(kind=kreal):: t2
      real(kind=kreal):: t3
      real(kind=kreal):: t4
      real(kind=kreal):: t5
      real(kind=kreal):: tt
  
      !
      ! phase I
      !
      call PTIME(t1)
      if ( Nstop>1 ) call INV3(Diag(1,1),Ir)
      l = Parent(1)
      Nch(l) = Nch(l) - 1
      Nch(1) = -1
      do ic = 2, Nstop - 1
        call S3UM(ic,Xlnzr,Colno,Zln,Diag,Nch,Parent,Temp,Indx)
      enddo
      !
      ! phase II
      !
      call PTIME(t2)
      do ic = Nstop, Neqns
        call S3UM1(ic,Xlnzr,Colno,Zln,Temp,Indx)
      enddo
      !
      ! phase III
      !
      call PTIME(t3)
      call S3UM2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
      !
      ! phase IV
      !
      call PTIME(t4)
      call S3UM3(Neqns-Nstop+1,Dsln,Diag(1,Nstop),Indx,Temp)
      call PTIME(t5)
      tt = t5 - t1
      t1 = t2 - t1
      t2 = t3 - t2
      t3 = t4 - t3
      t4 = t5 - t4
    end subroutine NUFCT3
  
    !======================================================================!
    !> @brief NUFCT6 performs cholesky factorization in row order
    !     (i) xlnzr,colno,zln,diag
    !         symbolicaly factorized
    !     (o) zln,diag,dsln
    !======================================================================!
    subroutine NUFCT6(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Parent(*)
      integer(kind=kint), intent(out):: Indx(*)
      integer(kind=kint), intent(inout):: Nch(*)
      real(kind=kreal), intent(inout):: Zln(36,*)
      real(kind=kreal), intent(inout):: Diag(21,*)
      real(kind=kreal), intent(out):: Temp(*)
      real(kind=kreal), intent(inout):: Dsln(36,*)
      !------
      integer(kind=kint):: ic
      integer(kind=kint):: Ir
      integer(kind=kint):: l
      real(kind=kreal):: t1
      real(kind=kreal):: t2
      real(kind=kreal):: t3
      real(kind=kreal):: t4
      real(kind=kreal):: t5
      real(kind=kreal):: tt
  
      !
      ! phase I
      !
      call PTIME(t1)
      if ( Nstop>1 ) call INV6(Diag(1,1),Ir)
      l = Parent(1)
      Nch(l) = Nch(l) - 1
      Nch(1) = -1
      do ic = 2, Nstop - 1
        call S6UM(ic,Xlnzr,Colno,Zln,Diag,Nch,Parent,Temp,Indx)
      enddo
      !
      ! phase II
      !
      call PTIME(t2)
      do ic = Nstop, Neqns
        call S6UM1(ic,Xlnzr,Colno,Zln,Temp,Indx)
      enddo
      !
      ! phase III
      !
      call PTIME(t3)
      call S6UM2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
      !
      ! phase IV
      !
      call PTIME(t4)
      call S6UM3(Neqns-Nstop+1,Dsln,Diag(1,Nstop),Indx,Temp)
      call PTIME(t5)
      tt = t5 - t1
      t1 = t2 - t1
      t2 = t3 - t2
      t3 = t4 - t3
      t4 = t5 - t4
    end subroutine NUFCT6
  
    !======================================================================!
    !> @brief NUFCTX performs cholesky factorization in row order
    !     (i) xlnzr,colno,zln,diag
    !         symbolicaly factorized
    !     (o) zln,diag,dsln
    !======================================================================!
    subroutine NUFCTX(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ndeg,Ndegl,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ndeg
      integer(kind=kint), intent(in):: Ndegl
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Parent(*)
      integer(kind=kint), intent(out):: Indx(*)
      integer(kind=kint), intent(inout):: Nch(*)
      real(kind=kreal), intent(inout):: Zln(Ndeg*Ndeg,*)
      real(kind=kreal), intent(inout):: Diag(Ndegl,*)
      real(kind=kreal), intent(out):: Temp(Ndeg*Ndeg,*)
      real(kind=kreal), intent(inout):: Dsln(Ndeg*Ndeg,*)
      !------
      integer(kind=kint):: ic
      integer(kind=kint):: Ir
      integer(kind=kint):: l
      real(kind=kreal):: t1
      real(kind=kreal):: t2
      real(kind=kreal):: t3
      real(kind=kreal):: t4
      real(kind=kreal):: t5
      real(kind=kreal):: tt
      real(kind=kreal):: zz(100)
      real(kind=kreal):: t(100)
  
      !
      ! phase I
      !
      call PTIME(t1)
      if ( Nstop>1 ) call INVX(Diag(1,1),Ndeg,Ir)
      l = Parent(1)
      Nch(l) = Nch(l) - 1
      Nch(1) = -1
      do ic = 2, Nstop - 1
        call SXUM(ic,Xlnzr,Colno,Zln,Diag,Nch,Parent,Temp,Indx,Ndeg,Ndegl,zz,t)
      enddo
      !
      ! phase II
      !
      call PTIME(t2)
      do ic = Nstop, Neqns
        call SXUM1(ic,Xlnzr,Colno,Zln,Temp,Indx,Ndeg,t)
      enddo
      !
      ! phase III
      !
      call PTIME(t3)
      call SXUM2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx,Ndeg,Ndegl)
      !
      ! phase IV
      !
      call PTIME(t4)
      call SXUM3(Neqns-Nstop+1,Dsln,Diag(1,Nstop),Indx,Temp,Ndeg,Ndegl,t)
      call PTIME(t5)
      tt = t5 - t1
      t1 = t2 - t1
      t2 = t3 - t2
      t3 = t4 - t3
      t4 = t5 - t4
    end subroutine NUFCTX
  
    !======================================================================!
    !> @brief NUSOL0 performs forward elimination and backward substitution
    !     (i/o)
    !           r_h_s    on entry     right hand side vector
    !                    on exit      solution vector
    !           iv       communication array
    !======================================================================!
    subroutine NUSOL0(R_h_s,NEQns,NDEg,NSTop,IPErm,XLNzr,COLno,DIAg,ZLN,DSLn,STAge,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: NEQns
      integer(kind=kint), intent(in):: NDEg
      integer(kind=kint), intent(in):: NSTop
      integer(kind=kint), intent(in):: IPErm(:)
      integer(kind=kint), intent(in):: XLNzr(:)
      integer(kind=kint), intent(in):: COLno(:)
      real(kind=kreal), intent(in):: DIAg(:)
      real(kind=kreal), intent(in):: ZLN(:)
      real(kind=kreal), intent(in):: DSLn(:)
      integer(kind=kint), intent(inout):: STAge
      integer(kind=kint), intent(out):: Ir
      real(kind=kreal), intent(inout):: R_h_s(*)
      !------
      integer(kind=kint):: ndegl
      integer(kind=kint):: ierror
      real(kind=kreal), pointer :: wk(:)
  
      if ( STAge/=30 .and. STAge/=40 ) then
        Ir = 50
        return
      else
        Ir = 0
      endif
  
      allocate (wk(NDEg*NEQns),stat=IERror)
      if ( IERror/=0 ) then
        write (*,*) "##Error: not enough memory"
        call HECMW_ABORT(HECMW_COMM_GET_COMM())
      endif
      !rmiv
      ndegl = NDEg*(NDEg+1)
      ndegl = ndegl/2
  
      select case( NDEg )
        case (1)
          call NUSOL1(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop)
        case (2)
          call NUSOL2(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop)
        case (3)
          call NUSOL3(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop)
        case (6)
          call NUSOLX(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop,NDEg,ndegl)
        case default
          call NUSOLX(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop,NDEg,ndegl)
      endselect
  
      STAge = 40
      deallocate (wk)
    end subroutine NUSOL0
  
    !======================================================================!
    !> @brief NUSOL1 performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOL1(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(*)
      real(kind=kreal), intent(in):: Diag(*)
      real(kind=kreal), intent(in):: Dsln(*)
      real(kind=kreal), intent(inout):: B(*)
      real(kind=kreal), intent(out):: Wk(*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
  
      ! forward
      do i = 1, Neqns
        Wk(i) = B(Iperm(i))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) Wk(i) = Wk(i) - SPDOT2(Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          Wk(i) = Wk(i) - DDOT(Wk(Nstop),Dsln(joc),i-Nstop)
          joc = joc + i - Nstop
        endif
      enddo
      do i = 1, Neqns
        Wk(i) = Wk(i)*Diag(i)
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            Wk(j) = Wk(j) - Wk(i)*Dsln(joc)
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            Wk(j) = Wk(j) - Wk(i)*Zln(k)
          enddo
        endif
      enddo
      ! permutaion
      do i = 1, Neqns
        B(Iperm(i)) = Wk(i)
      enddo
    end subroutine NUSOL1
  
    !======================================================================!
    !> @brief NUSOL2 performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOL2(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(4,*)
      real(kind=kreal), intent(in):: Diag(3,*)
      real(kind=kreal), intent(in):: Dsln(4,*)
      real(kind=kreal), intent(inout):: B(2,*)
      real(kind=kreal), intent(out):: Wk(2,*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
  
      ! forward
      do i = 1, Neqns
        Wk(1,i) = B(1,Iperm(i))
        Wk(2,i) = B(2,Iperm(i))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) call S2PDOT(Wk(1,i),Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          call D2SDOT(Wk(1,i),Wk(1,Nstop),Dsln(1,joc),i-Nstop)
          joc = joc + i - Nstop
        endif
      enddo
      do i = 1, Neqns
        Wk(2,i) = Wk(2,i) - Wk(1,i)*Diag(2,i)
        Wk(1,i) = Wk(1,i)*Diag(1,i)
        Wk(2,i) = Wk(2,i)*Diag(3,i)
        Wk(1,i) = Wk(1,i) - Wk(2,i)*Diag(2,i)
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Dsln(1,joc) - Wk(2,i)*Dsln(2,joc)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Dsln(3,joc) - Wk(2,i)*Dsln(4,joc)
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Zln(1,k) - Wk(2,i)*Zln(2,k)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Zln(3,k) - Wk(2,i)*Zln(4,k)
          enddo
        endif
      enddo
      ! permutaion
      do i = 1, Neqns
        B(1,Iperm(i)) = Wk(1,i)
        B(2,Iperm(i)) = Wk(2,i)
      enddo
    end subroutine NUSOL2
  
    !======================================================================!
    !> @brief NUSOL3 performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOL3(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(9,*)
      real(kind=kreal), intent(in):: Diag(6,*)
      real(kind=kreal), intent(in):: Dsln(9,*)
      real(kind=kreal), intent(inout):: B(3,*)
      real(kind=kreal), intent(out):: Wk(3,*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
  
      ! forward
      do i = 1, Neqns
        Wk(1,i) = B(1,Iperm(i))
        Wk(2,i) = B(2,Iperm(i))
        Wk(3,i) = B(3,Iperm(i))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) call S3PDOT(Wk(1,i),Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          call D3SDOT(Wk(1,i),Wk(1,Nstop),Dsln(1,joc),i-Nstop)
          joc = joc + i - Nstop
        endif
      enddo
      do i = 1, Neqns
        Wk(2,i) = Wk(2,i) - Wk(1,i)*Diag(2,i)
        Wk(3,i) = Wk(3,i) - Wk(1,i)*Diag(4,i) - Wk(2,i)*Diag(5,i)
        Wk(1,i) = Wk(1,i)*Diag(1,i)
        Wk(2,i) = Wk(2,i)*Diag(3,i)
        Wk(3,i) = Wk(3,i)*Diag(6,i)
        Wk(2,i) = Wk(2,i) - Wk(3,i)*Diag(5,i)
        Wk(1,i) = Wk(1,i) - Wk(2,i)*Diag(2,i) - Wk(3,i)*Diag(4,i)
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Dsln(1,joc) - Wk(2,i)*Dsln(2,joc) - Wk(3,i)*Dsln(3,joc)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Dsln(4,joc) - Wk(2,i)*Dsln(5,joc) - Wk(3,i)*Dsln(6,joc)
            Wk(3,j) = Wk(3,j) - Wk(1,i)*Dsln(7,joc) - Wk(2,i)*Dsln(8,joc) - Wk(3,i)*Dsln(9,joc)
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Zln(1,k) - Wk(2,i)*Zln(2,k) - Wk(3,i)*Zln(3,k)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Zln(4,k) - Wk(2,i)*Zln(5,k) - Wk(3,i)*Zln(6,k)
            Wk(3,j) = Wk(3,j) - Wk(1,i)*Zln(7,k) - Wk(2,i)*Zln(8,k) - Wk(3,i)*Zln(9,k)
          enddo
        endif
      enddo
      ! permutaion
      do i = 1, Neqns
        B(1,Iperm(i)) = Wk(1,i)
        B(2,Iperm(i)) = Wk(2,i)
        B(3,Iperm(i)) = Wk(3,i)
      enddo
    end subroutine NUSOL3
  
    !======================================================================!
    !> @brief NUSOLX performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOLX(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop,Ndeg,Ndegl)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ndeg
      integer(kind=kint), intent(in):: Ndegl
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(Ndeg,Ndeg,*)
      real(kind=kreal), intent(in):: Diag(Ndegl,*)
      real(kind=kreal), intent(in):: Dsln(Ndeg,Ndeg,*)
      real(kind=kreal), intent(inout):: B(Ndeg,*)
      real(kind=kreal), intent(out):: Wk(Ndeg,*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: joc1
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
      integer(kind=kint):: l
      integer(kind=kint):: loc1
      integer(kind=kint):: locd
      integer(kind=kint):: m
      integer(kind=kint):: n
  
      ! forward
      do l = 1, Ndeg
        Wk(l,1:Neqns) = B(l,Iperm(1:Neqns))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) call SXPDOT(Ndeg,Wk(1,i),Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          joc1 = i - Nstop
          call DXSDOT(Ndeg,Wk(1,i),Wk(1,Nstop),Dsln(1,1,joc),joc1)
          joc = joc + joc1
        endif
      enddo
      do i = 1, Neqns
        locd = 0
        do m = 1, Ndeg - 1
          locd = locd + m
          loc1 = locd + m
          do n = m + 1, Ndeg
            Wk(n,i) = Wk(n,i) - Wk(m,i)*Diag(loc1,i)
            loc1 = loc1 + n
          enddo
        enddo
        locd = 0
        do m = 1, Ndeg
          locd = locd + m
          Wk(m,i) = Wk(m,i)*Diag(locd,i)
        enddo
        do n = Ndeg, 2, -1
          locd = locd - 1
          do m = n - 1, 1, -1
            Wk(m,i) = Wk(m,i) - Wk(n,i)*Diag(locd,i)
            locd = locd - 1
          enddo
        enddo
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            do m = 1, Ndeg
              do n = 1, Ndeg
                Wk(m,j) = Wk(m,j) - Wk(n,i)*Dsln(n,m,joc)
              enddo
            enddo
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            do m = 1, Ndeg
              do n = 1, Ndeg
                Wk(m,j) = Wk(m,j) - Wk(n,i)*Zln(n,m,k)
              enddo
            enddo
          enddo
        endif
      enddo
      ! permutaion
      do l = 1, Ndeg
        B(l,Iperm(1:Neqns)) = Wk(l,1:Neqns)
      enddo
    end subroutine NUSOLX
  
    !======================================================================!
    !> @brief ZPIVOT
    !======================================================================!
    subroutine ZPIVOT(Neqns,Neqnsz,Nttbr,Jcol,Irow,Zpiv,Ir)
  
negative.code:       call INVXX(Zz(1,i),Zln(1,i),Diag(1,i),Ndeg)
      enddo
    end subroutine VXPROD
  
    !======================================================================!
    !> @brief VCOPY
    !======================================================================!
    subroutine VCOPY(A,C,N)
      implicit none
      !------
      integer(kind=kint), intent(in):: N
      real(kind=kreal), intent(in):: A(N)
      real(kind=kreal), intent(out):: C(N)
      !------
      C = A
    end subroutine VCOPY
  
    !======================================================================!
    !> @brief NUSOL0 performs forward elimination and backward substitution
    !     (i/o)
    !           r_h_s    on entry     right hand side vector
    !                    on exit      solution vector
    !           iv       communication array
    !======================================================================!
    subroutine NUSOL0(R_h_s,NEQns,NDEg,NSTop,IPErm,XLNzr,COLno,DIAg,ZLN,DSLn,STAge,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: NEQns
      integer(kind=kint), intent(in):: NDEg
      integer(kind=kint), intent(in):: NSTop
      integer(kind=kint), intent(in):: IPErm(:)
      integer(kind=kint), intent(in):: XLNzr(:)
      integer(kind=kint), intent(in):: COLno(:)
      real(kind=kreal), intent(in):: DIAg(:)
      real(kind=kreal), intent(in):: ZLN(:)
      real(kind=kreal), intent(in):: DSLn(:)
      integer(kind=kint), intent(inout):: STAge
      integer(kind=kint), intent(out):: Ir
      real(kind=kreal), intent(inout):: R_h_s(*)
      !------
      integer(kind=kint):: ndegl
      integer(kind=kint):: ierror
      real(kind=kreal), pointer :: wk(:)
  
      if ( STAge/=30 .and. STAge/=40 ) then
        Ir = 50
        return
      else
        Ir = 0
      endif
  
      allocate (wk(NDEg*NEQns),stat=IERror)
      if ( IERror/=0 ) then
        write (*,*) "##Error: not enough memory"
        call HECMW_ABORT(HECMW_COMM_GET_COMM())
      endif
      !rmiv
      ndegl = NDEg*(NDEg+1)
      ndegl = ndegl/2
  
      select case( NDEg )
        case (1)
          call NUSOL1(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop)
        case (2)
          call NUSOL2(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop)
        case (3)
          call NUSOL3(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop)
        case (6)
          call NUSOLX(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop,NDEg,ndegl)
        case default
          call NUSOLX(XLNzr,COLno,DSLn,ZLN,DIAg,IPErm,R_h_s,wk,NEQns,NSTop,NDEg,ndegl)
      endselect
  
      STAge = 40
      deallocate (wk)
    end subroutine NUSOL0
  
    !======================================================================!
    !> @brief NUSOL1 performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOL1(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(*)
      real(kind=kreal), intent(in):: Diag(*)
      real(kind=kreal), intent(in):: Dsln(*)
      real(kind=kreal), intent(inout):: B(*)
      real(kind=kreal), intent(out):: Wk(*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
  
      ! forward
      do i = 1, Neqns
        Wk(i) = B(Iperm(i))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) Wk(i) = Wk(i) - SPDOT2(Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          Wk(i) = Wk(i) - DDOT(Wk(Nstop),Dsln(joc),i-Nstop)
          joc = joc + i - Nstop
        endif
      enddo
      do i = 1, Neqns
        Wk(i) = Wk(i)*Diag(i)
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            Wk(j) = Wk(j) - Wk(i)*Dsln(joc)
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            Wk(j) = Wk(j) - Wk(i)*Zln(k)
          enddo
        endif
      enddo
      ! permutaion
      do i = 1, Neqns
        B(Iperm(i)) = Wk(i)
      enddo
    end subroutine NUSOL1
  
    !======================================================================!
    !> @brief NUSOL2 performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOL2(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(4,*)
      real(kind=kreal), intent(in):: Diag(3,*)
      real(kind=kreal), intent(in):: Dsln(4,*)
      real(kind=kreal), intent(inout):: B(2,*)
      real(kind=kreal), intent(out):: Wk(2,*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
  
      ! forward
      do i = 1, Neqns
        Wk(1,i) = B(1,Iperm(i))
        Wk(2,i) = B(2,Iperm(i))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) call S2PDOT(Wk(1,i),Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          call D2SDOT(Wk(1,i),Wk(1,Nstop),Dsln(1,joc),i-Nstop)
          joc = joc + i - Nstop
        endif
      enddo
      do i = 1, Neqns
        Wk(2,i) = Wk(2,i) - Wk(1,i)*Diag(2,i)
        Wk(1,i) = Wk(1,i)*Diag(1,i)
        Wk(2,i) = Wk(2,i)*Diag(3,i)
        Wk(1,i) = Wk(1,i) - Wk(2,i)*Diag(2,i)
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Dsln(1,joc) - Wk(2,i)*Dsln(2,joc)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Dsln(3,joc) - Wk(2,i)*Dsln(4,joc)
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Zln(1,k) - Wk(2,i)*Zln(2,k)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Zln(3,k) - Wk(2,i)*Zln(4,k)
          enddo
        endif
      enddo
      ! permutaion
      do i = 1, Neqns
        B(1,Iperm(i)) = Wk(1,i)
        B(2,Iperm(i)) = Wk(2,i)
      enddo
    end subroutine NUSOL2
  
    !======================================================================!
    !> @brief NUSOL3 performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOL3(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
      implicit none
      !------
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(9,*)
      real(kind=kreal), intent(in):: Diag(6,*)
      real(kind=kreal), intent(in):: Dsln(9,*)
      real(kind=kreal), intent(inout):: B(3,*)
      real(kind=kreal), intent(out):: Wk(3,*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
  
      ! forward
      do i = 1, Neqns
        Wk(1,i) = B(1,Iperm(i))
        Wk(2,i) = B(2,Iperm(i))
        Wk(3,i) = B(3,Iperm(i))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) call S3PDOT(Wk(1,i),Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          call D3SDOT(Wk(1,i),Wk(1,Nstop),Dsln(1,joc),i-Nstop)
          joc = joc + i - Nstop
        endif
      enddo
      do i = 1, Neqns
        Wk(2,i) = Wk(2,i) - Wk(1,i)*Diag(2,i)
        Wk(3,i) = Wk(3,i) - Wk(1,i)*Diag(4,i) - Wk(2,i)*Diag(5,i)
        Wk(1,i) = Wk(1,i)*Diag(1,i)
        Wk(2,i) = Wk(2,i)*Diag(3,i)
        Wk(3,i) = Wk(3,i)*Diag(6,i)
        Wk(2,i) = Wk(2,i) - Wk(3,i)*Diag(5,i)
        Wk(1,i) = Wk(1,i) - Wk(2,i)*Diag(2,i) - Wk(3,i)*Diag(4,i)
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Dsln(1,joc) - Wk(2,i)*Dsln(2,joc) - Wk(3,i)*Dsln(3,joc)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Dsln(4,joc) - Wk(2,i)*Dsln(5,joc) - Wk(3,i)*Dsln(6,joc)
            Wk(3,j) = Wk(3,j) - Wk(1,i)*Dsln(7,joc) - Wk(2,i)*Dsln(8,joc) - Wk(3,i)*Dsln(9,joc)
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            Wk(1,j) = Wk(1,j) - Wk(1,i)*Zln(1,k) - Wk(2,i)*Zln(2,k) - Wk(3,i)*Zln(3,k)
            Wk(2,j) = Wk(2,j) - Wk(1,i)*Zln(4,k) - Wk(2,i)*Zln(5,k) - Wk(3,i)*Zln(6,k)
            Wk(3,j) = Wk(3,j) - Wk(1,i)*Zln(7,k) - Wk(2,i)*Zln(8,k) - Wk(3,i)*Zln(9,k)
          enddo
        endif
      enddo
      ! permutaion
      do i = 1, Neqns
        B(1,Iperm(i)) = Wk(1,i)
        B(2,Iperm(i)) = Wk(2,i)
        B(3,Iperm(i)) = Wk(3,i)
      enddo
    end subroutine NUSOL3
  
    !======================================================================!
    !> @brief NUSOLX performs forward elimination and backward substitution
    !======================================================================!
    subroutine NUSOLX(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop,Ndeg,Ndegl)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ndeg
      integer(kind=kint), intent(in):: Ndegl
      integer(kind=kint), intent(in):: Neqns
      integer(kind=kint), intent(in):: Nstop
      integer(kind=kint), intent(in):: Xlnzr(*)
      integer(kind=kint), intent(in):: Colno(*)
      integer(kind=kint), intent(in):: Iperm(*)
      real(kind=kreal), intent(in):: Zln(Ndeg,Ndeg,*)
      real(kind=kreal), intent(in):: Diag(Ndegl,*)
      real(kind=kreal), intent(in):: Dsln(Ndeg,Ndeg,*)
      real(kind=kreal), intent(inout):: B(Ndeg,*)
      real(kind=kreal), intent(out):: Wk(Ndeg,*)
      !------
      integer(kind=kint):: i
      integer(kind=kint):: j
      integer(kind=kint):: joc
      integer(kind=kint):: joc1
      integer(kind=kint):: k
      integer(kind=kint):: ke
      integer(kind=kint):: ks
      integer(kind=kint):: l
      integer(kind=kint):: loc1
      integer(kind=kint):: locd
      integer(kind=kint):: m
      integer(kind=kint):: n
  
      ! forward
      do l = 1, Ndeg
        Wk(l,1:Neqns) = B(l,Iperm(1:Neqns))
      enddo
      joc = 1
      do i = 1, Neqns
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) call SXPDOT(Ndeg,Wk(1,i),Wk,Zln,Colno,ks,ke)
        if ( i>Nstop ) then
          joc1 = i - Nstop
          call DXSDOT(Ndeg,Wk(1,i),Wk(1,Nstop),Dsln(1,1,joc),joc1)
          joc = joc + joc1
        endif
      enddo
      do i = 1, Neqns
        locd = 0
        do m = 1, Ndeg - 1
          locd = locd + m
          loc1 = locd + m
          do n = m + 1, Ndeg
            Wk(n,i) = Wk(n,i) - Wk(m,i)*Diag(loc1,i)
            loc1 = loc1 + n
          enddo
        enddo
        locd = 0
        do m = 1, Ndeg
          locd = locd + m
          Wk(m,i) = Wk(m,i)*Diag(locd,i)
        enddo
        do n = Ndeg, 2, -1
          locd = locd - 1
          do m = n - 1, 1, -1
            Wk(m,i) = Wk(m,i) - Wk(n,i)*Diag(locd,i)
            locd = locd - 1
          enddo
        enddo
      enddo
      ! back ward
      do i = Neqns, 1, -1
        if ( i>=Nstop ) then
          do j = i - 1, Nstop, -1
            joc = joc - 1
            do m = 1, Ndeg
              do n = 1, Ndeg
                Wk(m,j) = Wk(m,j) - Wk(n,i)*Dsln(n,m,joc)
              enddo
            enddo
          enddo
        endif
        ks = Xlnzr(i)
        ke = Xlnzr(i+1) - 1
        if ( ke>=ks ) then
          do k = ks, ke
            j = Colno(k)
            do m = 1, Ndeg
              do n = 1, Ndeg
                Wk(m,j) = Wk(m,j) - Wk(n,i)*Zln(n,m,k)
              enddo
            enddo
          enddo
        endif
      enddo
      ! permutaion
      do l = 1, Ndeg
        B(l,Iperm(1:Neqns)) = Wk(l,1:Neqns)
      enddo
    end subroutine NUSOLX
  
    !======================================================================!
    !> @brief SPDOT2 performs inner product of sparse vectors
    !======================================================================!
    real(kind=kreal) function SPDOT2(B,Zln,Colno,Ks,Ke)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ke
      integer(kind=kint), intent(in):: Ks
      integer(kind=kint), intent(in):: Colno(*)
      real(kind=kreal), intent(in):: Zln(*)
      real(kind=kreal), intent(in):: B(*)
      !------
      integer(kind=kint):: j
      integer(kind=kint):: jj
      real(kind=kreal):: s
  
      s = 0.0D0
      do jj = Ks, Ke
        j = Colno(jj)
        s = s + Zln(jj)*B(j)
      enddo
      SPDOT2 = s
    end function SPDOT2
  
    !======================================================================!
    !> @brief S2PDOT performs inner product of sparse vectors
    !======================================================================!
    subroutine S2PDOT(Bi,B,Zln,Colno,Ks,Ke)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ke
      integer(kind=kint), intent(in):: Ks
      integer(kind=kint), intent(in):: Colno(*)
      real(kind=kreal), intent(in):: Zln(4,*)
      real(kind=kreal), intent(in):: B(2,*)
      real(kind=kreal), intent(inout):: Bi(2)
      !------
      integer(kind=kint):: j
      integer(kind=kint):: jj
  
      do jj = Ks, Ke
        j = Colno(jj)
        Bi(1) = Bi(1) - Zln(1,jj)*B(1,j) - Zln(3,jj)*B(2,j)
        Bi(2) = Bi(2) - Zln(2,jj)*B(1,j) - Zln(4,jj)*B(2,j)
      enddo
    end subroutine S2PDOT
  
    !======================================================================!
    !> @brief S3PDOT performs inner product of sparse vectors
    !======================================================================!
    subroutine S3PDOT(Bi,B,Zln,Colno,Ks,Ke)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ke
      integer(kind=kint), intent(in):: Ks
      integer(kind=kint), intent(in):: Colno(*)
      real(kind=kreal), intent(in):: Zln(9,*)
      real(kind=kreal), intent(in):: B(3,*)
      real(kind=kreal), intent(inout):: Bi(3)
      !------
      integer(kind=kint):: j
      integer(kind=kint):: jj
  
      do jj = Ks, Ke
        j = Colno(jj)
        Bi(1) = Bi(1) - Zln(1,jj)*B(1,j) - Zln(4,jj)*B(2,j) - Zln(7,jj)*B(3,j)
        Bi(2) = Bi(2) - Zln(2,jj)*B(1,j) - Zln(5,jj)*B(2,j) - Zln(8,jj)*B(3,j)
        Bi(3) = Bi(3) - Zln(3,jj)*B(1,j) - Zln(6,jj)*B(2,j) - Zln(9,jj)*B(3,j)
      enddo
    end subroutine S3PDOT
  
    !======================================================================!
    !> @brief S6PDOT  performs inner product of sparse vectors
    !======================================================================!
    subroutine S6PDOT(Bi,B,Zln,Colno,Ks,Ke)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ke
      integer(kind=kint), intent(in):: Ks
      integer(kind=kint), intent(in):: Colno(*)
      real(kind=kreal), intent(in):: Zln(36,*)
      real(kind=kreal), intent(in):: B(6,*)
      real(kind=kreal), intent(inout):: Bi(6)
      !------
      integer(kind=kint):: j
      integer(kind=kint):: jj
  
      do jj = Ks, Ke
        j = Colno(jj)
        Bi(1) = Bi(1) - Zln(1,jj)*B(1,j) - Zln(7,jj)*B(2,j) - Zln(13,jj)*B(3,j)&
          - Zln(19,jj)*B(4,j) - Zln(25,jj)*B(5,j) - Zln(31,jj)*B(6,j)
        Bi(2) = Bi(2) - Zln(2,jj)*B(1,j) - Zln(8,jj)*B(2,j) - Zln(14,jj)*B(3,j)&
          - Zln(20,jj)*B(4,j) - Zln(26,jj)*B(5,j) - Zln(32,jj)*B(6,j)
        Bi(3) = Bi(3) - Zln(3,jj)*B(1,j) - Zln(9,jj)*B(2,j) - Zln(15,jj)*B(3,j)&
          - Zln(21,jj)*B(4,j) - Zln(27,jj)*B(5,j) - Zln(33,jj)*B(6,j)
        Bi(4) = Bi(4) - Zln(4,jj)*B(1,j) - Zln(10,jj)*B(2,j) - Zln(16,jj)*B(3,j)&
          - Zln(22,jj)*B(4,j) - Zln(28,jj)*B(5,j) - Zln(34,jj)*B(6,j)
        Bi(5) = Bi(5) - Zln(5,jj)*B(1,j) - Zln(11,jj)*B(2,j) - Zln(17,jj)*B(3,j)&
          - Zln(23,jj)*B(4,j) - Zln(29,jj)*B(5,j) - Zln(35,jj)*B(6,j)
        Bi(6) = Bi(6) - Zln(6,jj)*B(1,j) - Zln(12,jj)*B(2,j) - Zln(18,jj)*B(3,j)&
          - Zln(25,jj)*B(4,j) - Zln(30,jj)*B(5,j) - Zln(36,jj)*B(6,j)
      enddo
    end subroutine S6PDOT
  
    !======================================================================!
    !> @brief SXPDOT performs inner product of sparse vectors
    !======================================================================!
    subroutine SXPDOT(Ndeg,Bi,B,Zln,Colno,Ks,Ke)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ke
      integer(kind=kint), intent(in):: Ks
      integer(kind=kint), intent(in):: Ndeg
      integer(kind=kint), intent(in):: Colno(*)
      real(kind=kreal), intent(in):: Zln(Ndeg,Ndeg,*)
      real(kind=kreal), intent(in):: B(Ndeg,*)
      real(kind=kreal), intent(inout):: Bi(Ndeg)
      !------
      integer(kind=kint):: j
      integer(kind=kint):: jj
      integer(kind=kint):: m
      integer(kind=kint):: n
  
      do jj = Ks, Ke
        j = Colno(jj)
        do m = 1, Ndeg
          do n = 1, Ndeg
            Bi(n) = Bi(n) - Zln(n,m,jj)*B(m,j)
          enddo
        enddo
      enddo
    end subroutine SXPDOT
  
    !======================================================================!
    !> @brief D2SDOT performs inner product of sparse vectors
    !======================================================================!
    subroutine D2SDOT(Wi,A,B,N)
      implicit none
      !------
      integer(kind=kint), intent(in):: N
      real(kind=kreal), intent(in):: A(2,*)
      real(kind=kreal), intent(in):: B(4,*)
      real(kind=kreal), intent(inout):: Wi(2)
      !------
      integer(kind=kint):: jj
  
      do jj = 1, N
        Wi(1) = Wi(1) - A(1,jj)*B(1,jj) - A(2,jj)*B(3,jj)
        Wi(2) = Wi(2) - A(1,jj)*B(2,jj) - A(2,jj)*B(4,jj)
      enddo
    end subroutine D2SDOT
  
    !======================================================================!
    !> @brief D3SDOT performs inner product of sparse vectors
    !======================================================================!
    subroutine D3SDOT(Wi,A,B,N)
      implicit none
      !------
      integer(kind=kint), intent(in):: N
      real(kind=kreal), intent(in):: A(3,*)
      real(kind=kreal), intent(in):: B(9,*)
      real(kind=kreal), intent(inout):: Wi(3)
      !------
      integer(kind=kint):: jj
  
      do jj = 1, N
        Wi(1) = Wi(1) - A(1,jj)*B(1,jj) - A(2,jj)*B(4,jj) - A(3,jj)*B(7,jj)
        Wi(2) = Wi(2) - A(1,jj)*B(2,jj) - A(2,jj)*B(5,jj) - A(3,jj)*B(8,jj)
        Wi(3) = Wi(3) - A(1,jj)*B(3,jj) - A(2,jj)*B(6,jj) - A(3,jj)*B(9,jj)
      enddo
    end subroutine D3SDOT
  
    !======================================================================!
    !> @brief DXSDOT performs inner product of sparse vectors
    !======================================================================!
    subroutine DXSDOT(Ndeg,Wi,A,B,N)
      implicit none
      !------
      integer(kind=kint), intent(in):: Ndeg
      real(kind=kreal), intent(in):: A(Ndeg,*)
      real(kind=kreal), intent(in):: B(Ndeg,Ndeg,*)
      real(kind=kreal), intent(inout):: Wi(Ndeg)
      integer(kind=kint), intent(inout):: N
      !------
      integer(kind=kint):: jj
      integer(kind=kint):: m
  
      do jj = 1, N
        do m = 1, Ndeg
          do N = 1, Ndeg
            Wi(N) = Wi(N) - B(N,m,jj)*A(m,jj)
          enddo
        enddo
      enddo
    end subroutine DXSDOT
  
  end module HECMW_SOLVER_DIRECT
  