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