Commit c24edfd2 authored by Lotshaw, Phil's avatar Lotshaw, Phil
Browse files

code for the simulations

parent fc3e61fa

Code/BFGS_OMP_mod.f90

0 → 100644
+277 −0
Original line number Diff line number Diff line 
module BFGS



implicit none



contains



SUBROUTINE dfpmin(p,n,gtol,iter,fret,func,dfunc,graph_number)

!from Numerical Recipes in FORTRAN 77: The art of Scientific Programming

!https://www.goodreads.com/book/show/27815.Numerical_Recipes_in_FORTRAN_77

implicit none



INTEGER iter,n,NMAX,ITMAX,graph_number

Double precision fret,gtol,p(n),func,EPS,STPMX,TOLX

PARAMETER (NMAX=50,ITMAX=200,STPMX=6.d0,EPS=3.d-10,TOLX=4.d0*EPS) 

EXTERNAL dfunc,func

! USES dfunc,func,lnsrch



!assumes all values are of order unity



!Given a starting point p(1:n) that is a vector of length n, 

!the Broyden-Fletcher-Goldfarb-Shanno variant of Davidon-Fletcher-Powell 

!minimization is performed on a function func, using its gradient as calculated 

!by a routine dfunc. The convergence requirement on zeroing the gradient is 

!input as gtol. Returned quantities are p(1:n) (the location of the mini- mum), 

!iter (the number of iterations that were performed), 

!and fret (the minimum value of the function). The routine lnsrch is called to 

!perform approximate line minimizations. 

!Parameters: NMAX is the maximum anticipated value of n; 

!ITMAX is the maximum allowed number of iterations; 

!STPMX is the scaled maximum step length allowed in line searches; 

!TOLX is the convergence criterion on x values.



INTEGER i,its,j

LOGICAL check

Double precision den,fac,fad,fae,fp,stpmax,sum,sumdg,sumxi,temp,test, &

& dg(NMAX),g(NMAX),hdg(NMAX),hessin(NMAX,NMAX), pnew(NMAX),xi(NMAX)

!Calculate starting function value and gradient,



fp=func(n,p,graph_number)

call dfunc(n,p,g,func,graph_number)

sum=0.d0



!and initialize the inverse Hessian to the unit matrix.



do i=1,n

	do j=1,n

		hessin(i,j)=0.d0 

	enddo 

	hessin(i,i)=1.d0

	!Initial line direction.

	xi(i)=-g(i)

	sum=sum+p(i)**2

enddo 



stpmax=STPMX*max(sqrt(sum),float(n))



!Main loop over the iterations.

do its=1,ITMAX  

	iter=its

	!The new function evaluation occurs in lnsrch; save the function value in fp for the next

	!line search. It is usually safe to ignore the value of check. 

	call lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,check,func,graph_number)

	fp=fret

	do i=1,n 

		!Update the line direction, and the current point.

		xi(i)=pnew(i)-p(i) 

		p(i)=pnew(i)

	enddo 

	test=0.d0

	!Test for convergence on ∆x.

	do i=1,n

		temp=abs(xi(i))/max(abs(p(i)),1.d0)

		if(temp.gt.test)test=temp 

	enddo

	if(test.lt.TOLX)return 

	!Save the old gradient,

	do i=1,n

		dg(i)=g(i) 

	enddo

	!and get the new gradient.

	call dfunc(n,p,g,func,graph_number)

	!Test for convergence on zero gradient.

	test=0.d0

	den=max(fret,1.d0) 

	do i=1,n

		temp=abs(g(i))*max(abs(p(i)),1.d0)/den

		if(temp.gt.test)test=temp 

	enddo

	if(test.lt.gtol)return 

	!Compute difference of gradients,

	do i=1,n

		dg(i)=g(i)-dg(i) 

	enddo

	!and difference times current matrix.

	do i=1,n 

		hdg(i)=0.d0

		do j=1,n 

			hdg(i)=hdg(i)+hessin(i,j)*dg(j)

		enddo 

	enddo

	!Calculate dot products for the denominators.

	fac=0.d0

	fae=0.d0

	sumdg=0.d0

	sumxi=0.d0

 	do i=1,n 

 		fac=fac+dg(i)*xi(i) 

 		fae=fae+dg(i)*hdg(i) 

 		sumdg=sumdg+dg(i)**2 

 		sumxi=sumxi+xi(i)**2

	enddo 

	!Skip update if fac not sufficiently positive.

	if(fac.gt.sqrt(EPS*sumdg*sumxi))then

		fac=1.d0/fac

		fad=1.d0/fae

		!The vector that makes BFGS different from DFP:

		do i=1,n 

			dg(i)=fac*xi(i)-fad*hdg(i) 

		enddo

		!The BFGS updating formula: 

		do i=1,n 

			do j=i,n

				hessin(i,j)=hessin(i,j)+fac*xi(i)*xi(j) -fad*hdg(i)*hdg(j)+fae*dg(i)*dg(j)

				hessin(j,i)=hessin(i,j) 

			enddo

		enddo 

	endif



	!Now calculate the next direction to go,

	do i=1,n  

		xi(i)=0.d0

		do j=1,n 

			xi(i)=xi(i)-hessin(i,j)*g(j)

		enddo 

	enddo

