anchor.diff:        enddo
       enddo outer
     end subroutine place_one_on_diag
  -
  -  subroutine blk_trimatmul_add(ndof, A, B, C, ABC)
  -    implicit none
  -    integer, intent(in) :: ndof
  -    real(kind=kreal), intent(in) :: A(:), B(:), C(:)
  -    real(kind=kreal), intent(inout) :: ABC(:)
  -    real(kind=kreal), allocatable :: AB(:)
  -    integer :: ndof2, i, j, k, i0, j0, ij, ik, jk
  -
  -    ndof2=ndof*ndof
  -    allocate(AB(ndof2))
  -    AB=0.d0
  -
  -    do i=1,ndof
  -      i0=(i-1)*ndof
  -      do j=1,ndof
  -        ij=i0+j
  -        j0=(j-1)*ndof
  -        do k=1,ndof
  -          ik=i0+k
  -          jk=j0+k
  -          AB(ik)=AB(ik)+A(ij)*B(jk)
  -        enddo
  -      enddo
  -    enddo
  -
  -    do i=1,ndof
  -      i0=(i-1)*ndof
  -      do j=1,ndof
  -        ij=i0+j
  -        j0=(j-1)*ndof
  -        do k=1,ndof
  -          ik=i0+k
  -          jk=j0+k
  -          ABC(ik)=ABC(ik)+AB(ij)*C(jk)
  -        enddo
  -      enddo
  -    enddo
  -
  -    deallocate(AB)
  -  end subroutine blk_trimatmul_add
   
     subroutine replace_hecmat(hecMAT, BTtKT)
       implicit none
  
positive.code:       if (hecMAT%itemU(i) <= 0) stop 'ERROR: negative itemU'
        if (hecMAT%itemU(i) > n) stop 'ERROR: too big itemU'
      enddo
  
  !!$    ! check diagonal
  !!$    do i=1,n
  !!$      write(*,*) hecMAT%D((i-1)*ndof2+1:i*ndof2)
  !!$      do idof=1,ndof
  !!$        idx=(i-1)*ndof2+(idof-1)*ndof+idof
  !!$        if (hecMAT%D(idx) == 0) then
  !!$          hecMAT%D(idx)=1.d0
  !!$        endif
  !!$      enddo
  !!$    enddo
    end subroutine replace_hecmat
  
    subroutine make_new_hecmat(hecMAT, BTtKT, hecTKT)
  
negative.code: 
    end subroutine trimatmul_TtKT
  
    subroutine blk_trimatmul_add(ndof, A, B, C, ABC)
      implicit none
      integer, intent(in) :: ndof
      real(kind=kreal), intent(in) :: A(:), B(:), C(:)
      real(kind=kreal), intent(inout) :: ABC(:)
      real(kind=kreal), allocatable :: AB(:)
      integer :: ndof2, i, j, k, i0, j0, ij, ik, jk
  
      ndof2=ndof*ndof
      allocate(AB(ndof2))
      AB=0.d0
  
      do i=1,ndof
        i0=(i-1)*ndof
        do j=1,ndof
          ij=i0+j
          j0=(j-1)*ndof
          do k=1,ndof
            ik=i0+k
            jk=j0+k
            AB(ik)=AB(ik)+A(ij)*B(jk)
          enddo
        enddo
      enddo
  
      do i=1,ndof
        i0=(i-1)*ndof
        do j=1,ndof
          ij=i0+j
          j0=(j-1)*ndof
          do k=1,ndof
            ik=i0+k
            jk=j0+k
            ABC(ik)=ABC(ik)+AB(ij)*C(jk)
          enddo
        enddo
      enddo
  
      deallocate(AB)
    end subroutine blk_trimatmul_add
  
    subroutine place_one_on_diag(BTtKT, iwS, num_lagrange)
      implicit none
      type (T_matrix), intent(inout) :: BTtKT
  