#include <stdio.h> #include <stdlib.h> #include <math.h> #include <string.h> #define N 3000 typedef unsigned char mtype; mtype mij0[N][N]; // for compressed row matrix storage int nnz[N]; unsigned short *jnz[N]; mtype *mij[N]; // Number of switch #define MAXQ 6 unsigned char spin[N]; unsigned char spin_bak[N]; int hist[MAXQ]; int e_total, e_best, e_bound; // For exp(-beta e) table #define EMAX 1024 int etab[EMAX]; // Probability table precision #define PMAX (1<<16) void init_etab(double beta) { int i; for(i=0;i<EMAX;i++) etab[i]=(int)(PMAX*exp(-beta*i)); } // Cal. effective field void cal_h(int i, int h[MAXQ]) { int j,q; unsigned short *jptr = jnz[i]; mtype *mptr = mij[i]; for(q=0;q<MAXQ;q++) h[q]= hist[q]; h[spin[i]]--; #if 0 for(j=0;j<N;j++)h[spin[j]] -= mij0[i][j]; #else // Hot spot while(*mptr != 0) { h[ spin[*jptr] ] -= *mptr; jptr++; mptr++; } #endif } // Heatbath update void do_hb(int i, int zero) { int j, q, q0, si_new; int h[MAXQ], hmin; int p[MAXQ], psum; cal_h(i, h); hmin=h[0]; q0=0; for(q=1;q<MAXQ;q++) { if(hmin>h[q]) { hmin=h[q]; q0=q; } } if(zero==1) { // Quench si_new=q0; } else { psum=0; for(q=0;q<MAXQ;q++) { int e = h[q]-hmin; p[q] = (e>=EMAX) ? 0:etab[e]; psum += p[q]; } psum = (int)(psum*1.0*rand()/RAND_MAX); q=0; while(psum >p[q]) { psum -= p[q]; q++; } si_new = q; } hist[spin[i]]--; hist[si_new]++; e_total += -h[spin[i]]+h[si_new]; spin[i]=si_new; if(e_total < e_best) { e_best = e_total; memcpy(spin_bak,spin,N); } } void read_file(char *fn) { char buf[N*2+100]; FILE *fp=fopen(fn,"r"); int i,j; fgets(buf,999,fp); e_bound = 0; for(i=0;i<N;i++) { int nz=0; fgets(buf,N+10,fp); for(j=0;j<N;j++) { mij0[i][j] = buf[j]-'0'; if(mij0[i][j]!=0 && i!=j) nz++; } nnz[i]=nz; mij0[i][i]=0; e_bound += nz; } fclose(fp); e_bound/=2; // Prepare compressed row data for(i=0;i<N;i++) { int nz=nnz[i]; jnz[i] = (unsigned short *)malloc(sizeof(unsigned short)*nz); mij[i] = (mtype *)malloc(sizeof(mtype)*(nz+1)); nz=0; for(j=0;j<N;j++) { if(mij0[i][j]!=0) { jnz[i][nz] = j; mij[i][nz] = mij0[i][j]; nz++; } } mij[i][nz] =0; } } int ene() { int c1=0; int c2=0; int i,j; for(i=0;i<N;i++) { for(j=0;j<i;j++) { c1 += (spin[i]==spin[j])? 0:mij0[i][j]; c2 += (spin[i]==spin[j])? 1:0; }} return c1+c2; } // Mean field critical beta double beta_mf(double m_ave) { int i,j; return MAXQ*1.0/N/m_ave; } void print_m() { int q; for(q=0;q<MAXQ;q++) printf(" %d",hist[q]); } int main(int argc, char **argv) { int i,j,q,step, maxstep; double beta, beta_min, beta_max, beta_inc; double beta_c; read_file(argv[1]); maxstep = atoi(argv[2]); beta_c = beta_mf(9.0); printf("#MF Beta_c = %f\n", beta_c); beta=beta_min=0.04; beta_max = 2.0; beta_inc = exp(log(beta_max/beta_min)/maxstep); for(q=0;q<MAXQ;q++) hist[q]=0; for(i=0;i<N;i++) { spin[i] = (rand()>>8)%MAXQ; hist[spin[i]]++; } e_total = e_best = ene(); // Main loop for(step=0;step<maxstep;step++) { init_etab(beta); for(i=0;i<N;i++) { do_hb(i,0); } if(step%10==0) { int k; printf("%e %d %d\n",beta,e_best, e_best-e_bound); } if (e_total == e_bound) break; beta *= beta_inc; } // Quench from the best config memcpy(spin,spin_bak,N); for(q=0;q<MAXQ;q++) hist[q]=0; for(i=0;i<N;i++) hist[spin[i]]++; for(step=0;step<40;step++) { for(i=0;i<N;i++) do_hb(i,1); } printf("#%d steps %d\n",maxstep,ene()); printf("# CSize "); print_m(); printf("\n"); //put_image("matrix.ppm", spin, mij0); }