anchor.diff:        NU=hecMAT%indexU(NN)
         NZ=NN*(ndof2+ndof)/2+NU*ndof2 &
           +fstrMAT%numU_lagrange*ndof
  -
  -      !!$#ifdef NEED_DIAG
  -      !!$       NZ=NZ+fstrMAT%num_lagrange
  -      !!$#endif
  -
       else
         NL=hecMAT%indexL(NN)
         NU=hecMAT%indexU(NN)
         NZ=(NN+NU+NL)*ndof2 &
           +(fstrMAT%numL_lagrange+fstrMAT%numU_lagrange)*ndof
       endif
  +    ! diagonal elements must be allocated for CRS format even if they are all zero
  +    if (spMAT%type==SPARSE_MATRIX_TYPE_CSR) NZ=NZ+fstrMAT%numL_lagrange
       call sparse_matrix_init(spMAT, N_loc, NZ)
       call sparse_matrix_hec_set_conv_ext(spMAT, hecMESH, ndof)
       if(fstrMAT%num_lagrange > 0) &
  
positive.code:           ls=fstrMAT%indexU_lagrange(i-1)+1
            le=fstrMAT%indexU_lagrange(i)
            do l=ls,le
              if (spMAT%type==SPARSE_MATRIX_TYPE_COO .and. spMAT%IRN(m)/=ii) &
                stop "ERROR: sparse_matrix_contact_set_a"
              if (spMAT%JCN(m)/=j0+fstrMAT%itemU_lagrange(l)) &
                stop "ERROR: sparse_matrix_contact_set_a"
              spMAT%A(m)=fstrMAT%AU_lagrange((l-1)*ndof+idof)
  
negative.code:     ! Lower Lagrange
      if (fstrMAT%num_lagrange > 0) then
        i0=spMAT%offset+ndof*hecMAT%N
        do i=1,fstrMAT%num_lagrange
          ii=i0+i
          if (spMAT%type==SPARSE_MATRIX_TYPE_CSR) spMAT%IRN(ii-spMAT%offset)=m
          if (.not. sparse_matrix_is_sym(spMAT)) then
            ls=fstrMAT%indexL_lagrange(i-1)+1
            le=fstrMAT%indexL_lagrange(i)
            do l=ls,le
  