!and go back for another iteration. 

enddo 

print *, "too many iterations in dfpmin"

return

END subroutine dfpmin





SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func,graph_number) 

!from Numerical Recipes in FORTRAN 77: The art of Scientific Programming

!https://www.goodreads.com/book/show/27815.Numerical_Recipes_in_FORTRAN_77

implicit none

INTEGER n, graph_number

LOGICAL check

DOUBLE PRECISION f,fold,stpmax,g(n),p(n),x(n),xold(n),func,ALF,TOLX 

PARAMETER (ALF=1.d-4,TOLX=1.d-7)

EXTERNAL func

! USESfunc

!Given an n-dimensional point xold(1:n), the value of the function and gradient there, 

!fold and g(1:n), and a direction p(1:n), finds a new point x(1:n) along the direction p 

!from xold where the function func has decreased “sufficiently.” The new function value is

!returned in f. stpmax is an input quantity that limits the length of the steps so that 

!you do not try to evaluate the function in regions where it is undefined or subject to 

!overflow. p is usually the Newton direction. The output quantity check is false on a normal 

!exit. It is true when x is too close to xold. In a minimization algorithm, this usually 

!signals convergence and can be ignored. However, in a zero-finding algorithm the calling 

!program should check whether the convergence is spurious.



!Parameters: ALF ensures sufficient decrease in function value; TOLX is the convergence

!criterion on ∆x.

INTEGER i

Double Precision a,alam,alam2,alamin,b,disc,f2,rhs1,rhs2,slope,sum,temp,test,tmplam

check=.false. 

sum=0.d0



do i=1,n

	sum=sum+p(i)*p(i) 

enddo

sum=sqrt(sum) 

!scale if attempted step is too big

if(sum.gt.stpmax) then

	do i=1,n 

		p(i)=p(i)*stpmax/sum

	enddo

endif

slope=0.d0 

do i=1,n

	slope=slope+g(i)*p(i) 

enddo

if (slope .ge. 0.d0) then

	print *, 'roundoff problem in lnsrch'

	print *, 'slope:',slope

	stop

endif

!Compute λmin.

test=0.d0 

do i=1,n

	temp=abs(p(i))/max(abs(xold(i)),1.d0)

	if(temp.gt.test)test=temp 

enddo

alamin=TOLX/test

!Always try full Newton step first.

alam=1.d0 

