anchor.diff:    !         symbolicaly factorized
     !     (o) zln,diag,dsln
     !======================================================================!
  -  SUBROUTINE NUFCT2(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: Neqns
  -    INTEGER, INTENT(IN):: Nstop
  -    INTEGER, INTENT(IN):: Xlnzr(*)
  -    INTEGER, INTENT(IN):: Colno(*)
  -    INTEGER, INTENT(IN):: Parent(*)
  -    INTEGER, INTENT(OUT):: Ir
  -    INTEGER, INTENT(OUT):: Indx(*)
  -    INTEGER, INTENT(OUT):: Nch(*)
  -    REAL(KIND=8), INTENT(OUT):: Zln(4,*)
  -    REAL(KIND=8), INTENT(OUT):: Diag(3,*)
  -    REAL(KIND=8), INTENT(OUT):: Temp(4,*)
  -    REAL(KIND=8), INTENT(OUT):: Dsln(4,*)
  -    !------
  -    INTEGER:: ic
  -    INTEGER:: l
  -    REAL(KIND=8):: t1
  -    REAL(KIND=8):: t2
  -    REAL(KIND=8):: t3
  -    REAL(KIND=8):: t4
  -    REAL(KIND=8):: t5
  -    REAL(KIND=8):: tt
  +  subroutine NUFCT2(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
  +    implicit none
  +    !------
  +    integer, intent(IN):: Neqns
  +    integer, intent(IN):: Nstop
  +    integer, intent(IN):: Xlnzr(*)
  +    integer, intent(IN):: Colno(*)
  +    integer, intent(IN):: Parent(*)
  +    integer, intent(OUT):: Ir
  +    integer, intent(OUT):: Indx(*)
  +    integer, intent(OUT):: Nch(*)
  +    real(KIND=8), intent(OUT):: Zln(4,*)
  +    real(KIND=8), intent(OUT):: Diag(3,*)
  +    real(KIND=8), intent(OUT):: Temp(4,*)
  +    real(KIND=8), intent(OUT):: Dsln(4,*)
  +    !------
  +    integer:: ic
  +    integer:: l
  +    real(KIND=8):: t1
  +    real(KIND=8):: t2
  +    real(KIND=8):: t3
  +    real(KIND=8):: t4
  +    real(KIND=8):: t5
  +    real(KIND=8):: tt
   
       !
       ! phase I
       !
       Ir = 0
  -    CALL PTIME(t1)
  -    IF ( Nstop>1 ) CALL INV2(Diag(1,1),Ir)
  +    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
  +    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
  +    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)
  +    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)
  +    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
  +  end subroutine NUFCT2
   
     !======================================================================!
     !> @brief NUFCT3 performs cholesky factorization in row order
  
positive.code:         ladp = lxleaf - neqns1
          last = ladp
          !rmem
          ALLOCATE (ADJncp(NEQns+1),STAT=IERror)
          IF ( IERror/=0 ) STOP "ALLOCATION ERROR, adjncp: SUB. matini.f"
          ALLOCATE (XLEaf(NEQns+1),STAT=IERror)
          IF ( IERror/=0 ) STOP "ALLOCATION ERROR, xleaf: SUB. matini.f"
          ALLOCATE (LEAf(NTTbr),STAT=IERror)
          IF ( IERror/=0 ) STOP "ALLOCATION ERROR, leaf: SUB. matini.f"
          CALL GNLEAF(IA,JA,INVp,IPErm,NCH,ADJncp,XLEaf,LEAf,NEQns,lnleaf)
          CALL FORPAR(NEQns,PARent,NCH,NSTop)
          !*Deallocation of work arrays
          DEALLOCATE (IA)
          DEALLOCATE (JA)
          DEALLOCATE (DEG)
          DEALLOCATE (MARker)
          DEALLOCATE (RCHset)
          DEALLOCATE (NBRhd)
          DEALLOCATE (QSIze)
          DEALLOCATE (QLInk)
          DEALLOCATE (ADJncy)
          DEALLOCATE (ZPIv)
          !*Nullify pointers
          NULLIFY (IA)
          NULLIFY (JA)
          NULLIFY (DEG)
          NULLIFY (MARker)
          NULLIFY (RCHset)
          NULLIFY (NBRhd)
          NULLIFY (QSIze)
          NULLIFY (QLInk)
          NULLIFY (ADJncy)
          NULLIFY (ZPIv)
          !
          !   build up xlnzr,colno  (this is the symbolic fct.)
          !
          maxl = lxleaf - (left+neqns1)
          ALLOCATE (XLNzr(NEQns+1),STAT=IERror)
          IF ( IERror/=0 ) STOP "ALLOCATION ERROR, xlnzr: SUB. matini.f"
          CALL PRE_GNCLNO(PARent,XLEaf,LEAf,XLNzr,NEQns,NSTop,lncol,ir1)
          ALLOCATE (COLno(lncol),STAT=IERror)
          IF ( IERror/=0 ) STOP "ALLOCATION ERROR, colno: SUB. matini.f"
          CALL GNCLNO(PARent,XLEaf,LEAf,XLNzr,COLno,NEQns,NSTop,lncol,ir1)
          !*Deallocate work arrays
          DEALLOCATE (PORdr)
          DEALLOCATE (ADJncp)
          DEALLOCATE (XLEaf)
          DEALLOCATE (LEAf)
          DEALLOCATE (BTRee)
          !*Nullify pointers
          NULLIFY (PORdr)
          NULLIFY (ADJncp)
          NULLIFY (XLEaf)
          NULLIFY (LEAf)
          NULLIFY (BTRee)
          !rmem
          left = (left+neqns1) + lncol
          !rmiv
  
negative.code:     Zln(9) = Zln(9)*Diag(6)
      Zln(6) = Zln(6) - Zln(9)*Diag(5)
      Zln(3) = Zln(3) - Zln(6)*Diag(2) - Zln(9)*Diag(4)
    end subroutine INV33
  
    !======================================================================!
    !> @brief INV6
    !======================================================================!
    subroutine INV6(Dsln,Ir)
      implicit none
      !------
      integer, intent(OUT):: Ir
      real(KIND=8), intent(OUT):: Dsln(21)
      !------
      real(KIND=8):: t(5)
  
      Ir = 0
      Dsln(1) = 1.0D0/Dsln(1)
  