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;
} |
