anchor.diff:    integer                        :: k_pairs, kl
     
     integer                        :: iter2
  -  double precision, allocatable  :: W(:,:),  U(:,:), R(:,:), S(:,:)
  +  double precision, allocatable  :: W(:,:),  U(:,:), S(:,:)
     double precision, allocatable  :: y(:,:), h(:,:), lambda(:), s2(:)
     double precision, allocatable  :: c(:), s_(:,:), s_tmp(:,:)
     double precision               :: diag_h_mat_elem
  
positive.code:   allocate(                                                          &
        W(sze_8,N_st_diag*itermax),                           &
        U(sze_8,N_st_diag*itermax),                           &
        R(sze_8,N_st_diag),                                            &
        S(sze_8,N_st_diag*itermax),                           &
        h(N_st_diag*itermax,N_st_diag*itermax),      &
        y(N_st_diag*itermax,N_st_diag*itermax),      &
  
negative.code:   
    converged = .False.
    
    double precision               :: r1, r2
    do k=N_st+1,N_st_diag-2,2
        do i=1,sze
          call random_number(r1)
          call random_number(r2)
          r1 = dsqrt(-2.d0*dlog(r1))
          r2 = dtwo_pi*r2
          u_in(i,k) = r1*dcos(r2)
          u_in(i,k+1) = r1*dsin(r2)
        enddo
    enddo
    do k=N_st_diag-1,N_st_diag
        do i=1,sze
          call random_number(r1)
          call random_number(r2)
          r1 = dsqrt(-2.d0*dlog(r1))
          r2 = dtwo_pi*r2
          u_in(i,k) = r1*dcos(r2)
        enddo
    enddo
  
  