anchor.diff:        ! Express eigenvectors of h in the determinant basis
         ! --------------------------------------------------
         
  -     !$OMP PARALLEL DEFAULT(NONE)                                     &
  -        !$OMP PRIVATE(k,i,l,iter2) &
  -        !$OMP SHARED(U,W,R,y,iter,lambda,N_st_diag,sze,to_print, &
  -        !$OMP residual_norm,nuclear_repulsion,N_st)
  +      call dgemm('N','N', sze, N_st_diag, N_st_diag*iter, &
  +       1.d0, U, size(U,1), y, size(y,1)*size(y,2), 0.d0, U(1,1,iter+1), size(U,1))
  +      call dgemm('N','N',sze,N_st_diag,N_st_diag*iter, &
  +       1.d0, W, size(W,1), y, size(y,1)*size(y,2), 0.d0, W(1,1,iter+1), size(W,1))
  +
  +      ! Compute residual vector
  +      ! -----------------------
  +      
         do k=1,N_st_diag
  -        !$OMP DO
  -        do i=1,sze
  -          U(i,k,iter+1) = 0.d0
  -          W(i,k,iter+1) = 0.d0
  -        enddo
  -        !$OMP END DO
  -         do iter2=1,iter
  -          do l=1,N_st_diag
  -            !$OMP DO
  -            do i=1,sze
  -              U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
  -              W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
  -            enddo
  -            !$OMP END DO NOWAIT
  -          enddo
  -        enddo
  -      
  -        ! Compute residual vector
  -        ! -----------------------
  -      
  -        !$OMP DO
           do i=1,sze
             R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1)
           enddo
  -        !$OMP END DO
  -        !$OMP SINGLE
           if (k <= N_st) then
             residual_norm(k) = u_dot_u(R(1,k),sze)
             to_print(1,k) = lambda(k) + nuclear_repulsion
             to_print(2,k) = residual_norm(k)
           endif
  -        !$OMP END SINGLE
         enddo
  -      !$OMP END PARALLEL
         
         write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))')  iter, to_print(:,1:N_st)
         call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged)
  
positive.code:       ! Gram-Schmidt
        ! ------------
        
        double precision               :: c(N_st_diag)
        do k=1,N_st_diag
          do iter2=1,iter
            call dgemv('T',sze,N_st_diag,1.d0,U(1,1,iter2),size(U,1),  &
                U(1,k,iter+1),1,0.d0,c,1)
            call dgemv('N',sze,N_st_diag,-1.d0,U(1,1,iter2),size(U,1), &
                c,1,1.d0,U(1,k,iter+1),1)
          enddo
          call dgemv('T',sze,N_st_diag,1.d0,U(1,1,iter+1),size(U,1),   &
              U(1,k,iter+1),1,0.d0,c,1)
          call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1),        &
              c,1,1.d0,U(1,k,iter+1),1)
          call normalize( U(1,k,iter+1), sze )
        enddo
  
      
        !DEBUG : CHECK OVERLAP
        !print *,  '==='
        !do k=1,iter+1
        !  do l=1,k
        !  c = u_dot_v(U(1,1,k),U(1,1,l),sze)
        !  print *,  k,l, c
        !  enddo
        !enddo
        !print *,  '==='
        !pause
        !END DEBUG
  
  
      enddo
  
  
negative.code:       ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
        ! -------------------------------------------
  
        call dgemm('T','N', N_st_diag*iter, N_st_diag, sze,               &
            1.d0, U, size(U,1), W(1,1,iter), size(W,1),   &
            0.d0, h(1,1,1,iter), size(h,1)*size(h,2))
  
        ! Diagonalize h
        ! -------------
        call lapack_diag(lambda,y,h,N_st_diag*davidson_sze_max,N_st_diag*iter)
  