anchor.diff:  
     if(lambda_type == 0) then
       get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint)
  +    get_dij_index = get_dij_index * rho_mrcc(i,s)
     else
       call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi)
       get_dij_index = HIi * lambda_mrcc(s, i)
  
positive.code:          print *, "residu ", k, res
         end if
         
         if(res < 1d-16) exit
       end do
       
       norm = 0.d0
  
negative.code:      
       if(AtA_size > size(AtA_val)) stop "SIZA"
       print *, "ATA SIZE", ata_size
       do i=1,nex
         x(i) = AtB(i)
       enddo
  