1 continue

	do i=1,n

		x(i)=xold(i)+alam*p(i) 

	enddo

	f=func(n,x,graph_number) 

	!Convergence on ∆x. For zero finding, the calling program should verify the convergence.

	if(alam.lt.alamin)then

		do i=1,n 

			x(i)=xold(i)

		enddo

		check=.true. 

		return

	!Sufficient function decrease.

	else if(f.le.fold+ALF*alam*slope)then 

		return

	!Backtrack. 

	else 

		!First time.

		if(alam.eq.1.d0)then

			tmplam=-slope/(2.d0*(f-fold-slope)) 

		!Subsequent backtracks.

		else

			rhs1=f-fold-alam*slope 

			rhs2=f2-fold-alam2*slope 

			a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2) 

			b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/(alam-alam2)

			if(a.eq.0.d0)then

				tmplam=-slope/(2.d0*b) 

			else

				disc=b*b-3.d0*a*slope 

				if (disc .lt. 0.d0) then 

					tmplam=0.5d0*alam

				else if(b.le.0.d0)then 

					tmplam=(-b+sqrt(disc))/(3.d0*a)

				else 

					tmplam=-slope/(b+sqrt(disc))

				endif 

			endif

			if(tmplam.gt.0.5d0*alam)tmplam=0.5d0*alam

		endif

	endif

	alam2=alam

	f2=f

	alam=max(tmplam,0.1d0*alam)

goto 1

end subroutine lnsrch







subroutine dfunc(n,x,grad,func,graph_number)

!from Numerical Recipes in FORTRAN 77: The art of Scientific Programming

!https://www.goodreads.com/book/show/27815.Numerical_Recipes_in_FORTRAN_77

	implicit none



	integer :: n,i,graph_number

	double precision :: x(n), x_dif_plus(n), x_dif_min(n)

	double precision :: nominal_delta = 1.d-6

	double precision :: h,temp

	double precision :: grad(n)

	double precision :: func



	External func



	!set h exact to numerical precision

	!so h can be represented exactly in base 2

	temp =  x(1) + nominal_delta

 	h = temp-x(1)



 	do i = 1,n

  		x_dif_plus = x

 		x_dif_min = x

 		x_dif_plus(i) = x(i) + h

 		x_dif_min(i) = x(i) - h

 		grad(i) = (func(n,x_dif_plus,graph_number) - func(n,x_dif_min,graph_number))/(2.d0*h)

	enddo



end subroutine dfunc





end module BFGS
 
 No newline at end of file

Code/QAOA_BFGS_OMP.f90

0 → 100644
+148 −0
Original line number Diff line number Diff line 
include'QAOA_parameters_OMP_mod.f90'

include'QAOA_subroutines_OMP_mod.f90'

include"BFGS_OMP_mod.f90"

Program QAOA

use parameters

use QAOA_subroutines

use BFGS

use omp_lib

implicit none



integer :: i,j,k,x,y,z,n,m,nits

integer :: iter,even,odd,tid

double precision :: fret

double precision :: timef,time0

double precision :: C,p_cz_max,C_test,angles_test(2*p_max)

double precision :: C_max_c, p_max_c, beta_max_c(p_max), gamma_max_c(p_max)

double precision :: prob_cmax

double precision :: angles(2*p_max)

double precision :: tmp,C_graph_avg,C_graph_stdev

double precision :: C_graph_avg_omp, C_graph_stdev_omp



integer :: BFGS_loops=0

integer :: n_threads



double precision :: grad(2*p_max), norm

character*1 :: Mixer='X'



