146 printf(
"%.2lf\t",m.
m[i][j]);
160 printf(
"%.2lf\t",v.
v[i]);
169 #define clearV(vv) clearVInternal((vv).v)
170 #define clearVInternal(v) \
171 ((v)[0] = (v)[1] = (v)[2] = 0.0 )
176 #define cpyV(r,a) cpyVInternal((r).v,(a).v)
177 #define cpyVInternal(r,a) \
178 {(r)[0] = (a)[0]; (r)[1] = (a)[1]; (r)[2] = (a)[2]; }
183 #define addV(r,a,b) addVInternal((r).v,(a).v,(b).v)
184 #define addVInternal(r,a,b) \
185 {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1]; (r)[2] = (a)[2] + (b)[2]; }
191 #define subV(r,a,b) subVInternal((r).v,(a).v,(b).v)
192 #define subVInternal(r,a,b) \
193 {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1]; (r)[2] = (a)[2] - (b)[2]; }
198 #define norm(vv) normInternal((vv).v)
199 #define normInternal(v) \
200 (sqrt( pow((v)[0],2) + pow((v)[1],2) + pow((v)[2],2) ))
205 #define multSV(r,s,vv) multSVInternal((r).v,s,(vv).v)
206 #define multSVInternal(r,s,v) \
207 {(r)[0] = (s) * (v)[0]; (r)[1] = (s) * (v)[1]; (r)[2] = (s) * (v)[2]; }
212 #define clearM(mm) clearMInternal((mm).m)
213 #define clearMInternal(m) \
214 ((m)[0][0] = (m)[0][1] = (m)[0][2] = \
215 (m)[1][0] = (m)[1][1] = (m)[1][2] = \
216 (m)[2][0] = (m)[2][1] = (m)[2][2] = 0.0 )
228 #define det(mm) detInternal((mm).m)
229 #define detInternal(m) \
230 ( (m)[0][0] * (m)[1][1] * (m)[2][2] \
231 + (m)[0][1] * (m)[1][2] * (m)[2][0] \
232 + (m)[0][2] * (m)[1][0] * (m)[2][1] \
233 - (m)[0][0] * (m)[1][2] * (m)[2][1] \
234 - (m)[0][1] * (m)[1][0] * (m)[2][2] \
235 - (m)[0][2] * (m)[1][1] * (m)[2][0] )
240 #define init_matrix(mm,a,b,c,d,e,f,g,h,i) init_matrixInternal((mm).m,a,b,c,d,e,f,g,h,i)
241 #define init_matrixInternal(m,a,b,c,d,e,f,g,h,i) \
242 {(m)[0][0] = a; (m)[0][1] = b; (m)[0][2] = c; \
243 (m)[1][0] = d; (m)[1][1] = e; (m)[1][2] = f; \
244 (m)[2][0] = g; (m)[2][1] = h; (m)[2][2] = i; }
249 #define multSM(r,s,mm) multSMInternal((r).m,s,(mm).m)
250 #define multSMInternal(r,s,m) \
251 {(r)[0][0] = (s) * (m)[0][0]; (r)[0][1] = (s) * (m)[0][1]; (r)[0][2] = (s) * (m)[0][2]; \
252 (r)[1][0] = (s) * (m)[1][0]; (r)[1][1] = (s) * (m)[1][1]; (r)[1][2] = (s) * (m)[1][2]; \
253 (r)[2][0] = (s) * (m)[2][0]; (r)[2][1] = (s) * (m)[2][1]; (r)[2][2] = (s) * (m)[2][2]; }
258 #define multSI(r,s) multSIInternal((r).m,s)
259 #define multSIInternal(r,s) \
260 {(r)[0][0] = (s); (r)[0][1] = 0; (r)[0][2] = 0; \
261 (r)[1][0] = 0; (r)[1][1] = (s); (r)[1][2] = 0; \
262 (r)[2][0] = 0; (r)[2][1] = 0; (r)[2][2] = (s); }
268 #define addM(r,a,b) addMInternal((r).m,(a).m,(b).m)
269 #define addMInternal(r,a,b) \
270 {(r)[0][0] = (a)[0][0] + (b)[0][0]; (r)[0][1] = (a)[0][1] + (b)[0][1]; (r)[0][2] = (a)[0][2] + (b)[0][2]; \
271 (r)[1][0] = (a)[1][0] + (b)[1][0]; (r)[1][1] = (a)[1][1] + (b)[1][1]; (r)[1][2] = (a)[1][2] + (b)[1][2]; \
272 (r)[2][0] = (a)[2][0] + (b)[2][0]; (r)[2][1] = (a)[2][1] + (b)[2][1]; (r)[2][2] = (a)[2][2] + (b)[2][2]; }
278 #define multVVt(r,a,b) multVVtInternal((r).m,(a).v,(b).v)
279 #define multVVtInternal(r,a,b) \
280 {(r)[0][0] = (a)[0] * (b)[0]; (r)[0][1] = (a)[0] * (b)[1]; (r)[0][2] = (a)[0] * (b)[2]; \
281 (r)[1][0] = (a)[1] * (b)[0]; (r)[1][1] = (a)[1] * (b)[1]; (r)[1][2] = (a)[1] * (b)[2]; \
282 (r)[2][0] = (a)[2] * (b)[0]; (r)[2][1] = (a)[2] * (b)[1]; (r)[2][2] = (a)[2] * (b)[2]; }
288 #define multMV(r,mm,vv) multMVInternal((r).v,(mm).m,(vv).v)
289 #define multMVInternal(r,m,v) \
290 {(r)[0] = (m)[0][0] * (v)[0] + (m)[0][1] * (v)[1] + (m)[0][2] * (v)[2]; \
291 (r)[1] = (m)[1][0] * (v)[0] + (m)[1][1] * (v)[1] + (m)[1][2] * (v)[2]; \
292 (r)[2] = (m)[2][0] * (v)[0] + (m)[2][1] * (v)[1] + (m)[2][2] * (v)[2]; }
298 #define inv(r,mm) invInternal((r).m,(mm).m)
299 #define invInternal(r,m) \
300 {(r)[0][0] = (m)[1][1] * (m)[2][2] - (m)[1][2] * (m)[2][1]; \
301 (r)[0][1] = - (m)[0][1] * (m)[2][2] + (m)[0][2] * (m)[2][1]; \
302 (r)[0][2] = (m)[0][1] * (m)[1][2] - (m)[0][2] * (m)[1][1]; \
303 (r)[1][0] = - (m)[1][0] * (m)[2][2] + (m)[1][2] * (m)[2][0]; \
304 (r)[1][1] = (m)[0][0] * (m)[2][2] - (m)[0][2] * (m)[2][0]; \
305 (r)[1][2] = - (m)[0][0] * (m)[1][2] + (m)[0][2] * (m)[1][0]; \
306 (r)[2][0] = (m)[1][0] * (m)[2][1] - (m)[1][1] * (m)[2][0]; \
307 (r)[2][1] = - (m)[0][0] * (m)[2][1] + (m)[0][1] * (m)[2][0]; \
308 (r)[2][2] = (m)[0][0] * (m)[1][1] - (m)[0][1] * (m)[1][0]; \
309 multSMInternal(r,1/detInternal(m),r); }
314 #define multMM(r,a,b) multMMInternal((r).m,(a).m,(b).m)
315 #define multMMInternal(r,a,b) \
316 {(r)[0][0] = (a)[0][0] * (b)[0][0] + (a)[0][1] * (b)[1][0] + (a)[0][2] * (b)[2][0]; \
317 (r)[0][1] = (a)[0][0] * (b)[0][1] + (a)[0][1] * (b)[1][1] + (a)[0][2] * (b)[2][1]; \
318 (r)[0][2] = (a)[0][0] * (b)[0][2] + (a)[0][1] * (b)[1][2] + (a)[0][2] * (b)[2][2]; \
319 (r)[1][0] = (a)[1][0] * (b)[0][0] + (a)[1][1] * (b)[1][0] + (a)[1][2] * (b)[2][0]; \
320 (r)[1][1] = (a)[1][0] * (b)[0][1] + (a)[1][1] * (b)[1][1] + (a)[1][2] * (b)[2][1]; \
321 (r)[1][2] = (a)[1][0] * (b)[0][2] + (a)[1][1] * (b)[1][2] + (a)[1][2] * (b)[2][2]; \
322 (r)[2][0] = (a)[2][0] * (b)[0][0] + (a)[2][1] * (b)[1][0] + (a)[2][2] * (b)[2][0]; \
323 (r)[2][1] = (a)[2][0] * (b)[0][1] + (a)[2][1] * (b)[1][1] + (a)[2][2] * (b)[2][1]; \
324 (r)[2][2] = (a)[2][0] * (b)[0][2] + (a)[2][1] * (b)[1][2] + (a)[2][2] * (b)[2][2]; }
352 for(iter=0; iter<
ITER_SJ; iter++){
355 if( A.
m[i][i] == 0 )
continue;
360 if(i != j) x_1.
v[i] += A.
m[i][j] * x.
v[j];
362 x_1.
v[i] = (b.
v[i] - x_1.
v[i]) / A.
m[i][i];
387 double detA =
det(A);
434 printf(
"%s [-n NEWTON METHOD ITERATIONS] [-j JACOBI METHOD ITERATIONS] FILE \n",name);
436 printf(
" -n NEWTON METHOD ITERATIONS number of iterations (default 10)\n");
437 printf(
" -j JACOBI METHOD ITERATIONS number of iterations (default 100)\n");
438 printf(
" FILE input graph file\n");
460 double number = (double) random();
468 if ((random() % 100) <
P_FIXED){
495 while ((c = (
char) getopt (argc, argv,
"hn:j:")) != -1)
498 sscanf(optarg,
"%d",&
iter1);
501 sscanf(optarg,
"%d",&
iter2);
512 if(optind == (argc-1)){
515 char * graph_file = argv[optind];
517 printf(
"# Graph: %s\n",graph_file);
538 printf(
"# N fixed: %d\n",nfixed);
582 Vector v4 = {{0,100,200}};
584 Vector v5 = {{100,100,0}};
586 Vector v6 = {{200,100,0}};
588 Vector v7 = {{300,100,0}};
624 norm += pow(Gi.
v[d],2);
638 int main(
int argc,
char ** argv) {
657 for(i=0;i<
iter1;i++){
658 printf(
"# iter: %d/%d\n",i,iter1);
673 for(j=0;j<
iter2;j++){
741 double scalar = -
K * ((
L - n) / n);
767 double scalar_a, scalar_b;
780 scalar_a = -((
L-n)/(n));
782 scalar_b = (
L/(pow(n,3)));
813 Vector gij =
g(vertex,nghb_index);
817 Matrix Wij =
W(vertex,nghb_index);
824 Matrix Hij =
W(nghb_index,vertex);
857 if(nghb_index == vertex)
continue;
881 new_ei =
solve(Hii,sum);
void print_help(char *name)
#define hit_tileNewType(baseType)
#define hit_bShapeAddEdge2(shape, x, y)
#define hit_fileHBReadBitmap(hbfile)
#define hit_bShapeEdgeTarget(s, edge)
#define hit_bShapeEdgeIterator(var, shape, vertex)
void hit_gbTileCopyVertices(void *destP, void *srcP)
void hit_gbTileClearVertices(void *varP)
void solve_system_iter(int vertex)
#define hit_gbTileVertexAt(var, vertex)
Vector solve(Matrix A, Vector b)
#define hit_bShapeVertexIterator(var, shape)
Vector solve_jacobi(Matrix A, Vector b)
void hit_comInit(int *pargc, char **pargv[])
#define hit_bShapeNvertices(shape)
void random_coordinates()
#define hit_gbTileEdgeAt(var, pos1, pos2)
int main(int argc, char *argv[])
HitShape HIT_BITMAP_SHAPE_NULL
void init_graph(HitTile_double graph, HitShape shape_global)
void calculate_GH(int vertex)
#define multMV(r, mm, vv)
#define hit_tileFree(var)
void hit_shapeFree(HitShape s)
#define hit_clockGetSeconds(c)
#define hit_gbTileDomainShapeAlloc(var, baseType, shape, allocOpts)
#define hit_gbTileEdgeIteratorAt(var, vertex, edge_index)
#define hit_clockStart(c)