1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
| unsigned long graphLayout(Graph * graph, double * coords, unsigned long maxiter)
{
assert(graph);
assert(coords);
/* main simulation loop */
double energy = 0.0, _energy = +INFINITY;
double maxmove = -INFINITY;
double * _coords = malloc(graph->nvtxs*DIM*sizeof(double)); assert(_coords);
double * position = coords;
double * _position = _coords;
memcpy(_position, position, graph->nvtxs*DIM*sizeof(double));
int converged = 0;
double step = 1.0;
unsigned long iter = 0;
if(maxiter == 0) maxiter = ULONG_MAX;
if(period == 0) period = ULONG_MAX;
/* time loop */
while(!converged && iter < maxiter) {
int i;
energy = 0.0;
maxmove = -INFINITY;
/* for all particles xi, compute forces applied to it */
#pragma omp parallel for private(i)
for(int i = 0 ; i < graph->nvtxs ; i++) {
/* position of particle i */
double * _xi = _position+i*DIM;
/* reset force applied to particle i */
double force[DIM];
for(int d = 0 ; d < DIM ; d++) force[d] = 0.0;
/* compute repulsive forces (electrical: f=-C.K^2/|xi-xj|.Uij) */
for(int j = 0 ; j < graph->nvtxs ; j++) {
if(i == j) continue;
double * _xj = _position+j*DIM;
double dist = DISTANCE(_xi,_xj);
// power used for repulsive force model (standard is 1/r, 1/r^2 works well)
// double coef = 0.0; -C*K*K/dist; // power 1/r
double coef = -C*K*K*K/(dist*dist); // power 1/r^2
#pragma omp reduction (+:force)
for(int d = 0 ; d < DIM ; d++) force[d] += coef*(_xj[d]-_xi[d])/dist;
}
/* compute attractive forces (spring: f=|xi-xj|^2/K.Uij) */
for(int k = graph->xadj[i] ; k < graph->xadj[i+1] ; k++) {
int j = graph->adjncy[k]; /* edge (i,j) */
double * _xj = _position+j*DIM;
double dist = DISTANCE(_xi,_xj);
double coef = dist*dist/K;
for(int d = 0 ; d < DIM ; d++) force[d] += coef*(_xj[d]-_xi[d])/dist;
}
/* update particle position & energy */
double * xi = position+i*DIM;
double nforce2 = NORM2(force);
double move[DIM];
for(int d = 0 ; d < DIM ; d++) move[d] = step*force[d]/sqrt(nforce2);
for(int d = 0 ; d < DIM ; d++) xi[d] = _xi[d] + move[d];
double nmove = NORM(move);
/* reduction */
if(nmove > maxmove) maxmove = nmove; // ???
//#pragma omp reduction (+:energy)
energy += nforce2;
}
// energy /= graph->nvtxs; // normalize energy
/* update step length */
if(energy < _energy/5.0) step /= STEPCOEF; /* step length increased, if energy strongly reduced */
else if(energy < _energy) step = step; /* step length unchanged, if energy reduced */
else step *= STEPCOEF; /* step length reduced, if energy increased */
#pragma omp single
/* check convergence */
// if((energy < _energy) && ((_energy-energy)/energy < 0.0001)) converged = 1;
if(maxmove < K*TOL) converged = 1;
_energy = energy;
if((iter % period == 0) || converged) {
char vtkfile[255];
sprintf(vtkfile, "/tmp/output-%06ld.vtk", iter);
saveVTKGraph(vtkfile, graph, position);
printf("iter %ld: energy %.2f, maxmove %.6f\n", iter, energy, maxmove);
}
iter++;
/* swap former and newer position */
double * tmp = _position; _position = position; position = tmp;
} // end of while loop
free(_coords);
return iter;
} |
Partager