https://github.com/janverschelde/PHCpack
Tip revision: fa6da5f66dd0a2c5843939c9118759669a8117fe authored by janv@uic.edu on 31 July 2015, 21:14:37 UTC
added description of a session with phc -s and updated conf.py for version 2.3.99
added description of a session with phc -s and updated conf.py for version 2.3.99
Tip revision: fa6da5f
mpi2fac_d.c
/* The program "mpi2fac_d" is a parallel implementation of the factorization
* of a pure positive dimensional solution set into irreducible components,
* using a dynamic load distribution.
* The program prompts the user for an embedded system (system, slices,
* and witness set) and for the maximal number of monodromy loops. */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "mpi.h"
#include "../Lib/phcpack.h"
#include "../Lib/syscon.h"
#include "../Lib/solcon.h"
#include "parallel_phcpack.h"
#include "parallel_monodromy.h"
#include "manage_components.h"
#define v 1 /* verbose flag:
0 no output during the computations,
1 only one line messages on progress of computations
2 display of data, like systems and solutions */
#define tf 0 /* testing flag:
0 no interactive test is necessary;
1 some interactive tests will be run */
int initialize_slices ( int myid, int n, int dim, int deg, int nb_slices,
int numprocs, double *trackwtime );
/*
* DESCRIPTION :
* Computes as many witness sets as the value of nb_slices + 2,
* to decompose a solution set of degree deg and dimension dim
* in n-space into irreducible factors.
*
* ON ENTRY :
* myid identification of the node;
* n ambient dimension;
* dim dimension of the solution set;
* deg degree of the solution set;
* nb_slices maximal number of slices to be used;
* numprocs number of processors.
*
* ON RETURN :
* trackwtime time node myid has been tracking paths;
* the return value of the function is either 0 or 1:
* 0 if the trace grid and all slices are well,
* 1 otherwise: if the trace grid is not good. */
void start_walk(int sliceA, int start_label, int sliceB, SLAVE_STATE* state);
/*
* DESCRIPTION :
* This is the main job scheduling function. */
void run_master(int n, int s, int max_s, int np, double* breakup_time);
/*
* DESCRIPTION :
* This function is executed by the root node 0. */
void run_slave(int n, double* trackwtime);
/*
* DESCRIPTION :
* This function is executed by the computational nodes. */
void pick_slices_random(int s, int* i, int* j);
/*
* DESCRIPTION :
* Returns two random numbers "i","j" from 0 to "s"-1. */
void pick_slices_by_rejection_rate(int s, int max_s, ISLICE_STATS* stats, int* i, int* j);
/*
* DESCRIPTION :
* Picks a pair of slices with the total lowest rejection rate. */
void pick_slices_by_truncated_rejection_history(int s, int max_s, ISLICE_STATS* stats,
int* i, int* j);
/*
* DESCRIPTION :
* Picks a pair of slices with the lowest rejection rate in the short past. */
int main ( int argc, char *argv[] )
{
int myid,numprocs,n,dim,deg,nb_slices,fail;
double startwtime,trackwtime,wtime,*mytime,init_time,breakup_time;
MPI_Status status;
adainit();
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
srand(time(NULL)); /* for different random gamma constant */
if(myid == 0)
{
mytime = (double*) calloc(numprocs, sizeof(double));
fail = read_witness_set(&n,&dim,°);
fail = define_output_file();
printf("\nGive the number of slices : "); scanf("%d",&nb_slices);
startwtime = MPI_Wtime();
}
else
trackwtime = MPI_Wtime();
MPI_Barrier(MPI_COMM_WORLD); /* everybody must wait for node 0 */
dimension_broadcast(myid,&n);
MPI_Bcast(&nb_slices,1,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&dim,1,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(°,1,MPI_INT,0,MPI_COMM_WORLD);
/* initialize slices */
init_time = MPI_Wtime();
fail = initialize_slices(myid,n,dim,deg,nb_slices,numprocs,&trackwtime);
init_time = MPI_Wtime()-init_time;
MPI_Barrier(MPI_COMM_WORLD);
if(myid == 0)
run_master(deg,nb_slices+1,nb_slices+1,numprocs, &breakup_time);
else {
trackwtime = 0; //counting slaves time only after this point
run_slave(deg,&trackwtime);
}
MPI_Barrier(MPI_COMM_WORLD);
if(myid == 0)
wtime = MPI_Wtime() - startwtime;
else
wtime = trackwtime;
MPI_Gather(&wtime,1,MPI_DOUBLE,mytime,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
if(myid == 0) {
printf("Initialization time = %f seconds.\n", init_time);
printf("Breakup time = %f seconds.\n", breakup_time);
print_time(mytime,numprocs);
}
MPI_Finalize();
adafinal();
return 0;
}
int initialize_slices ( int myid, int n, int dim, int deg, int nb_slices,
int numprocs, double *trackwtime )
{
int fail;
fail = initialize_monodromy(nb_slices,deg,dim);
fail = witness_set_distribute(myid,n,dim,deg,numprocs);
MPI_Barrier(MPI_COMM_WORLD); /* initialization of sampler done */
all_slices_broadcast(myid,nb_slices,dim,n);
fail = build_trace_grid(myid,n,dim,deg,numprocs,trackwtime);
MPI_Barrier(MPI_COMM_WORLD); /* wait for finish build_trace_grid */
trace_grid_broadcast(myid,deg,n);
MPI_Barrier(MPI_COMM_WORLD); /* everybody waits for trace grid */
if(tf > 0)
{
if(myid == 0) printf("\nTesting linear trace operations ...\n");
interactive_trace_test(myid);
}
fail = track_paths_to_new_slices
(myid,nb_slices,deg,dim,n,numprocs,trackwtime);
MPI_Barrier(MPI_COMM_WORLD);
if(tf > 0)
{
if(myid == 0) printf("\nRunning the in_slice_test ...\n");
interactive_in_slice_test(myid);
}
MPI_Barrier(MPI_COMM_WORLD);
if(tf > 0)
{
if(myid == 0) printf("\nTesting sample_loop interactively ...\n");
interactive_track_to_make_loops(myid);
}
return fail;
}
void start_walk (int sliceA, int start_label, int sliceB, SLAVE_STATE* state)
{
/* change state */
state->status = TRACKING;
//(is_in_slice(start_label, sliceB)) ? TRACKING : NEW_POINT;
state->sliceA = sliceA;
state->start_label = start_label;
state->sliceB = sliceB;
/* send the job to slave */
if (v>=3) printf("M: sending a job to slave %d\n", state->num);
MPI_Send(state, sizeof(SLAVE_STATE), MPI_CHAR,
state->num, /* the number of the slave */
SEND_TRACK_JOB, MPI_COMM_WORLD);
}
void run_master(int n, int s, int max_s, int np, double* breakup_time)
{
int i,j,k,sl,t,c,c_nodes,pair_index;
double endwtime, startwtime;
POINT pt;
int n_paths = 0;
int n_rejected = 0;
ISLICE_STATS* stats = (ISLICE_STATS*)calloc(max_s*max_s,sizeof(ISLICE_STATS));
SLAVE_STATE* state;
int* ind;
/* ind[i] = the number of the component that i-th point belongs to */
int n_chunks;
COMPONENT* part;
MPI_Status status;
int flag, ept, spt;
startwtime = MPI_Wtime();
if (v>1) printf("runmaster: n=%d, s=%d, max_s=%d, np=%d\n", n,s,max_s,np);
for (i=0; i<max_s; i++) for (j=i+1; j<max_s; j++)
stats[i*max_s+j].n_rejected = stats[i*max_s+j].n_paths = 0;
/* records what slave is doing */
state = (SLAVE_STATE*) calloc(np, sizeof(SLAVE_STATE));
for (sl=1; sl<np; sl++) {
state[sl].num = sl;
state[sl].status = FREE;
}
ind = (int*)calloc(n+1,sizeof(int));
n_chunks = n;
part = (COMPONENT*)calloc(n,sizeof(COMPONENT));
/* partition: list of lists */
for(i=0; i<n; i++){
ind[i+1] = i; /* each point is a singleton in the beginning,
point labels start from 1 */
makeCOMPONENT(&part[i], i, 1);
part[i].pts[0] = i+1;
if (is_irreducible(&part[i])) {
part[i].is_irred = 1;
n_chunks--;
}
part[i].n_nodes_tracking = 0;
}
/* MAIN LOOP */
#define MAX_PATHS 10000 /*!!!*/
while (n_chunks>0
&& n_paths<MAX_PATHS) {
/* make busy the slaves */
while(free_slave(state, np)>0) {
/* pick a pair of slices */
//pick_slices_random(s,&i,&j);
//pick_slices_by_rejection_rate(s, max_s, stats,&i,&j);
pick_slices_by_truncated_rejection_history(s, max_s, stats,&i,&j);
/* find the smallest reducible component
with the smallest number of nodes assigned to it
(!!! we should implement sorting instead!!!)*/
c = -1; c_nodes = 0;
for (t=1; t<n; t++){
if (part[t].size!=0
&& !part[t].is_irred
&& (c==-1
|| part[t].size<part[c].size
|| (part[t].size==part[c].size
&& part[c].n_nodes_tracking < c_nodes)
))
{
c_nodes = part[c].n_nodes_tracking;
c = t;
}
}
start_walk(i, part[c].pts[0], j, &state[free_slave(state, np)]);
part[c].n_nodes_tracking++; //one more slave is assigned to this component
n_paths++;
}
/* collect results */
MPI_Iprobe(MPI_ANY_SOURCE, SEND_END_POINT, MPI_COMM_WORLD, &flag, &status);
while (flag) {
sl = status.MPI_SOURCE; /* the number of slave who sent this */
if (v>=3) printf("M: received a message from slave %d\n", sl);
MPI_Recv(&ept, /* end point */
1, MPI_INT, sl,
SEND_END_POINT, MPI_COMM_WORLD, &status);
spt = state[sl].start_label; /* start point */
part[ind[spt]].n_nodes_tracking--; // one fewer slave is assigned to this component
pair_index = state[sl].sliceA*max_s+state[sl].sliceB; /* index for "stats"[][] */
if (v>=2) printf("M: there is a path from %d to %d\n", spt, ept);
if (ept>0 && ept<=n
&& ind[ept] != ind[spt] /* if not in the same component... */
&& !part[ind[ept]].is_irred)
{
record_path_outcome(&stats[pair_index], PATH_LED_TO_MERGER);
mergeCOMPONENT(&part[ind[ept]],&part[ind[spt]],ind);
n_chunks--;
if (is_irreducible(&part[ind[ept]])) {
part[ind[ept]].is_irred = 1;
n_chunks--;
}
if (v>=3) {
printf("#chunks = %d, index: ", n_chunks);
for (i=1; i<=n; i++)
printf("[%d]", ind[i]);
printf("\n");
}
} else {
n_rejected++;
if (ept<=0 || ept>n) t=PATH_DIVERGED;
else if (ind[ept] != ind[spt] && part[ind[ept]].is_irred) t=PATH_CROSSED;
else t=PATH_STAYED_WITHIN_COMPONENT;
record_path_outcome(&stats[pair_index],t);
}
if (v>=2) print_partition(n, part); /* print the updated partition */
state[sl].status = FREE;
MPI_Iprobe(MPI_ANY_SOURCE, SEND_END_POINT, MPI_COMM_WORLD, &flag, &status);
}// end "collect results"
}
/* TERMINATE SLAVES */
for (sl=1; sl<np; sl++)
MPI_Send(&sl, 1, MPI_INT,
sl, /* the number of the slave */
SEND_QUIT, MPI_COMM_WORLD);
/* PRINT STATISTICS */
printf(
"*** Statistics ***\n\
Summary: \n\
paths tracked\n\
%d -- for trace grid,\n\
%d -- for slices initialization;\n\
after that \n\
%d more paths tracked,\n\
%d rejected.\n", 2*n, n*(s-1), n_paths, n_rejected);
printf("For pairs of slices:\n");
for (i=0; i<s; i++) for (j=i+1; j<s; j++){
int total;
printf(" (%d)<->(%d): \n\
history [", i,j);
total = MIN(stats[i*max_s+j].n_paths, N_PATHS_TO_REMEMBER);
for(k=0; k<total; k++){
switch(stats[i*max_s+j].path_outcome[k]) {
case PATH_LED_TO_MERGER: printf("M"); break;
case PATH_DIVERGED: printf("d"); break;
case PATH_STAYED_WITHIN_COMPONENT: printf("c"); break;
case PATH_CROSSED: printf("x"); break;
};
}
printf("]\n\
%d paths tracked,\n\
%d rejected.\n",
stats[i*max_s+j].n_paths, stats[i*max_s+j].n_rejected);
}
print_partition(n,part);
/* CLEAN-UP */
free(ind);
for(i=0; i<n; i++){
killCOMPONENT(&part[i]);
}
free(part);
free(stats);
endwtime = MPI_Wtime();
*breakup_time += (endwtime - startwtime);
}//end{run_master}
void run_slave(int n, double* trackwtime)
{
MPI_Status status;
SLAVE_STATE st;
int target_label,fail,flag = 0;
int sl;
double endwtime, startwtime;
while(!flag) /* while there is no "quit" message */
{
MPI_Iprobe(MPI_ANY_SOURCE, SEND_TRACK_JOB, MPI_COMM_WORLD, &flag, &status);
if (flag) {
MPI_Recv(&st, sizeof(SLAVE_STATE), MPI_CHAR, status.MPI_SOURCE,
SEND_TRACK_JOB, MPI_COMM_WORLD, &status);
if (v>=3) printf("S_%d: tracking label %d from slice %d to slice %d\n",
st.num, st.start_label, st.sliceA, st.sliceB);
startwtime = MPI_Wtime();
fail = sample_loop(st.sliceA,st.sliceB,st.start_label,&target_label);
endwtime = MPI_Wtime();
*trackwtime += (endwtime - startwtime);
// end_point = (rand() % n) + 1; // simulate the discovery of the endpoint
if (v>=3) printf("S_%d: done tracking, returning label %d\n", st.num, target_label);
MPI_Send(&target_label, 1, MPI_INT,
0, /*status.MPI_SOURCE, /* send to where the job was received from */
SEND_END_POINT, MPI_COMM_WORLD);
}
MPI_Iprobe(MPI_ANY_SOURCE, SEND_QUIT, MPI_COMM_WORLD, &flag, &status);
}
/* slave receives its number and quits */
MPI_Recv(&sl, 1, MPI_INT, status.MPI_SOURCE,
SEND_QUIT, MPI_COMM_WORLD, &status);
if (v>=3) printf("S_%d: finita la comedia\n", sl);
}
void pick_slices_random(int s, int* i, int* j)
{
int a,b;
a = rand()%s;
b = (a + rand()%(s-1)+1)%s;
if (a<b) {*i=a; *j=b;} else {*i=b; *j=a;};
}
void pick_slices_by_rejection_rate(int s, int max_s, ISLICE_STATS* stats, int* i, int* j)
{
int ii,jj; double rr = 2, tr;
while (rr>1) {
for (ii=0; ii<s; ii++)
for (jj=ii+1; jj<s; jj++)
if (isGood(&stats[ii*max_s+jj])
// && ii>0 // avoid the base slice for now !!!
) {
tr = (stats[ii*max_s+jj].n_paths==0) ? 0:
(double) stats[ii*max_s+jj].n_rejected / stats[ii*max_s+jj].n_paths;
if (tr<rr)
{ rr = tr; *i = ii; *j = jj; }
}
};
}
void pick_slices_by_truncated_rejection_history(int s, int max_s, ISLICE_STATS* stats,
int* i, int* j)
{
int ii,jj,t,total,c;
double rate,top = -1;
for (ii=0; ii<s; ii++)
for (jj=ii+1; jj<s; jj++)
{
c = 0;
total = MIN(stats[ii*max_s+jj].n_paths, N_PATHS_TO_REMEMBER);
for (t=0; t<total; t++)
if (stats[ii*max_s+jj].path_outcome[t] == PATH_LED_TO_MERGER)
c++; // count the number of successful paths
rate = (2*total<N_PATHS_TO_REMEMBER) ? 1 : (double)c/total;
/* there should be at least a half of history record,
before a heuristic choice is made */
if (rate>top)
{ top = rate; *i = ii; *j = jj; }
}
}