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:       call random_number(r1)
        call random_number(r2)
        u_in(i,k) = dsqrt(-2.d0*dlog(r1))*dcos(dtwo_pi*r2)
  !      call RANDOM_NUMBER(u_in(i,k))
  !      u_in(i,k) = u_in(i,k) - 0.5d0
      enddo
      
      ! Gram-Schmidt
  
negative.code:           U(i,k,iter+1) = 0.d0
            W(i,k,iter+1) = 0.d0
          enddo
        enddo
  !      do k=1,N_st_diag
  !         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
  