anchor.diff:        ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
         ! -------------------------------------------
   
  -!
  +
   !      do l=1,N_st_diag
   !        do k=1,N_st_diag
   !          do iter2=1,iter-1
  
positive.code:           U(i,k,iter+1) = 0.d0
            W(i,k,iter+1) = 0.d0
          enddo
           do iter2=1,iter
            do l=1,N_st_diag
              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
            enddo
          enddo
        enddo
  
  
  !      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
  
negative.code: !        enddo
  !      enddo
  
        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
        ! -------------
  