anchor.diff:          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)
  +        call STAIJ1(0,j,i,val,FCT%NEQns,FCT%NDEg,FCT%LEN_colno,FCT%NSTop,FCT%LEN_dsln,FCT%INVp,FCT%XLNzr,FCT%COLno,&
  +             FCT%DIAg,FCT%ZLN,FCT%DSLn,FCT%RALoc,FCT%STAge,Ir)
         enddo
       enddo
   
  
positive.code:   !                20  building up matrix
    !                30  after LU decomposition
    !                40  after solving
    !======================================================================!
    subroutine MATINI(NEQns,NTTbr,JCOl,IROw,LEN_colno,NSTop,LEN_dsln,IPErm,INVp,PARent,NCH,XLNzr,COLno,IALoc,RALoc,STAge,Ir)
      implicit none
      !------
      integer(kind=kint), intent(in):: NEQns
      integer(kind=kint), intent(in):: NTTbr
      integer(kind=kint), intent(in) :: jcol(:)
      integer(kind=kint), intent(in) :: irow(:)
      integer(kind=kint), intent(out):: LEN_colno
      integer(kind=kint), intent(out):: NSTop
      integer(kind=kint), intent(out):: LEN_dsln
      integer(kind=kint), pointer :: IPErm(:)
      integer(kind=kint), pointer :: INVp(:)
      integer(kind=kint), pointer :: PARent(:)
      integer(kind=kint), pointer :: NCH(:)
      integer(kind=kint), pointer :: XLNzr(:)
      integer(kind=kint), pointer :: COLno(:)
      integer(kind=kint), intent(out):: IALoc
      integer(kind=kint), intent(out):: RALoc
      integer(kind=kint), intent(out):: STAge
      integer(kind=kint), intent(out):: Ir
      !------
      !*Work arrays
  
negative.code:   !           r_h_s    on entry     right hand side vector
    !                    on exit      solution vector
    !           iv       communication array
    !  updates STAge
    !======================================================================!
    subroutine NUSOL0(R_h_s,FCT,Ir)
      implicit none
      !------
      real(kind=kreal), intent(inout):: R_h_s(:)
      type (cholesky_factor), intent(inout):: FCT
      integer(kind=kint), intent(out):: Ir
      !------
      integer(kind=kint):: ndegl
      integer(kind=kint):: ierror
      real(kind=kreal), pointer :: wk(:)
  
      if ( FCT%STAge/=30 .and. FCT%STAge/=40 ) then
        Ir = 50
        return
      else
        Ir = 0
      endif
  
      allocate (wk(FCT%NDEg*FCT%NEQns),stat=IERror)
      if ( IERror/=0 ) then
        write (*,*) "##Error: not enough memory"
        call HECMW_ABORT(HECMW_COMM_GET_COMM())
      endif
      !rmiv
      ndegl = FCT%NDEg*(FCT%NDEg+1)
      ndegl = ndegl/2
  
      select case( FCT%NDEg )
        case (1)
          call NUSOL1(FCT%XLNzr,FCT%COLno,FCT%DSLn,FCT%ZLN,FCT%DIAg,FCT%IPErm,R_h_s,wk,FCT%NEQns,FCT%NSTop)
        case (2)
          call NUSOL2(FCT%XLNzr,FCT%COLno,FCT%DSLn,FCT%ZLN,FCT%DIAg,FCT%IPErm,R_h_s,wk,FCT%NEQns,FCT%NSTop)
        case (3)
          call NUSOL3(FCT%XLNzr,FCT%COLno,FCT%DSLn,FCT%ZLN,FCT%DIAg,FCT%IPErm,R_h_s,wk,FCT%NEQns,FCT%NSTop)
        case (6)
          call NUSOLX(FCT%XLNzr,FCT%COLno,FCT%DSLn,FCT%ZLN,FCT%DIAg,FCT%IPErm,R_h_s,wk,FCT%NEQns,FCT%NSTop,FCT%NDEg,ndegl)
        case default
          call NUSOLX(FCT%XLNzr,FCT%COLno,FCT%DSLn,FCT%ZLN,FCT%DIAg,FCT%IPErm,R_h_s,wk,FCT%NEQns,FCT%NSTop,FCT%NDEg,ndegl)
      endselect
  
      FCT%STAge = 40
      deallocate (wk)
    end subroutine NUSOL0
  
  