The text of the mpC N-body program presented in Sections 7.1 and 9.2 is broken down into two source files. The first file contains mpC-specific code of the application. The name of this file must have extension mpc (e.g., nbody.mpc). The contents of this file are as follows:

#include <stdio.h> #include <stdlib.h> #include <mpc.h> typedef double Triplet[3]; typedef struct{Triplet p; Triplet v; double m;} Body; nettype Nbody(int m, int k, int n[m]) {    coord I=m;    node { I>=0: bench*((n[I]/k)*(n[I]/k));};    link { I>0: length*(n[I]*sizeof(Body)) [I]->[0];};    parent [0];    scheme {       int i;       par (i=0; i<m; i++) 100%%[i];       par (i=1; i<m; i++) 100%%[i]->[0];    }; }; repl M;         // The number of groups double [host]t;   // The time variable repl int tgsize; // The number of bodies in a test group #include "serial_code.h" void [net Nbody(m, l, n[m]) g] ShareMasses (void *masses) {    double mass;    repl i, j;    typedef double (*pArray)[m];    mass = (*(pArray)masses)[I coordof i];    for(i=0; i<m; i++)       for(j=0; j<m; j++)          [g:I==i](*(pArray)masses)[j] = [g:I==j]mass; } void [net Nbody(m, l, n[m]) g] ShareCenters(void *centers) {    Triplet center;    repl i, j;    typedef Triplet (*pArray)[m];    center[] = (*(pArray)centers)[I coordof i][];    for(i=0; i<m; i++)       for(j=0; j<m; j++)          [g:I==i](*(pArray)centers)[j][] =           [g:I==j]center[]; } void [*]main(int [host]argc, char **[host]argv) {    // Get the number of groups on the host-process    // and broadcast it to all processes    M = [host]GetNumberOfGroups(argv[1]);    {       repl N[M]; // Array of group sizes       Body *[host]Groups[[host]M];       // Initialize N and Groups on the abstract        host-processor       [host]: InputGroups(argv[1], M, &N, &Groups);       // Broadcast the group sizes across all processes       N[] = [host]N[];       // Set the size of test group to the size of the smallest       // group in the system of bodies       tgsize = [?<]N[];       {          Body OldTestGroup[tgsize], TestGroup[tgsize];          typedef Body (*pTestGroup)[tgsize];          TestGroup[] = (*(pTestGroup)Groups[0])[];          OldTestGroup[] = TestGroup[];          recon UpdateGroup(tgsize, &OldTestGroup,                    &TestGroup, 1, NULL, NULL, 0); } {    net Nbody(M, tgsize, N) g;    int [g]myN, // Size of my group          [g]mycoord; // My coordinate in network g    mycoord = I coordof g;    myN = ([g]N)[mycoord];    {       double [g]Masses[[g]M]; // Masses of all groups       Triplet [g]Centers[[g]M]; // Centers of mass of all groups       Body [g]OldGroup[myN], // State of my group at time t             [g]Group[myN]; // State of my group at time              t+DELTA       repl   [g]gc, // Counter of groups in the system of bodies             [g]bc, // Counter of bodies in group             [g]gsize; // Size of group       //Scatter groups       for(gc=0; gc<[g]M; gc++)       {          [host]: gsize = N[gc];          {             typedef Body (*pGroup)[[host]gsize];             [g:I==gc]Group[] = (*(pGroup)Groups[gc])[];       }    } // Visualize the system of bodies on the computer display // associated with the abstract host-processor [host]: FirstDrawGroups(argc, argv, M, N, Groups); // Compute masses of the groups in parallel for(bc=0, Masses[mycoord]=0.0; bc<myN; bc++)    Masses[mycoord] += Group[bc].m; // Communicate to share the masses among abstract // processors of network g ([([g]M, [g]tgsize, [g]N)g])ShareMasses(Masses); // Main loop of the parallel algorithm do {    OldGroup[] = Group[];    // Compute centers of masses of the groups in parallel    Centers[] = 0.0;    for(bc=0; bc<myN; bc++)       Centers[mycoord][] +=       (OldGroup[bc].m/Masses[mycoord])*(OldGroup[bc].p)[]; // Communicate to share the centers among abstract // processors of network g ([([g]M, [g]tgsize, [g]N)g])ShareCenters(Centers); // Update the groups of bodies in parallel ([g]UpdateGroup)(myN, &OldGroup, &Group, [g]M,                   &Centers, &Masses, mycoord);             // Increment the time variable t             t += DELTA;             // Gather all groups to the abstract host-processor             for(gc=0; gc<[g]M; gc++)             {                gsize = [host]N[[host]gc];                {                   typedef Body (*pGroup)[gsize];                   (*(pGroup)Groups[gc])[] = [g:I==gc]Group[];             }          }             // Visualise the groups on the computer display             // associated with the abstract host-processor             if(DrawGroups([host]M, [host]N, Groups)<0)                MPC Exit(-1);             } while(1);             [host]: DrawBye(M, N, Groups);          }       }    } }

