anchor.diff:    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),      &
  
positive.code:       call dgemm('N','N', sze, N_st_diag, shift2,                    &
            1.d0, S, size(S,1), y, size(y,1), 0.d0, S(1,shift2+1), size(S,1))
  
        ! Compute residual vector
        ! -----------------------
        
        do k=1,N_st_diag
          do i=1,sze
            R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) )      &
                * (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz)
          enddo
          if (k <= N_st) then
            residual_norm(k) = u_dot_u(R(1,k),sze)
            to_print(1,k) = lambda(k) + nuclear_repulsion
            to_print(2,k) = s2(k)
            to_print(3,k) = residual_norm(k)
  
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
  
  