anchor.diff:        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
  
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:       ! 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
  