anchor.diff:  
         T(5) = T(5) + A(3,jj)*B(2,jj) + A(6,jj)*B(5,jj) + A(9,jj)*B(8,jj)
         T(6) = T(6) + A(3,jj)*B(3,jj) + A(6,jj)*B(6,jj) + A(9,jj)*B(9,jj)
  -    ENDDO
  -  END SUBROUTINE D3DOTL
  +    enddo
  +  end subroutine D3DOTL
   
     !======================================================================!
     !> @brief D3SDOT performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE D3SDOT(Wi,A,B,N)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: N
  -    REAL(KIND=8), INTENT(IN):: A(3,*)
  -    REAL(KIND=8), INTENT(IN):: B(9,*)
  -    REAL(KIND=8), INTENT(OUT):: Wi(3)
  -    !------
  -    INTEGER:: jj
  -
  -    DO jj = 1, N
  +  subroutine D3SDOT(Wi,A,B,N)
  +    implicit none
  +    !------
  +    integer, intent(IN):: N
  +    real(KIND=8), intent(IN):: A(3,*)
  +    real(KIND=8), intent(IN):: B(9,*)
  +    real(KIND=8), intent(OUT):: Wi(3)
  +    !------
  +    integer:: 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
  +    enddo
  +  end subroutine D3SDOT
   
     !======================================================================!
     !> @brief DDOT performs inner product of sparse vectors
     !======================================================================!
  -  REAL(KIND=8) FUNCTION DDOT(A,B,N)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: N
  -    REAL(KIND=8), INTENT(IN):: A(N)
  -    REAL(KIND=8), INTENT(IN):: B(N)
  -    !------
  -    INTEGER:: i
  -    REAL(KIND=8):: s
  +  real(KIND=8) function DDOT(A,B,N)
  +    implicit none
  +    !------
  +    integer, intent(IN):: N
  +    real(KIND=8), intent(IN):: A(N)
  +    real(KIND=8), intent(IN):: B(N)
  +    !------
  +    integer:: i
  +    real(KIND=8):: s
   
       s = 0.0D0
  -    DO i = 1, N
  +    do i = 1, N
         s = s + A(i)*B(i)
  -    ENDDO
  +    enddo
       DDOT = s
  -  END FUNCTION DDOT
  +  end function DDOT
   
     !======================================================================!
     !> @brief DXDOT performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE DXDOT(Ndeg,T,A,B,L)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: L
  -    INTEGER, INTENT(IN):: Ndeg
  -    REAL(KIND=8), INTENT(IN):: A(Ndeg,Ndeg,*)
  -    REAL(KIND=8), INTENT(IN):: B(Ndeg,Ndeg,*)
  -    REAL(KIND=8), INTENT(OUT):: T(Ndeg,Ndeg)
  -    !------
  -    INTEGER:: jj
  -    INTEGER:: k
  -    INTEGER:: m
  -    INTEGER:: n
  -
  -    DO n = 1, Ndeg
  -      DO m = 1, Ndeg
  +  subroutine DXDOT(Ndeg,T,A,B,L)
  +    implicit none
  +    !------
  +    integer, intent(IN):: L
  +    integer, intent(IN):: Ndeg
  +    real(KIND=8), intent(IN):: A(Ndeg,Ndeg,*)
  +    real(KIND=8), intent(IN):: B(Ndeg,Ndeg,*)
  +    real(KIND=8), intent(OUT):: T(Ndeg,Ndeg)
  +    !------
  +    integer:: jj
  +    integer:: k
  +    integer:: m
  +    integer:: n
  +
  +    do n = 1, Ndeg
  +      do m = 1, Ndeg
           T(n,m) = 0.0D0
  -        DO k = 1, Ndeg
  -          DO jj = 1, L
  +        do k = 1, Ndeg
  +          do jj = 1, L
               T(n,m) = T(n,m) + A(n,k,jj)*B(m,k,jj)
  -          ENDDO
  -        ENDDO
  -      ENDDO
  -    ENDDO
  -  END SUBROUTINE DXDOT
  +          enddo
  +        enddo
  +      enddo
  +    enddo
  +  end subroutine DXDOT
   
     !======================================================================!
     !> @brief DXDOTL performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE DXDOTL(Ndeg,T,A,B,L)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: L
  -    INTEGER, INTENT(IN):: Ndeg
  -    REAL(KIND=8), INTENT(IN):: A(Ndeg,Ndeg,*)
  -    REAL(KIND=8), INTENT(IN):: B(Ndeg,Ndeg,*)
  -    REAL(KIND=8), INTENT(OUT):: T(Ndeg,Ndeg)
  -    !------
  -    INTEGER:: jj
  -    INTEGER:: k
  -    INTEGER:: m
  -    INTEGER:: n
  -
  -    DO n = 1, Ndeg
  -      DO m = 1, n
  +  subroutine DXDOTL(Ndeg,T,A,B,L)
  +    implicit none
  +    !------
  +    integer, intent(IN):: L
  +    integer, intent(IN):: Ndeg
  +    real(KIND=8), intent(IN):: A(Ndeg,Ndeg,*)
  +    real(KIND=8), intent(IN):: B(Ndeg,Ndeg,*)
  +    real(KIND=8), intent(OUT):: T(Ndeg,Ndeg)
  +    !------
  +    integer:: jj
  +    integer:: k
  +    integer:: m
  +    integer:: n
  +
  +    do n = 1, Ndeg
  +      do m = 1, n
           T(n,m) = 0.0D0
  -        DO k = 1, Ndeg
  -          DO jj = 1, L
  +        do k = 1, Ndeg
  +          do jj = 1, L
               T(n,m) = T(n,m) + A(n,k,jj)*B(m,k,jj)
  -          ENDDO
  -        ENDDO
  -      ENDDO
  -    ENDDO
  -  END SUBROUTINE DXDOTL
  +          enddo
  +        enddo
  +      enddo
  +    enddo
  +  end subroutine DXDOTL
   
     !======================================================================!
     !> @brief DXSDOT performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE DXSDOT(Ndeg,Wi,A,B,N)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: Ndeg
  -    REAL(KIND=8), INTENT(IN):: A(Ndeg,*)
  -    REAL(KIND=8), INTENT(IN):: B(Ndeg,Ndeg,*)
  -    REAL(KIND=8), INTENT(OUT):: Wi(Ndeg)
  -    INTEGER, INTENT(INOUT):: N
  -    !------
  -    INTEGER:: jj
  -    INTEGER:: m
  -
  -    DO jj = 1, N
  -      DO m = 1, Ndeg
  -        DO N = 1, Ndeg
  +  subroutine DXSDOT(Ndeg,Wi,A,B,N)
  +    implicit none
  +    !------
  +    integer, intent(IN):: Ndeg
  +    real(KIND=8), intent(IN):: A(Ndeg,*)
  +    real(KIND=8), intent(IN):: B(Ndeg,Ndeg,*)
  +    real(KIND=8), intent(OUT):: Wi(Ndeg)
  +    integer, intent(INOUT):: N
  +    !------
  +    integer:: jj
  +    integer:: 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
  +        enddo
  +      enddo
  +    enddo
  +  end subroutine DXSDOT
   
     !======================================================================!
     !> @brief S3PDOT performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE S3PDOT(Bi,B,Zln,Colno,Ks,Ke)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: Ke
  -    INTEGER, INTENT(IN):: Ks
  -    INTEGER, INTENT(IN):: Colno(*)
  -    REAL(KIND=8), INTENT(IN):: Zln(9,*)
  -    REAL(KIND=8), INTENT(IN):: B(3,*)
  -    REAL(KIND=8), INTENT(OUT):: Bi(3)
  -    !------
  -    INTEGER:: j
  -    INTEGER:: jj
  -
  -    DO jj = Ks, Ke
  +  subroutine S3PDOT(Bi,B,Zln,Colno,Ks,Ke)
  +    implicit none
  +    !------
  +    integer, intent(IN):: Ke
  +    integer, intent(IN):: Ks
  +    integer, intent(IN):: Colno(*)
  +    real(KIND=8), intent(IN):: Zln(9,*)
  +    real(KIND=8), intent(IN):: B(3,*)
  +    real(KIND=8), intent(OUT):: Bi(3)
  +    !------
  +    integer:: j
  +    integer:: 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
  +    enddo
  +  end subroutine S3PDOT
   
     !======================================================================!
     !> @brief S2PDOT performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE S2PDOT(Bi,B,Zln,Colno,Ks,Ke)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: Ke
  -    INTEGER, INTENT(IN):: Ks
  -    INTEGER, INTENT(IN):: Colno(*)
  -    REAL(KIND=8), INTENT(IN):: Zln(4,*)
  -    REAL(KIND=8), INTENT(IN):: B(2,*)
  -    REAL(KIND=8), INTENT(OUT):: Bi(2)
  -    !------
  -    INTEGER:: j
  -    INTEGER:: jj
  -
  -    DO jj = Ks, Ke
  +  subroutine S2PDOT(Bi,B,Zln,Colno,Ks,Ke)
  +    implicit none
  +    !------
  +    integer, intent(IN):: Ke
  +    integer, intent(IN):: Ks
  +    integer, intent(IN):: Colno(*)
  +    real(KIND=8), intent(IN):: Zln(4,*)
  +    real(KIND=8), intent(IN):: B(2,*)
  +    real(KIND=8), intent(OUT):: Bi(2)
  +    !------
  +    integer:: j
  +    integer:: 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
  +    enddo
  +  end subroutine S2PDOT
   
     !======================================================================!
     !> @brief S6PDOT  performs inner product of sparse vectors
     !======================================================================!
  -  SUBROUTINE S6PDOT(Bi,B,Zln,Colno,Ks,Ke)
  -    IMPLICIT NONE
  -    !------
  -    INTEGER, INTENT(IN):: Ke
  -    INTEGER, INTENT(IN):: Ks
  -    INTEGER, INTENT(IN):: Colno(*)
  -    REAL(KIND=8), INTENT(IN):: Zln(36,*)
  -    REAL(KIND=8), INTENT(IN):: B(6,*)
  -    REAL(KIND=8), INTENT(OUT):: Bi(6)
  -    !------
  -    INTEGER:: j
  -    INTEGER:: jj
  -
  -    DO jj = Ks, Ke
  +  subroutine S6PDOT(Bi,B,Zln,Colno,Ks,Ke)
  +    implicit none
  +    !------
  +    integer, intent(IN):: Ke
  +    integer, intent(IN):: Ks
  +    integer, intent(IN):: Colno(*)
  +    real(KIND=8), intent(IN):: Zln(36,*)
  +    real(KIND=8), intent(IN):: B(6,*)
  +    real(KIND=8), intent(OUT):: Bi(6)
  +    !------
  +    integer:: j
  +    integer:: 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)
  