The name of the second source file must be source_code.h. This file contains serial routine code and is included in the program with the include directive. In particular, this file contains GUI code implementing visualization of the system of bodies. The GUI code is based on the Xlib library. The contents of this file are as follows:

#include <stdio.h> #include <math.h> #include <string.h> //The number of bodies in groups of default system of bodies #define N0 200 #define N1 100 #define N2 400 #define N3 200 #define N4 600 #define X 0 #define Y 1 #define Z 2 #define Vx 3 #define Vy 4 #define Vz 5 #define MASS 6 #define ALPHA 0.5 #define BETA 0.5 #define GAMMA 6.67e-11 #define POINTSIZE 1.e7 #define MASSUNIT 1.e23 #define XMAX 800. #define YMAX 800. #define ZMAX 800. #define BOXSIZE 120. #define EPS 0.001 #define CRITICAL_R 10. #define INTERVAL 1 #define MAXGROUPS 9 /*9*/ #define MAXBODIES 600 /*600*/ #define DELTA 360. double [host]wtime; void Merging(), SingleBody(); void UpdateGroup(int n, Body (*OldGroup)[n], Body (*Group)[n],                    int m, double (*Centers) [m][3],                    double (*Masses)[m], int mycoord) {    int j, k;    double sigma[3], gma[3], F[3], Fabs, r, r2, dr[3];    double dV[n][n][3];    for(j=0; j<n; j++)       if((*OldGroup)[j].m>0.)       {          sigma[]=0.;          for(k=0; k<n; k++)             if(j<k&&(*OldGroup)[k].m>0.)             {                Body jbody;                dr[] = ((*OldGroup)[k].p)[]-((*OldGroup)[j].p)[];                r2 = [+](dr[]*dr[]);                r = sqrt(r2);                if(r/POINTSIZE < sqrt(((*OldGroup)[k].m+                   (*OldGroup)[j].m)/MASSUNIT)+CRITICAL_R)                {                (jbody.p)[] = dr[];                (jbody.v)[] =-                   ((*OldGroup)[k].v)[]+((*OldGroup)[j].v)[];                jbody.m = (*OldGroup)[k].m;                SingleBody(&jbody);                dV[j][k][] = (jbody.v)[];                sigma[] += dV[j][k][]/DELTA;                }                else                {                Fabs = GAMMA*(*OldGroup)[k].m/r2;                sigma[] += Fabs*(dr[]/r);    } } else if(j>k&&(*OldGroup)[k].m>0.) {    dr[] = ((*OldGroup)[k].p)[]-((*OldGroup)[j].p)[];    r2 = [+](dr[]*dr[]);    r = sqrt(r2);    if(r/POINTSIZE<sqrt(((*OldGroup)[k].m+                     (*OldGroup)[j].m)/MASSUNIT)+CRITICAL_R)    sigma[] -= ((*OldGroup)[k].m/(*OldGroup)[j].m)                   *dV[k][j][]/DELTA;       else       {          Fabs = GAMMA*(*OldGroup)[k].m/r2;          sigma[] += Fabs*(dr[]/r);       }    }    gma[] = 0.;    for(k=0; k<m; k++)       if(k!=mycoord)       {          dr[]=(*Centers)[k][]-((*OldGroup)[j].p)[];          r2=[+](dr[]*dr[]);          Fabs=GAMMA*(*Masses)[k]/r2;          r=sqrt(r2);          gma[]+=Fabs*(dr[]/r);       }       sigma[] += gma[];       ((*Group)[j].v)[] += sigma[]*DELTA;       ((*Group)[j].p)[] += (sigma[]*DELTA/2.+(ALPHA*                   ((*Group)[j].v)[]+BETA*                   ((*OldGroup)[j].v)[]))*DELTA;       ((*Group)[j].p)[Z] -= BETA*sigma[Z]*DELTA*DELTA;    }    Merging(*Group, n); } void Merging(Body *rs, int n) {    int i, j;    double dr[3], r2, r;    for(i=0; i<n; i++)       if(rs[i].m>0.)          for(j=i+1; j<n; j++)             if(rs[j].m>0.)             {             dr[] = (rs[j].p)[]-(rs[i].p)[];             r2 = [+](dr[]*dr[]);             r = sqrt(r2);             if((int)(r/POINTSIZE)<=(int)(sqrt(rs[i].m/MASSUNIT)/2)+                (int)(sqrt(rs[j].m/MASSUNIT)/2))             {             double m;             m = rs[i].m+rs[j].m;             (rs[i].p)[] =  ((rs[j].p)[]*rs[j].m+ (rs[i].p)[]*rs[i].m)/m;             (rs[i].v)[] = ((rs[j].v)[]*rs[j].m+(rs[i].v)[]                *rs[i].m)/m;             rs[i].m = m;             rs[j].m = 0.;             (rs[j].p)[] = 0.;             (rs[j].v)[] = 0.;             }       } } int GetNumberOfGroups(char *fname) {    FILE *pf;    int i, M;    if(fname!=NULL)    {       pf = fopen(fname, "r");       if(pf==NULL)       {          printf("Can’t open file ‘%s’\n", fname);          return -1;    }    if((i = fscanf(pf, " Number of groups = %d",&M))<1 ||       M<=0 || M>MAXGROUPS)    {       fclose(pf);       if(i<1)          printf("Cannot read number of groups from file ‘%s’\n",                 fname);    else       printf("Number of groups = %d\n", M);       return -1;    }    fclose(pf);    }    else       return 5; } int InputGroups(char *fname, int m, int (*N)[m],                   Body *(*Groups)[m]) {    FILE *pf;    int i, j;    double x0, y0, z0;    double r, phi, xi;    #define PI2 6.2831854    if(fname!=NULL)    {       pf = fopen(fname, "r");       if(pf==NULL)       {          printf("Can’t open file ‘%s’\n", fname);          return -1;       }       fscanf(pf, " Number of groups = %*d", &m);       if((i = fscanf(pf, " Group sizes ="))<0)       {          fclose(pf);          printf("Cannot read groups’ sizes from file ‘%s’\n",                    fname);          return -1;    }       for(i=0; i<m; i++)          if((j = fscanf(pf, " %d", *N+i))<1 ||             (*N)[i]<=0 || (*N)[i]>MAXBODIES)          {          fclose(pf);          if(j<1)             printf("Cannot read size of group #%d from ‘%s’\n",                    i, fname);          else             printf("Size of group #%d = %d\n", i, N[i]);          return -1;       }       fclose(pf);    }    else    {       (*N)[0] = N0;       (*N)[1] = N1;       (*N)[2] = N2;       (*N)[3] = N3;       (*N)[4] = N4;    }    for(i=0; i<m; i++)    {       (*Groups)[i] = calloc((*N)[i], sizeof(Body));       if((*Groups)[i]==NULL)       {          printf("Cannot allocate storage for group #%d\n", i);          return -1;    } } for(i=0; i<m; i++) {    switch(i)          {          case 0: x0 = y0 = BOXSIZE/2.; break;          case 1: x0 = BOXSIZE*1.5; y0 = YMAX-BOXSIZE/2.;           break;          case 2: x0 = XMAX/2.; y0 = YMAX/2.; break;          case 3: x0 = XMAX-BOXSIZE/2.; y0 = BOXSIZE/2.;           break;          case 4: x0 = XMAX-BOXSIZE/2.; y0 = YMAX-          BOXSIZE/2.; break;          case 5: x0 = BOXSIZE*1.5; y0 = YMAX/2.; break;          case 6: x0 = XMAX-BOXSIZE/2.; y0 = YMAX/2.;           break;          case 7: x0 = XMAX/2.; y0 = YMAX-BOXSIZE/2.;           break;          case 8: x0 = XMAX/2.; y0 = BOXSIZE/2.; break;       }       for(j=0; j<(*N)[i]; j++)       {          r = (BOXSIZE/2./RAND_MAX)*rand();          phi = (PI2/RAND_MAX)*rand();          ((*Groups)[i][j].p)[X] = (x0+r*cos(phi))*POINTSIZE;          ((*Groups)[i][j].p)[Y] = (y0+r*sin(phi))*POINTSHZE;          ((*Groups)[i][j].p)[Z] = 0.;          ((*Groups)[i][j].v)[] = 0.;          (*Groups)[i][j].m = (j==(*N)[i]/2)?50.*MASSUNIT                   :0.5*MASSUNIT;       }    }    return 1; } void SingleBody(Body *body) {    double r2, r, a, dt=DELTA, t, k1[3], k2[3], vk1[3],     vk2[3],       ak[3], ar[3];    int i=1;    r2 = [+]((body->p)[]*(body->p)[]);    a = GAMMA*body->m/r2;    r = sqrt(r2);    ar[] = a*((body->p)[]/r);    vk1[] = (body->v)[]+ar[]*dt;    k1[] = (body->p)[]+vk1[]*dt/2.+(body->v)[]*dt; label:    vk2[] = (body->v)[];    k2[] = (body->p)[];    ak[] = ar[];    for(dt/=2., t=0.; t<DELTA; t+=dt)    {       i++;       vk2[] += ak[]*dt;       k2[] += vk2[]*dt/2.+(ALPHA*vk2[]+BETA*(vk2[]-       ak[]*dt))*dt;       r2 = [+](k2[]*k2[]);       r = sqrt(r2);       a = GAMMA*body->m/r2;       ak[] = a*(k2[]/r);    }    if(fabs(vk2[X]-vk1[X])<=EPS*fabs(vk2[X]) &&       fabs(vk2[Y]-vk1[Y])<=EPS*fabs(vk2[Y]) &&       fabs(vk2[Z]-vk1[Z])<=EPS*fabs(vk2[Z]) &&       fabs(k2[X]-k1[X])<=EPS*fabs(k2[X]) &&       fabs(k2[Y]-k1[Y])<=EPS*fabs(k2[Y]) &&       fabs(k2[Z]-k1[Z])<=EPS*fabs(k2[Z]))    {       (body->p)[] = k2[];       (body->v)[] = vk2[]-(body->v)[];       return;    }    else    {       vk1[] = vk2[];       k1[] = k2[];       goto label;    } } #pragma keywords ANSI #include <X11/Xlib.h> #include <X11/Xutil.h> #include <X11/Xos.h> #include <X11/bitmaps/icon> #pragma keywords SHORT #define WIDTH XMAX #define HEIGHT YMAX Display *[host]display; int [host]screen; Window [host]win; GC [host]gc; XEvent [host]report; XFontStruct *[host]font_info; char *[host]fontname="9x15"; void [host]DrawBody(); void [host]DrawText(); void [host]SayGoodBye(); int [host]FirstDrawGroups(int argc, char **argv, int m,                    int *N, Body **Groups)    unsigned width=WIDTH, height=HEIGHT;    int x=0, y=0;    unsigned border_width=4;    unsigned display_width, display_height;    char *window_name="Groups Demo";    char *icon_name="galaxy";    Pixmap icon_pixmap;    XSizeHints size_hints;    unsigned long valuemask=0;    XGCValues values;    char *display_name=NULL;    int i, j, k;    wtime = MPC_Wtime();    if((display = XOpenDisplay(display_name)) == NULL)    {       printf( "FirstDrawGroups: cannot connect to Xserver %s\n",                 XDisplayName(display_name));       return (-1);    }    screen = XDefaultScreen(display);    win = XCreateSimpleWindow(display, XRootWindow                   (display,screen), x, y,                    width, height, border_width,                    XBlackPixel(display, screen),                   XBlackPixel(display, screen)); icon_pixmap = XCreateBitmapFromData(display, win,                    icon_bits,                   icon_width,                    icon_height); size_hints.flags = PPosition|PSize|PMinSize; size_hints.x = x; size_hints.y = y; size_hints.width = width; size_hints.height = height; size_hints.min_width = width; size_hints.min_height = height; XSetStandardProperties(display, win, window_name,                    icon_name, icon_pixmap, argv,                    argc, &size_hints); XSelectInput(display, win, ExposureMask|KeyPressMask|                   ButtonPressMask|StructureNotifyMask); if((font_info = XLoadQueryFont(display, fontname)) == NULL) {    printf("FirstDrawGroups: cannot open %s font\n", fontname);    return -1; } gc=XCreateGC(display, win, valuemask, &values); XSetForeground(display, gc, XWhitePixel(display, screen)); XSetLineAttributes(display, gc, 1, LineSolid,                   CapNotLast, JoinMiter); XSetFillStyle(display, gc, FillSolid); XSetFont(display, gc, font_info->fid); XMapWindow(display, win);          while(1)          {             XNextEvent(display, &report);             switch(report.type)             {             case Expose:                while(XCheckTypedEvent(display, Expose, &report));                for(i=0; i<m; i++)                   for(j=0; j<N[i]; j++)                      if(Groups[i][j].m>0.)                      DrawBody(&gc, (int)(Groups[i][j].m/MASSUNIT),                      (int)((Groups[i][j].p)[X]/POINTSIZE),                      (int)((Groups[i][j].p)[Y]/POINTSIZE),                      (int)((Groups[i][j].p)[Z]/POINTSIZE));                      DrawText(t);                      XFlush(display);                      return 1;                case ConfigureNotify:                      break;                case ButtonPress:                case KeyPress:                   return 1;                default:                   break;       }    } } int [host]DrawGroups (int m, int *N, Body **Groups) {    int i, j, k, size;    XEvent event;    event.type = Expose;    XSendEvent(display, win, False, 0L, &event);    while(1)    {       XNextEvent(display, &report);       switch(report.type)       {          case Expose             while (XCheckTypedEvent (display, Expose, &report))                ;             XClearWindow(display, win);             for(i = 0; i<m; i++)                for(j = 0; j<N[i]; j++)                   if(Groups[i][j].m>0.)                   DrawBody(&gc,(int)(Groups[i][j].m/MASSUNIT),                      (int)((Groups[i][j].p) [X]/POINTSIZE),                      (int)((Groups[i][j].p) [Y]/POINTSIZE),                      (int)((Groups[i][j].p) [Z]/POINTSIZE));             DrawText(t);             XFlush(display);             XSendEvent(display, win, False, 0L, &event);       return 1;          case ConfigureNotify:          break;          case ButtonPress:          case KeyPress:             XFreeGC(display, gc);             XCloseDisplay(display);             return 1;          default:             break;       }    } } void [host]DrawBody(GC *gc, int m, int x, int y, int z) {    int xx, yy, r;    r = (int)(sqrt((double)m)/2);    for(xx=x-r; xx<=x+r; xx++)       for(yy=y-r; yy<=y+r; yy++)          if((xx-x)*(xx-x)+(yy-y)*(yy-y)<=r*r)             XDrawPoint(display, win, *gc, xx, yy); } void [host]DrawText(double t) {    char text[25];    int len, wid;    double gt, wt;    sprintf(text, "Real time = %.1f hours", gt=t/3600.);    len = strlen(text);    wid = XTextWidth(font_info, text, len);    XDrawString(display, win, gc, (int)((WIDTH-wid)/2),                   20, text, len);    sprintf( text, "Wall time = %.1f seconds", _wt=MPC_Wtime()-wtime);    len = strlen(text);    wid = XTextWidth(font_info, text, len);    XDrawString(display, win, gc, (int)((WIDTH-wid)/2),                   40, text, len);    sprintf(text, "Rate = %.1f sec/h", gt?wt/gt:0.);    len = strlen(text);    wid = XTextWidth(font_info, text, len);    XDrawString(display, win, gc, (int)((WIDTH-wid)/2),                   60, text, len); } void [host]SayGoodBye(void) {    char text[50];    int len, wid;    sprintf(text, "Good bye! See you soon.");    len = strlen(text);    wid = XTextWidth(font_info, text, len);    XDrawString(display, win, gc, (int)((WIDTH-wid)/2),                   (int)(HEIGHT/3), text, len);    XFlush(display); }

Parallel Computing on Heterogeneous Networks
Parallel Computing on Heterogeneous Networks (Wiley Series on Parallel and Distributed Computing)
Year: 2005
Pages: 95

flylib.com © 2008-2017.
If you may any questions please contact us: flylib@qtcs.net