https://github.com/janverschelde/PHCpack
Tip revision: 15b934b91f644c5c82f821311af3396c048880b0 authored by Jan Verschelde on 25 September 2014, 17:22:46 UTC
release 2.3.94 exports the new option phc -o for the contents of the symbol table after reading an input system
release 2.3.94 exports the new option phc -o for the contents of the symbol table after reading an input system
Tip revision: 15b934b
parallel_pieri.c
#include<stdio.h>
#include<stdlib.h>
#include "mpi.h"
#include "parallel_tree.h"
#include "queue.h"
#define SEND_PIVOT 100 /* tag for sending pivots */
#define SEND_SSOL 101 /* tag for sending start solutions */
#define SEND_PIVOT_B 102 /* tag for sending pivot back to the root node */
#define SEND_TSOL 103 /* tag for sending target solutions */
#define SEND_RES 104 /* tag for sending residuals */
#define v 0 /* verbose flag */
extern void adainit();
extern int _ada_use_c2pieri ( int job, int *a, int *b, double *c );
extern void adafinal();
void dimension_broadcast ( int myid, int *m, int *p, int *q, int *nd_tp );
/* Root node reads the dimension information and the pivot method to use, then
broadcasts to all the nodes */
void input_planes_broadcast( int myid, int m, int p, int q );
/* Root node initializes the machine with m*p + q*(m+p) random input m-planes,
* then broadcasts to all the nodes */
void interpolation_points_broadcast(int myid, int m, int p, int q);
/* Root node initializes the machine with m*p + q*(m+p) random interpolation
points, then broadcasts to all the nodes */
int sever_distribute ( int m, int p, int q, int numprocs, int nd_tp );
/* Server generates jobs with a tree structure, then distributes jobs to other
nodes with a queue. Returns the total number of solutions */
void client_compute ( int m, int p, int q, int nd_tp );
/* Clients do the computation jobs received from the server with Ada program */
int main(int argc, char* argv[])
{
int numprocs,myid,fail=-1;
int m,p,q,sum,nd_tp;
double startwtime, endwtime;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Comm_rank(MPI_COMM_WORLD,&myid);
srand(time(NULL));
adainit();
dimension_broadcast(myid,&m,&p,&q,&nd_tp);
input_planes_broadcast(myid,m,p,q);
interpolation_points_broadcast(myid,m,p,q);
if(myid==0)
{
startwtime = MPI_Wtime();
sum = server_distribute(m,p,q,numprocs,nd_tp);
endwtime = MPI_Wtime();
printf("\nTotal wall time = %lf seconds on %d processors\n",
endwtime-startwtime, numprocs);
}
else
client_compute(m,p,q,nd_tp);
adafinal();
MPI_Finalize();
return 0;
}
void dimension_broadcast ( int myid, int *m, int *p, int *q, int *nd_tp )
{
int fail=-1,mpq[3];
int *b;
double *c;
if(myid==0)
{
printf("Give p, dimension of the solution planes :");
scanf("%d", p);
printf("Give m, the co-dimension so that n = m+p :");
scanf("%d", m);
printf("Give q, the degree of the maps :");
scanf("%d", q);
while(1)
{
printf("Menu of pivot methods:\n");
printf(" 1. top pivots.\n");
printf(" 2. bottom pivots.\n");
printf(" 3. mixed pivots.\n");
printf("Type 1, 2, 3 to choose: ");
scanf("%d", nd_tp);
if(*nd_tp<1 || *nd_tp>3)
printf("Please choose 1~3.\n");
else
break;
}
}
MPI_Bcast(p, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(m, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(q, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(nd_tp, 1, MPI_INT, 0, MPI_COMM_WORLD);
mpq[0] = *m;
mpq[1] = *p;
mpq[2] = *q;
fail = _ada_use_c2pieri(1,mpq,b,c);
if(fail!=0)
{
if(myid==0)
printf("Invalid choice. Please quit the program and try again...\n");
}
if(v==1)
printf("No. %d processor:p=%d,m=%d,q=%d\n",myid,*p,*m,*q);
}
void input_planes_broadcast( int myid, int m, int p, int q )
{
int n = m*p + q*(m+p);
int i,nbcff = 2*(m+p)*m*n;
int fail,mpq[3];
double c[nbcff];
mpq[0] = m; mpq[1] = p; mpq[2] = q;
if(myid==0)
{
for (i=0; i<nbcff;i++) c[i] = ((double) rand())/RAND_MAX;
}
MPI_Bcast(c, nbcff, MPI_DOUBLE, 0, MPI_COMM_WORLD);
fail = _ada_use_c2pieri(2,mpq,&n,c);
}
void interpolation_points_broadcast(myid,m,p,q)
{
int n = m*p + q*(m+p);
int i,nbcff = 2*(m+p)*m*n;
int fail,mpq[3];
double c[nbcff];
mpq[0] = m; mpq[1] = p; mpq[2] = q;
nbcff = 2*n;
if(myid==0)
for (i=0; i<nbcff;i++) c[i] = ((double) rand())/RAND_MAX;
MPI_Bcast(c, nbcff, MPI_DOUBLE, 0, MPI_COMM_WORLD);
fail = _ada_use_c2pieri(3,mpq,&n,c);
}
void package(int p,int *start,int *target,int length,int address,int pack[4*p+2])
/* Saves top and bottom pivots of the start solution in the first 2*p elements,
appends top and bottom pivots of the target solution in the second 2*p
elements,length of solution in the (4p+1)th position,address of the node in
the (4p+2)th position */
{
int i;
for(i=0; i<2*p; i++)
pack[i] = start[i];
for(i=0; i<2*p; i++)
pack[i+2*p] = target[i];
pack[4*p] = length;
pack[4*p+1] = address;
}
void print_solution(int length, double *sol)
/* Prints out the solutions */
{
int i;
for(i=0; i<length; i=i+2)
{
if(sol[i]>0) printf(" ");
printf("%.15e\t",sol[i]);
if(sol[i+1]>0) printf(" ");
printf("%.15e\n",sol[i+1]);
}
}
int server_distribute ( int m, int p, int q, int numprocs, int nd_tp )
{
int length,buffer[4*p+2],counter,sol_num=0,/* the total number of solutions*/
dest=1,source,client_number = numprocs-1,finish=-1,diff,flag=1;
int i,plength=4*p+2,pack[plength],dim = m*p+q*(m+p),wait_cpu[numprocs-1];
Node *leave,*nd,*top_leave,*bottom_leave;
double *sol=NULL,res;
queue_node *nq_front=NULL, *nq_end=NULL;
/* nq_front and nq_end point the front and the end of the queue for nodes */
queue *front=NULL, *end=NULL;
/* front and end points the front and the end of the queue for jobs */
job *j;
MPI_Status status;
if(nd_tp==1)
{
nd=Top_Root(m,p,q);
leave=Top_Leave(m,p,q);
}
else if(nd_tp==2)
{
nd=Bottom_Root(p);
leave=Bottom_Leave(m,p,q);
}
else
{
leave=Mixed_Leave(m,p,q);
top_leave=Top_Leave(m,p,q);
bottom_leave=Bottom_Leave(m,p,q);
Mixed_Root(m,p,q,leave,&nq_front,&nq_end);
}
for(i=0; i<numprocs-1; i++)
wait_cpu[i] = i+1;
/* At the beginning all the client processors are waiting for jobs */
while(1) /* distributes computation tasks to slave processors */
{
if(nd_tp!=3)
{
if(nd->level==dim)
/* if nd is a leave,prints out the solution and put the
client ID in an array */
{
sol_num++;
printf("Soltion No. %d\n",sol_num);
print_solution(dim*2,sol);
MPI_Recv(&res,1,MPI_DOUBLE,MPI_ANY_SOURCE,SEND_RES,
MPI_COMM_WORLD,&status);
printf("The residual is : %.15e\n\n",res);
Deallocate(nd,p);
/*deallocate memory for the leave*/
if(v==1)
printf("client_number=%d\n",client_number);
if(client_number==(numprocs-1)) /* All of the clients returned leaves */
break;
}
else
{
if(nd_tp==1)
{
if(v==1)
{
printf("Create top parent for node ");
Print_Localization(p,nd->top);
}
Q_Create_Top_Parent(nd,m,p,q,leave,sol,&front,&end);
}
else if(nd_tp==2)
{
if(v==1)
{
printf("Create bottom parent for node ");
Print_Localization(p,nd->bottom);
}
Q_Create_Bottom_Parent(nd,m,p,q,leave,sol,&front,&end);
}
}
}
else /* the node type is mixed */
{
if((flag==1)&&(leave->level%2 == 1))
{
diff=1;
flag=0;
}
else
diff=2;
/*if the level of the leave is odd, add one degree of freedom
on the first level to save computation, add two degree of
freedom for the other levels.*/
if(nq_front!=NULL) /*generates jobs originating from the forest roots */
{
while(nq_front!=NULL)
/* finds all the parents for the node in the queue */
{
nq_front=pop_node(nq_front,&nd);
if(v==1)
{
printf("Create mixed parent for node: \n");
Print_Localization(p,nd->top);
Print_Localization(p,nd->bottom);
}
Q_Create_Mixed_Parent(nd,m,p,q,top_leave,bottom_leave,diff,
sol,&front,&end);
}
}
else
{
if(nd->level==dim)
/* if nd is a leave,prints out the solution and put the
client ID in an array */
{
sol_num++;
printf("Solution No. %d\n",sol_num);
print_solution(dim*2,sol);
MPI_Recv(&res,1,MPI_DOUBLE,MPI_ANY_SOURCE,SEND_RES,
MPI_COMM_WORLD,&status);
printf("The residual is : %.15e\n\n",res);
Deallocate(nd,p);
/*deallocate memory for the leave*/
if(v==1)
printf("client_number=%d\n",client_number);
if(client_number==(numprocs-1))
/* All of the clients returned leaves */
break;
}
else
{
if(v==1)
{
printf("Create mixed parent for node: \n");
Print_Localization(p,nd->top);
Print_Localization(p,nd->bottom);
}
Q_Create_Mixed_Parent(nd,m,p,q,top_leave,bottom_leave,diff,
sol,&front,&end);
}
}
}
if(sol!=NULL)
free(sol);
while((front!=NULL)&&(client_number>0))
/* Distributes one job to each of the client processors at the beginning */
{
front = pop(front, &j);
package(p,j->start_pivot,j->target_pivot,j->length,j->address,pack);
MPI_Send(pack,plength,MPI_INT,wait_cpu[client_number-1],
SEND_PIVOT,MPI_COMM_WORLD);
MPI_Send(j->start_sol,j->length,MPI_DOUBLE,wait_cpu[client_number-1],
SEND_SSOL,MPI_COMM_WORLD);
if(v==1)
{
printf("Sending a job to No. %d processor: ",wait_cpu[client_number-1]);
Print_Localization(plength, pack);
printf("The start solution:\n");
print_solution(j->length,j->start_sol);
}
Deallocate_job(j);
client_number--;
}
MPI_Recv(buffer,plength,MPI_INT,MPI_ANY_SOURCE,SEND_PIVOT_B,
MPI_COMM_WORLD,&status);
length = buffer[4*p];
sol = (double*)calloc(length, sizeof(double));
source = status.MPI_SOURCE;
MPI_Recv(sol,length,MPI_DOUBLE,source,SEND_TSOL,
MPI_COMM_WORLD,&status);
wait_cpu[client_number++] = source;
if(v==1)
{
printf("Received pivot array from No. %d processor: ", source);
Print_Localization(4*p+2, buffer);
printf("Received solution from No. %d processor:\n", source);
for(i=0; i<length; i=i+2)
{
if(sol[i]>0) printf(" ");
printf("%.14e\t",sol[i]);
if(sol[i+1]>0) printf(" ");
printf("%.14e\n",sol[i+1]);
}
}
nd = (Node*) buffer[4*p+1]; /* Assigns the node address to nd */
}
/* Tell all the clients no more jobs to be done */
for(i=1; i<numprocs; i++)
MPI_Send(&finish,1,MPI_INT,i,SEND_PIVOT,MPI_COMM_WORLD);
Deallocate(leave,p);
if(nd_tp==3)
{
Deallocate(top_leave,p);
Deallocate(bottom_leave,p);
}
printf("The number of roots : %d\n", sol_num);
return sol_num;
}
void client_compute ( int m, int p, int q, int nd_tp )
{
int length,plength=4*p+2,buffer[plength],diff,
i,dim=m*p+q*(m+p),fail,*a,*b,deg_freedom;
double ssol[dim*2],tsol[dim*2],*c,res;
MPI_Status status;
while(1)
{
MPI_Recv(buffer,plength,MPI_INT,0,SEND_PIVOT,MPI_COMM_WORLD,&status);
if(buffer[0]==-1) break;
length = buffer[4*p];
MPI_Recv(ssol,length,MPI_DOUBLE,0,SEND_SSOL,MPI_COMM_WORLD,&status);
fail = _ada_use_c2pieri(4,&p,buffer,c);
/* store pivots of pattern at start solution curve */
fail = _ada_use_c2pieri(5,&p,&buffer[2*p],c);
/* store pivots of pattern at target solution curve */
deg_freedom = length/2;
fail = _ada_use_c2pieri(6,°_freedom,b,ssol);
/* store coefficients of start solution curve */
if(v==1)
{
printf("The start curve is:\n");
for(i=0;i<2*deg_freedom;i=i+2)
{
if(ssol[i]>0) printf(" ");
printf("%.14e\t",ssol[i]);
if(ssol[i+1]>0) printf(" ");
printf("%.14e\n",ssol[i+1]);
}
}
fail = _ada_use_c2pieri(8,a,b,tsol);
/* track solution path without intermediate output */
if(nd_tp==3)
{
if((buffer[4*p]==0)&&((dim%2)==1))
diff = 1;
else
diff = 2;
}
else
diff = 1;
deg_freedom+=diff;
fail = _ada_use_c2pieri(7,°_freedom,b,tsol);
/* retrieve coefficients of target solution curve */
buffer[4*p] = buffer[4*p]+2*diff;
MPI_Send(buffer,plength,MPI_INT,0,SEND_PIVOT_B,MPI_COMM_WORLD);
MPI_Send(tsol,length+2*diff,MPI_DOUBLE,0,SEND_TSOL,MPI_COMM_WORLD);
if(deg_freedom==dim)
{
fail = _ada_use_c2pieri(10,a,b,&res);
/* verify intersection conditions without output */
MPI_Send(&res,1,MPI_DOUBLE,0,SEND_RES,MPI_COMM_WORLD);
}
}
}