anchor.diff:        kl_pairs,                                                      &
         W, residual_norm,                                              &
         U, overlap,                                                    &
  -      R,                                                             &
  +      R, c,                                                          &
         h,                                                             &
         y,                                                             &
         lambda                                                         &
  
positive.code:       y(N_st_diag,davidson_sze_max,N_st_diag,davidson_sze_max),      &
        residual_norm(N_st_diag),                                      &
        overlap(N_st_diag,N_st_diag),                                  &
        lambda(N_st_diag*davidson_sze_max))
    
    ASSERT (N_st > 0)
  
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)
  