positive.code:           zz(36) = zz(36) - Temp(6,j)*Zln(6,jj) - Temp(12,j)*Zln(12,jj)&
                - Temp(18,j)*Zln(18,jj) - Temp(24,j)*Zln(24,jj)&
                - Temp(30,j)*Zln(30,jj) - Temp(36,j)*Zln(36,jj)
          ENDIF
        ENDDO
        CALL INV66(Zln(1,k),zz,Diag(1,jc))
        Temp(1:36,jc) = zz(1:36)
  
        t(1) = t(1) + zz(1)*Zln(1,k) + zz(7)*Zln(7,k) + zz(13)*Zln(13,k)&
  
negative.code: 
      left = iv1 + 5*neqns1
  
      allocate (IPErm(NEQns),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, iperm: SUB. matini"
      allocate (INVp(NEQns),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, invp: SUB. matini"
      allocate (DEG(NEQns+1),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, deg: SUB. matini"
      allocate (MARker(NEQns+1),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, marker: SUB. matini"
      allocate (RCHset(0:NEQns+1),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, rchset: SUB. matini"
      allocate (NBRhd(NEQns+1),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, nbrhd: SUB. matini"
      allocate (QSIze(NEQns+1),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, qsize: SUB. matini"
      allocate (QLInk(NEQns+1),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, qlink: SUB. matini"
      allocate (ADJncy(2*NTTbr),STAT=IERror)
      if ( IERror/=0 ) stop "ALLOCATION ERROR, adjncy: SUB. matini"
      call GENQMD(neqnsz,IA,JA,IPErm,INVp,DEG,MARker,RCHset,NBRhd,QSIze,QLInk,NOFsub,ADJncy)
      do
        !   build up the parent vector parent vector will be saved in
        !   work2 for a while
        call GENPAQ(IA,JA,INVp,IPErm,MARker,NEQns,RCHset)
        !
        !   build up the binary tree
        !
        lbtree = lwk3 - 2*NEQns
        last = lbtree
        !rmem
        allocate (BTRee(2*(NEQns+1)),STAT=IERror)
        if ( IERror/=0 ) stop "ALLOCATION ERROR, btree: SUB. matini"
        call GENBTQ(INVp,MARker,BTRee,ZPIv,izz,NEQns)
        !
        !   rotate the binary tree to avoid a zero pivot
        !
        if ( izz==0 ) then
          !
          !   post ordering
          !
          lpordr = last - neqns1
          last = lpordr
          !rmem
          allocate (PARent(NEQns),STAT=IERror)
          if ( IERror/=0 ) stop "ALLOCATION ERROR, parent: SUB. matini.f"
          allocate (NCH(NEQns+1),STAT=IERror)
          if ( IERror/=0 ) stop "ALLOCATION ERROR, nch: SUB. matini.f"
          allocate (PORdr(NEQns+1),STAT=IERror)
          if ( IERror/=0 ) stop "ALLOCATION ERROR, pordr: SUB. matini.f"
          call POSORD(PARent,BTRee,INVp,IPErm,PORdr,NCH,NEQns,DEG,MARker,RCHset)
          !
          !   generate skelton graph
          !
  