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:     ! Lower Lagrange
      if (fstrMAT%num_lagrange > 0) then
        i0=spMAT%offset+ndof*hecMAT%N
        if (.not. sparse_matrix_is_sym(spMAT)) then
          do i=1,fstrMAT%num_lagrange
            ii=i0+i
            if (spMAT%type == SPARSE_MATRIX_TYPE_CSR .and. spMAT%IRN(ii-spMAT%offset)/=m) &
              stop "ERROR: sparse_matrix_contact_set_a"
            ls=fstrMAT%indexL_lagrange(i-1)+1
            le=fstrMAT%indexL_lagrange(i)
            do l=ls,le
              j=fstrMAT%itemL_lagrange(l)
              j0=spMAT%offset+ndof*(j-1)
              do jdof=1,ndof
                if (spMAT%type==SPARSE_MATRIX_TYPE_COO .and. spMAT%IRN(m)/=ii) &
                  stop "ERROR: sparse_matrix_contact_set_a"
                if (spMAT%JCN(m)/=j0+jdof) &
                  stop "ERROR: sparse_matrix_contact_set_a"
                spMAT%A(m)=fstrMAT%AL_lagrange((l-1)*ndof+jdof)
                m=m+1
              enddo
            enddo
          enddo
          !!$#ifdef NEED_DIAG
          !!$      else
          !!$        do i=1,fstrMAT%num_lagrange
          !!$          ii=i0+i
          !!$          if (spMAT%type==SPARSE_MATRIX_TYPE_CSR .and. spMAT%IRN(ii-spMAT%OFFSET)/=m) &
            !!$               stop "ERROR: sparse_matrix_contact_set_a"
          !!$          if (spMAT%type==SPARSE_MATRIX_TYPE_COO .and. spMAT%IRN(m)/=ii) &
            !!$               stop "ERROR: sparse_matrix_contact_set_a"
          !!$          if (spMAT%JCN(m)/=ii) &
            !!$               stop "ERROR: sparse_matrix_contact_set_a"
          !!$          spMAT%A(m)=0.0d0
          !!$          m=m+1
          !!$        enddo
          !!$#endif
        endif
      endif
      if (spMAT%type == SPARSE_MATRIX_TYPE_CSR .and. spMAT%IRN(ii+1-spMAT%offset)/=m) &
        stop "ERROR: sparse_matrix_contact_set_a"
      if (m-1 /= spMAT%NZ) stop "ERROR: sparse_matrix_contact_set_a"
  
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
  