call system('mkdir '//trim(save_folder))

call system('cp QAOA_BFGS_OMP.f90 '//trim(save_folder)//'QAOA_BFGS.f90')

call system('cp QAOA_parameters_OMP_mod.f90 '//trim(save_folder)//'QAOA_parameters_mod.f90')

call system('cp QAOA_subroutines_OMP_mod.f90 '//trim(save_folder)//'QAOA_subroutines_mod.f90')

call system('cp BFGS_OMP_mod.f90 '//trim(save_folder)//'BFGS_mod.f90')



n_threads=omp_get_max_threads()

call OMP_SET_NUM_THREADS(n_threads)

print *, 'n_threads:',n_threads

time0=omp_get_wtime()



!set up the basis for calculations

call enumerate_basis

call calc_pair_list_fast

call init_random_seed()



open(1,file=trim(graph_file),status='old')

call Count_num_graphs

close(1)



allocate(C_opt(graph_num_tot),Beta_opt(graph_num_tot,p_max),Gamma_opt(graph_num_tot,p_max), &

		& C0(graph_num_tot),cz_max_save(graph_num_tot),p_C_max(graph_num_tot), &

		& p0_C_max(graph_num_tot),good_loops(graph_num_tot),vertex_degrees(graph_num_tot,0:n_qubits-1), &

		& all_odd_degree(graph_num_tot),all_even_degree(graph_num_tot),edges(graph_num_tot,0:n_edges_max-1,0:1), &

		& n_edges(graph_num_tot),cz_vec(graph_num_tot,0:dim-1))



edges=0

n_edges=0

cz_vec=0

vertex_degrees=0

open(1,file=trim(graph_file),status='old')

call find_all_edges

close(1)





good_loops=0

vertex_degrees=0

C_opt=0.d0



open(123,file=trim(save_folder)//'C_graph_avg_BFGS',status='unknown')

open(579,file=trim(save_folder)//'degeneracies',status='unknown')



do BFGS_loops=1,min_good_loops



	graph_num=0

	c_graph_avg=0.d0

	C_graph_stdev=0.d0

	

	more_graphs=.true.



!$OMP PARALLEL PRIVATE(angles,fret,C,prob_Cmax,C_graph_avg_omp,C_graph_stdev_omp,tid,graph_num)

tid=omp_get_thread_num()

C_graph_avg_omp=0.d0

C_graph_stdev_omp=0.d0

!$OMP DO SCHEDULE(dynamic)

	do graph_num = 1,graph_num_tot

		if (BFGS_loops .eq. 1) call Calc_cz_vec(BFGS_loops,graph_num)



		call Generate_angles(2*p_max,angles,graph_num)



		call dfpmin(angles,2*p_max,gtol,iter,fret,QAOA_ExpecC,dfunc,graph_num)

		C=-fret*cz_max_save(graph_num)



		call Check_optimum_degeneracies_simple(2*p_max,angles,C,graph_num,tid)

		

		if (cz_max_save(graph_num) .gt. 1.d-8) then

			C_graph_avg_omp= C_graph_avg_omp + C_opt(graph_num)/cz_max_save(graph_num)

			C_graph_stdev_omp = C_graph_stdev_omp + ( C_opt(graph_num) / cz_max_save(graph_num) )**2

		endif

	enddo

!$OMP enddo



!$OMP CRITICAL



	C_graph_avg = C_graph_avg + C_graph_avg_omp

	C_graph_stdev = C_graph_stdev + C_graph_stdev_omp

!$OMP END CRITICAL

!$OMP END PARALLEL



	if (mod(BFGS_loops,10) .eq. 0) then

		timef=omp_get_wtime()

		print *, 'BFGS_loops,time:',BFGS_loops,timef-time0

	endif



	C_graph_avg=C_graph_avg/dble(graph_num_tot)

	C_graph_stdev=C_graph_stdev/dble(graph_num_tot)

	C_graph_stdev = dsqrt( c_graph_stdev - ( (c_graph_avg)**2 ) )



	write(123,*) BFGS_loops, C_graph_avg,C_graph_stdev



enddo



open(2,file=trim(save_folder)//'QAOA_dat',status='unknown')

do graph_num = 1,graph_num_tot

	if (all_even_degree(graph_num)) then

		even = 1

	else

		even = 0

	endif

	if (all_odd_degree(graph_num)) then

		odd = 1

	else

		odd = 0

	endif

	write(2,*) graph_num, cz_max_save(graph_num), C0(graph_num), C_opt(graph_num), p_C_max(graph_num), p_max, &

		& (beta_opt(graph_num,j)/pi,j=1,p_max), (gamma_opt(graph_num,j)/pi,j=1,p_max),even, &

		& odd

enddo



close(2)





print *, 'BFGS_loops:', BFGS_loops

timef=omp_get_wtime()

print*, 'elapsed time:', timef-time0



close(123)

close(579)



end program QAOA

Code/QAOA_parameters_OMP_mod.f90

0 → 100644
+33 −0
Original line number Diff line number Diff line 
module parameters



implicit none



character*200, parameter :: save_folder="Test_Graph6/"

character*200, parameter :: graph_file='Graphs/graph6c.txt'

integer, parameter :: n_qubits=6, p_max=3, smaller_p=1

integer, parameter :: loops_param_array(3) = (/ 50, 100, 500 /)

integer, parameter :: min_good_loops =loops_param_array(p_max)



integer, parameter :: dim=2**n_qubits, z_max=2**n_qubits-1,n_edges_max=n_qubits*(n_qubits-1)/2

double precision, parameter :: pi=dacos(-1.d0)

double precision, parameter :: gtol=1.d-10, degen_tol=1.d-7,convergence_tolerance=1.d-7

logical, parameter :: force_angle_ordering=.true.

integer :: graph_num_tot,graph_num

integer :: degeneracy_cz_max



double precision, parameter :: gamma_max = pi, beta_max=pi/4.d0

integer, parameter :: min_it_gamma=0, max_it_gamma=200!for brute force grid search

integer, parameter :: min_it_beta=0, max_it_beta=50



logical :: more_graphs=.true.



integer :: pairs_list(0:n_qubits-1,0:dim/2-1,0:1),basis(0:dim-1,0:n_qubits-1)



integer, allocatable :: n_angles_ma(:)

double precision :: beta(p_max),gamma(p_max),cz_max

double precision, allocatable :: cz_max_save(:), C0(:), C_opt(:), p_C_max(:), beta_opt(:,:), gamma_opt(:,:),p0_C_max(:)

integer, allocatable :: good_loops(:), vertex_degrees(:,:), edges(:,:,:),n_edges(:),cz_vec(:,:)

logical, allocatable :: all_odd_degree(:),all_even_degree(:)





end module parameters
 
 No newline at end of file

Code/QAOA_subroutines_OMP_mod.f90

0 → 100644
+363 −0
Original line number Diff line number Diff line 
include"QAOA_subroutines_UnitaryEv_mod.f90"

include"QAOA_subroutines_graphs_mod.f90"

module QAOA_subroutines

use QAOA_subroutines_UnitaryEv

use QAOA_subroutines_graphs

implicit none



contains





subroutine simplify_angles_even_odd_simple(n_angles,angles,result,graph_number)



	use parameters



	implicit none



	integer :: graph_number

	integer :: n_angles

	double precision :: angles(n_angles),angles_test(n_angles)

	double precision :: result,result_test

	integer :: i,j,k

	logical :: dif_angles,pass



	dif_angles=.false.

	angles_test = angles



	!assumes beta_1 < 0



	!if all even degree, set -pi/2 \leq \gamma \leq pi/2

	if (all_even_degree(graph_number)) then



		do i = 1,n_angles/2

			if (dabs(angles_test(n_angles/2+i)) .gt. pi/2.d0) then

				if (angles_test(n_angles/2+i) .gt. 0.d0) then

					angles_test(n_angles/2+i) = angles_test(n_angles/2+i) - pi

				else

					angles_test(n_angles/2+i) = angles_test(n_angles/2+i) + pi

				endif

				dif_angles=.true.

			endif

		enddo



	!if all odd degree, set -pi/2 \leq \gamma \leq pi/2

	elseif (all_odd_degree(graph_number)) then



		do i = 1,n_angles/2

			!map the gamma angle to have |gamma_1| < \pi/2 while maintaining beta_1 < 0 if (i .eq. 1)

			!this is a conjuction of the odd symemtry and the time reversal symmetry

			if (dabs(angles_test(n_angles/2+i)) .gt. pi/2.d0) then

				if (angles_test(n_angles/2+i) .lt. 0.d0) then

					angles_test(n_angles/2+i) = angles_test(n_angles/2+i) + pi

				else

					angles_test(n_angles/2+i) = angles_test(n_angles/2+i) - pi

				endif

				do j = i,n_angles/2

					angles_test(j) = -angles_test(j)

				enddo

				if (i .eq. 1) angles_test=-angles_test

				dif_angles=.true.

			endif

		enddo



	endif



	if (dif_angles) then

		call Check_new_angle_symmetry(angles_test,result,result_test,pass,graph_number)

		if (pass) then

			angles = angles_test

		else

			print *, 'failed even/odd symmetry'	

			print *, 'result, result_test:',result,result_test

			stop

		endif

	endif



end subroutine simplify_angles_even_odd_simple





subroutine simplify_angles_simple(angles,result,graph_number)



	use parameters

	implicit none



	integer :: graph_number

	integer :: tmp_int

	double precision :: angles(2*p_max),angles_test(2*p_max)

	double precision :: result,result_test

	integer :: i

	logical :: pass



	!put angles into the intervals beta \in [-pi/4,pi/4], gamma \in [-pi,pi]

	angles_test=angles

	do i = 1,p_max

		if ( (angles_test(i) .lt. -0.25d0*pi) .or. (angles_test(i) .gt. 0.25d0*pi) ) then

			if (angles_test(i) .lt. -0.25d0*pi) then

				tmp_int = floor( dabs( (angles_test(i)-0.25*pi)/(pi/2.d0) ) )

				angles_test(i) = angles_test(i) + dble(tmp_int)*(pi/2.d0)

			else

				tmp_int = floor( dabs( (angles_test(i)+0.25*pi)/(pi/2.d0) ) )

				angles_test(i) = angles_test(i) - dble(tmp_int)*(pi/2.d0)

			endif

			if ( (angles_test(i) .lt. -0.25d0*pi) .or. (angles_test(i) .gt. 0.25d0*pi) ) then

				print *, 'error in angles_test!'

				print *, 'angles_test(i)/pi',angles_test(i)/pi

				stop

			endif

		endif

		if ( (angles_test(p_max+i) .lt. -pi) .or. (angles_test(p_max+i) .gt. pi)) then

			if (angles_test(p_max+i) .lt. -pi) then

				tmp_int = floor( dabs( (angles_test(p_max+i)-pi)/(pi*2.d0) ) )

				angles_test(p_max+i) = angles_test(p_max+i) + dble(tmp_int)*(pi*2.d0)

			else

				tmp_int = floor( dabs( (angles_test(p_max+i)+pi)/(pi*2.d0) ) )

				angles_test(p_max+i) = angles_test(p_max+i) - dble(tmp_int)*(pi*2.d0)

			endif

			if ( (angles_test(p_max+i) .lt. -pi) .or. (angles_test(i) .gt. pi) ) then

				print *, 'error in angles_test!'

				print *, 'angles_test(p_max+i)', angles_test(p_max+i)

				stop

			endif

		endif

	enddo



	if (angles_test(1) .gt. 0.d0) angles_test=-angles_test



	call Simplify_angles_even_odd_simple(2*p_max,angles_test,result,graph_number)



	!test these give the same result

	call Check_new_angle_symmetry(angles_test,result,result_test,pass,graph_number)



	if (.not. pass) then

		print *, 'failed to simplify angles'

		print *, 'result,result_test',result,result_test

		print *,'angles, angles_test:',angles,angles_test

		stop 

	else

		angles = angles_test

	endif



end subroutine simplify_angles_simple





subroutine Check_optimum_degeneracies_simple(n_angles,angles,result,graph_number,tid)



	use parameters 



	implicit none



	integer :: graph_number,tid

	integer :: i,j,n_angles

	logical :: degenerate_opt,identical_angles

	double precision :: result, angles(2*p_max)

	integer :: counts_opt,counts_new

	identical_angles=.true.

	degenerate_opt=.false.





	call simplify_angles_simple(angles,result,graph_number)



	if ( (all_even_degree(graph_number) .or. all_odd_degree(graph_number)) .and. &

			& dabs(angles(n_angles/2+1)) .gt. pi/2.d0) then

		print *, 'graph_number, all_even, all_odd:',graph_number, all_even_degree(graph_number), all_odd_degree(graph_number)

		print *, 'angles/pi:',(angles(i)/pi,i=1,n_angles)

		print *, 'angles > pi/2!'

		stop

	endif



	!check for degeneracy

	if (dabs(result - C_opt(graph_number)) .lt. degen_tol) degenerate_opt=.true.

	!if degenerate then try simplifying, check if the angles are the same

	if (degenerate_opt) then

		do i = 1,p_max

			if ( (dabs(angles(i) - beta_opt(graph_number,i)) .gt. 1.d-4) .or. &

					& ( dabs(angles(i+p_max) - gamma_opt(graph_number,i)) .gt. 1.d-4) ) then

				identical_angles=.false.

			endif

		enddo



		if (force_angle_ordering) then

			if ( (.not. identical_angles) ) then

				!choose the angles that best line up with the angle patterns

				!do this by counting how many go into the octant of parameter space 

				!-0.25pi < beta < 0, -pi/2 < gamma < 0

				counts_opt=0

				counts_new=0

				do i = 1,p_max

					if (beta_opt(graph_number,i) .le. 0.d0 .and. gamma_opt(graph_number,i) .le. 0.d0 &

						& .and. gamma_opt(graph_number,i) .ge. -pi/2.d0) then

						counts_opt = counts_opt+1

					endif

					if (angles(i) .le. 0.d0 .and. angles(n_angles/2+i) .le. 0.d0 .and. angles(n_angles/2+i) .ge. -pi/2.d0) then

						counts_new = counts_new+1

					endif

				enddo

				if (counts_new .gt. counts_opt) then

					write(579,*) graph_number,C_opt(graph_number),&

						& (beta_opt(graph_number,j)/pi,j=1,p_max),(gamma_opt(graph_number,j)/pi,j=1,p_max), &

						result, (angles(j)/pi,j=1,2*p_max)

					print *, 'graph_number,C_opt,C',graph_number,C_opt(graph_number),result

					print *, 'degenerate angles/pi:'

					print *, 'counts_old, counts_new:',counts_opt,counts_new

					print *, (beta_opt(graph_number,j)/pi,j=1,p_max),(gamma_opt(graph_number,j)/pi,j=1,p_max)

					print *, (angles(j)/pi,j=1,2*p_max)

					print *, ' '

					beta_opt(graph_number,:) = angles(1:p_max)

					gamma_opt(graph_number,:) = angles(p_max+1:2*p_max)

					C_opt(graph_number)=result

					p_C_max(graph_number) = -QAOA_pmax(2*p_max,angles,graph_number)

				endif

			endif

		endif

	endif



	if (result .gt. (C_opt(graph_number)+convergence_tolerance) ) then



    	C_opt(graph_number) = result

    	Beta_opt(graph_number,1:p_max) = angles(1:p_max)

    	Gamma_opt(graph_number,1:p_max) = angles(p_max+1:2*p_max)

    	p_C_max(graph_number) = -QAOA_pmax(2*p_max,angles,graph_number)

    	good_loops(graph_number) = 0

    else

    	good_loops(graph_number) = good_loops(graph_number)+1

    endif



end subroutine Check_optimum_degeneracies_simple





subroutine Generate_angles(n_angles,angles,graph_number)



	use parameters

	implicit none



	integer :: i

	integer :: n_angles,graph_number

	double precision :: angles(n_angles)

	double precision :: tmp



	do i = 1,p_max

		call random_number(tmp)

		angles(i) = pi/2.d0*(tmp-0.5d0)

		call random_number(tmp)

		angles(p_max+i) = 2.d0*pi*(tmp-0.5d0)

	enddo



end subroutine Generate_angles





subroutine set_angles(x,y,angles)



	!for brute force search of angles:

	!expand x and y in base max_it_gammma and max_it_beta 

	!to map x,y to a p-digit number

	!The p-digit number is then mapped to a point on mesh

	!for the numerical search over angles

	use parameters

	implicit none



	integer :: x,y

	integer :: y0,x0

	integer :: i,j

	double precision :: angles(2*p_max)



	beta=0.d0

	gamma=0.d0

	!write beta(i), gamma(i) on a grid of points

	do i= p_max, 1, -1

		y0=0.d0

		x0=0.d0

		if (i .lt. p_max) then

			do j = p_max, i, -1

				y0 = y0 + beta(j)*(max_it_beta**(j-1))

				x0 = x0 + gamma(j)*(max_it_gamma**(j-1))

			enddo

		endif

		beta(i) = (y-y0)/(max_it_beta**(i-1))

		gamma(i) = ((x-x0)/max_it_gamma**(i-1))

	enddo



	!map the beta(i),gamma(i) grid to the the QAOA angle intervals

	do i = 1,p_max

		beta(i)=pi/2*(dble(beta(i))/dble(max_it_beta-1)-0.5d0)

		gamma(i)=pi*2*(dble(gamma(i))/dble(max_it_gamma-1)-0.5d0)

		angles(i) = beta(i)

		angles(p_max+i) = gamma(i)

	enddo





end subroutine set_angles







SUBROUTINE init_random_seed()

            INTEGER :: i, n, clock

            INTEGER, DIMENSION(:), ALLOCATABLE :: seed

          

            CALL RANDOM_SEED(size = n)

            ALLOCATE(seed(n))

          

            CALL SYSTEM_CLOCK(COUNT=clock)

          

            seed = clock + 37 * (/ (i - 1, i = 1, n) /)

            CALL RANDOM_SEED(PUT = seed)

          

            DEALLOCATE(seed)

END SUBROUTINE init_random_seed





subroutine Calc_cz_vec(BFGS_loops,graph_number)

    !calculate a vector c_z

    !with entries c_z[z] that equal the cost function evaluation for

    !the binary expansion of z

    use parameters

    implicit none

    integer :: z,m

    integer :: BFGS_loops,graph_number

    complex*16 :: state(0:z_max)



    do z =0,2**n_qubits-1

        do m = 0,n_edges(graph_number)-1

            if (basis(z,edges(graph_number,m,0)) .ne. basis(z,edges(graph_number,m,1))) then

            	cz_vec(graph_number,z) = cz_vec(graph_number,z) + 1

            endif

        enddo

    enddo



	if (BFGS_loops .eq. 1) then

		!calculate C0, the <C> for the initial state

		state=dcmplx(1.d0/dsqrt(dble(dim)),0.d0)

		call Calc_expec_C(state,C0(graph_number),p0_C_max(graph_number),graph_number)

		cz_max_save(graph_number) = maxval(cz_vec(graph_number,:))



		call calc_vertex_degrees(graph_number)



		if (all_odd_degree(graph_number)) print *, graph_number,'all odd degree'

		if (all_even_degree(graph_number)) print *, graph_number,'all even degree'



	endif



end subroutine Calc_cz_vec





subroutine Check_new_angle_symmetry(angles,result,result_test,pass,graph_number)



	use parameters

	implicit none



	integer :: graph_number

	double precision :: angles(2*p_max)

	double precision :: result, result_test

	logical :: pass



	result_test = -QAOA_ExpecC(2*p_max,angles,graph_number)*cz_max_save(graph_number)



	if (dabs(result-result_test) .gt. degen_tol) then

		pass = .false.

	else

		pass = .true.

	endif



end subroutine Check_new_angle_symmetry





end module QAOA_subroutines
 
 No newline at end of file

Code/QAOA_subroutines_UnitaryEv_mod.f90

0 → 100644
+132 −0

File added.

Preview size limit exceeded, changes collapsed.