anchor.diff:    use m_hecmw_comm_f
     use hecmw_matrix_contact
     use hecmw_matrix_misc
  -  !$ use omp_lib
   
     private
   
  
positive.code:     call hecmw_sainv_make_u_nn(hecMAT)
  
    end subroutine hecmw_precond_nn_SAINV_setup
  
    subroutine hecmw_sainv_lu_nn()
      implicit none
      integer(kind=kint) :: i,j,js,je,in
      real(kind=kreal) :: X1, X2, X3, X(NDOF)
  
      do i=1, N
        write(*,*)NDOF
        SAINVD(9*i-5) = SAINVD(9*i-5)*SAINVD(9*i-4)
        SAINVD(9*i-2) = SAINVD(9*i-2)*SAINVD(9*i  )
        SAINVD(9*i-1) = SAINVD(9*i-1)*SAINVD(9*i  )
      enddo
  
      do i=1, N
        js = inumFI1L(i-1)+1
        je = inumFI1L(i)
        do j= js,je
          in= FI1L(j)
          X1= SAINVD(9*i-8)
          X2= SAINVD(9*i-4)
          X3= SAINVD(9*i  )
          SAINVL(9*j-8) = SAINVL(9*j-8)*X1
          SAINVL(9*j-7) = SAINVL(9*j-7)*X1
          SAINVL(9*j-6) = SAINVL(9*j-6)*X1
          SAINVL(9*j-5) = SAINVL(9*j-5)*X2
          SAINVL(9*j-4) = SAINVL(9*j-4)*X2
          SAINVL(9*j-3) = SAINVL(9*j-3)*X2
          SAINVL(9*j-2) = SAINVL(9*j-2)*X3
          SAINVL(9*j-1) = SAINVL(9*j-1)*X3
          SAINVL(9*j  ) = SAINVL(9*j  )*X3
        enddo
      enddo
  
    end subroutine hecmw_sainv_lu_nn
  
    subroutine hecmw_precond_nn_SAINV_apply(R, ZP)
      implicit none
      real(kind=kreal), intent(inout)  :: ZP(:)
      real(kind=kreal), intent(in)  :: R(:)
      integer(kind=kint) :: in, i, j, isL, ieL, isU, ieU, k, iold, rcm
      real(kind=kreal) :: SW1, SW2, SW3, X1, X2, X3
  
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP&PRIVATE(i,X1,X2,X3,SW1,SW2,SW3,j,in,isL,ieL,isU,ieU) &
  !$OMP&SHARED(N,SAINVD,SAINVL,SAINVU,inumFI1U,FI1U,inumFI1L,FI1L,R,T,ZP)
  !$OMP DO
        !C-- FORWARD
        do i= 1, N
          SW1= 0.0d0
          SW2= 0.0d0
          SW3= 0.0d0
  
          isL= inumFI1L(i-1)+1
          ieL= inumFI1L(i)
          do j= isL, ieL
            in= FI1L(j)
            X1= R(3*in-2)
            X2= R(3*in-1)
            X3= R(3*in  )
            SW1= SW1 + SAINVL(9*j-8)*X1 + SAINVL(9*j-7)*X2 + SAINVL(9*j-6)*X3
            SW2= SW2 + SAINVL(9*j-5)*X1 + SAINVL(9*j-4)*X2 + SAINVL(9*j-3)*X3
            SW3= SW3 + SAINVL(9*j-2)*X1 + SAINVL(9*j-1)*X2 + SAINVL(9*j  )*X3
          enddo
  
          X1= R(3*i-2)
          X2= R(3*i-1)
          X3= R(3*i  )
  
          T(3*i-2)= (X1 + SW1)*SAINVD(9*i-8)
          T(3*i-1)= (X2 + SAINVD(9*i-7)*X1 + SW2)*SAINVD(9*i-4)
          T(3*i  )= (X3 + SAINVD(9*i-6)*X1 + SAINVD(9*i-3)*X2 + SW3)*SAINVD(9*i  )
        enddo
  !$OMP END DO
  !$OMP DO
        !C-- BACKWARD
        do i= 1, N
          SW1= 0.0d0
          SW2= 0.0d0
          SW3= 0.0d0
  
          isU= inumFI1U(i-1) + 1
          ieU= inumFI1U(i)
          do j= isU, ieU
            in= FI1U(j)
            X1= T(3*in-2)
            X2= T(3*in-1)
            X3= T(3*in  )
            SW1= SW1 + SAINVU(9*j-8)*X1 + SAINVU(9*j-7)*X2 + SAINVU(9*j-6)*X3
            SW2= SW2 + SAINVU(9*j-5)*X1 + SAINVU(9*j-4)*X2 + SAINVU(9*j-3)*X3
            SW3= SW3 + SAINVU(9*j-2)*X1 + SAINVU(9*j-1)*X2 + SAINVU(9*j  )*X3
          enddo
  
          X1= T(3*i-2)
          X2= T(3*i-1)
          X3= T(3*i  )
  
          ZP(3*i-2)= X1 + SW1 + SAINVD(9*i-7)*X2 + SAINVD(9*i-6)*X3
          ZP(3*i-1)= X2 + SW2 + SAINVD(9*i-3)*X3
          ZP(3*i  )= X3 + SW3
        enddo
  !$OMP END DO
  !$OMP END PARALLEL
  
    end subroutine hecmw_precond_nn_SAINV_apply
  
  
negative.code:     implicit none
      type (hecmwST_matrix)     :: hecMAT
  
      integer(kind=kint) :: i, j, k, jS, jE, in, itr, NP, idof, jdof, iitr 
      real(kind=krealp) :: dd, dtmp(hecMAT%NDOF), YV(hecMAT%NDOF), X(hecMAT%NDOF)
  
      real(kind=krealp) :: FILTER, SIGMA_DIAG
      real(kind=krealp), allocatable :: zz(:), vv(:)
  
  