C.2. USER S GUIDE

C.1. SOURCE CODE

The text of the mpC parallel adaptive quadrature routine presented in Section 9.3.4 is broken down into five source files.

The first file contains the most principal code of the application. The name of this file must have extension mpc (e.g., pquanc8.mpc). The contents of this file are as follows:

#include <stdlib.h> #include <stdlib.h> #include <mpc.h> #include <math.h> #include "quanc8.h" #include "partition.h" typedef struct Interval {    double *x, *f;    double qprev;    struct Interval *next; } *IntervalList; #include "interval.h" double [*]pquanc8(double (*fun)(double), double [host]a,           double [host]b, double [host]abserr,           double [host]relerr, double *[host]errest, int *[host]nofun,           int *[host]flag, double *[host] badpoint) {    double [host]area=0., [*]result=0., [host]absarea,              [host]stone, [host]step, [host]minstep,              [host]cor11=0., [host]esterr, [host]tollerr,              [host]qleft, [host]qright, [host]qnow,              [host]qdiff, [host]prevleft=0.;    int [host]repeat=0, [host]maxrepeat=2;    double *[host]f, *[host]x;    IntervalList [host]ints, [host]nextints;    repl int p;    repl double *speeds;    if(a==b)       return result;    // Initialization for first interval    [host]:    {       *flag = 0;       *errest = 0.;       *nofun = 0;       stone = (b-a)/16.;       minstep = stone*relerr+abserr;       InitFirstInterval(fun, a, b, &x, &f);       *nofun=9;    }    MPC_Processors_static_info(&p, &speeds);    {       net SimpleNet(p) s;       int [s]myn;       repl [s]nints, *[s]ns, *[s]dsp, [s]src=0;       double *[s]nodes, *[host]times;       double *[s]mynodes, *[s]myvalues, [s]mytime;       [s]:       {          ns = malloc(p*sizeof(int));          dsp = malloc(p*sizeof(int));       }       [host]: times = malloc(p*sizeof(double));       ints = ([host]CreateInterval)(x, f, 0., NULL);       nints = 1;       //Main loop of the algorithm       for(step=stone; ints!=NULL; ints=nextints,              step/=2.0)       {          // Calculate points, where integrand should be evaluated          // at this iteration, proceeding from list ints of nints          // subintervals of width step, contribution of which into          // the integral is still to be estimated. Store the points          // in array nodes.          [host]:          {             int j, k;             struct Interval *i;             nodes = malloc(nints*8*sizeof(double));             for(i=ints, k=0; i!=NULL; i=i->next)                for (j=1; j<16; j+=2, k++)                   nodes[k] = i->x[j] = (i->x[j-1]+i->x[j+1])/2.;          }          nints = [host]nints;          // Become alert if the left end point of left subinterval          // does not change from iteration to iteration          [host]:             if(fabs(prevleft-nodes[0])<relerr*fabs(prevleft)                  || fabs(prevleft-nodes[0])<abserr)             repeat++;             else             {                prevleft = nodes[0];                repeat = 0;             }          // Convergence test          if(repeat>maxrepeat && step<minstep)          {             *flag = 1;             *badpoint = [host]nodes[0];             [s]: break;          }          [s]:          {             int k;             // Calculate the number of points assigned to each             // processor and store it in array ns; the number is             // proportional to the processor’s relative speed             Partition(p, speeds, ns, nints*8);             // Distribute the points across processors             myn = ns[I coordof myn];             mynodes = malloc(myn*sizeof(double));             myvalues = malloc(myn*sizeof(double));             for(k=1, dsp[0]=0; k < p; k++)                dsp[k] = ns[k-1]+dsp[k-1];             ([(p)s])MPC_Scatter(&src, nodes, dsp, ns, myn, mynodes);             // Evaluate the integrand at the points in parallel and             // measure the time taken by each processor             if(myn)             {                mytime = MPC_Wtime();                for(k=0; k<myn; k++)                   myvalues[k] = (*fun)(mynodes[k]);                mytime = MPC_Wtime()-mytime;             }             else                mytime = 0.;          }          // Calculate the actual relative speed of processors          // and store it in array speeds          times[] = mytime;          [host]:          {             int j;             double temp;             for(j=0, temp=0; j<p; j++)                if(ns[j] && times[j] > temp*ns[j])                   temp = times[j]/ns[j];             for(j=0; j<p; j++)                if(ns[j])                   speeds[j] = temp/times[j]*1000.*ns[j];          }          [s]: ([(p)s])MPC_Bcast(&src, speeds, 1, p, speeds, 1);          // Gather the function values to host-processor storing          // them in array nodes          ([([s]p)s])MPC_Gather( &src, nodes, dsp, ns, myn, myvalues);          // Complete computation of all elements of list ints that          // represents a set of subintervals of width step, whose          // contribution into the integral is still to be estimated          [host]:          {             int j, k;             struct Interval *i;             for(i=ints, k=0; i!=NULL; i=i->next)                for (j=1; j<16; j+=2, k++)                   i->f[j] = nodes[k];          }          // Increase number of function values          *nofun += ([host]nints)*8;          [s]free(mynodes);          [s]free(myvalues);          [host]: free(nodes);          // Calculate an estimate of the integral over each          // interval in list ints. If the interval is acceptable,          // add its contribution to the result. If not, bisect it          // and add 2 new subintervals to a list of subintervals to          // be processed at next iteration of the main loop          [host]:          {             struct Interval *i;             for (i=ints, nextints=NULL, nints=0; i!=NULL; i=i->next)             {                CentralQuanc8Calculation(i, step, &qleft, &qright);                qnow = qleft+qright;                qdiff = qnow - i->qprev;                area = area+qdiff;                esterr = fabs(qdiff)/1023.0;                absarea = fabs(area);                tollerr = fmax(abserr, relerr*absarea)*(step/stone);                if(esterr <= tollerr) //Interval convergence test                { // Interval converged. Add contributions into sums.                   result += qnow;                   *errest += esterr;                   cor11 += qdiff/1023.0;                   free(i->x);                   free(i->f);                }                else                { // No convergence. Bisect interval i and add its                   // subintervals to list nextints of intervals to be                   // processed at next iteration of the main loop                   nextints = BisectInterval(i, nextints, qleft, qright);                   nints += 2;                }             }          }          [host]FreeList(ints);       }    }    [host]:    {       double temp;       result += cor11;       temp = fabs(result)+*errest;       while (*errest!=0. && temp == fabs(result))       {          *errest *= 2.;          temp = fabs(result)+*errest;       }    }    return result; } #include "function.h" int [*] main() {    double [host]a=0., [host]b=2., [host]abserr=0.,             [host]relerr=1.e-11, [host]errest, [host]badpoint,             [host]sflag, [host]time, result;    int [host]nofun, [host]pflag;    time = ([host]MPC_Wtime)();    result = pquanc8( fun, a, b, abserr, relerr, &errest,  &nofun, &pflag, &badpoint);    time = ([host]MPC_Wtime)()-time;    [host]:       if(pflag)          printf("pquanc8:\nBad singularity at %e\n", badpoint);       else          printf("pquanc8:\nresult=%e errest=%e nofun=%d time=%es\n",                   result, errest, nofun, time);    time = ([host]MPC_Wtime)();    [host]: result = quanc8(fun, a, b, abserr, relerr,                                &errest, &nofun, &sflag);    time = ([host]MPC_Wtime)()-time;    [host]:       if(sflag)          printf("quanc8:\nBad singularity at %e\n",                   b-modf(sflag,NULL)*(b-a));       else          printf("quanc8:\nresult=%e errest=%e nofun=%d time=%es\n",                   result, errest, nofun, time);    }

The name of the second source file must be partition.h. This file contains serial routine code for calculation of the number of points to be assigned to each processor and is included in the program with an include directive. The contents of this file are as follows:

// Calculates distribution of n units of computation across // p processors proportional to their powers    void Partition (int p, double *powers, int *ns, int n) {    int k, sum, imaxfrac, max, imax;    double totalpower, *fracs, maxfrac;    fracs = calloc(p, sizeof(double));    for(k=0, totalpower=0.0; k<p; k++)       totalpower += powers[k];    for(k=0, sum=0; k<p; sum+=ns[k++])    {       ns[k] = (int)((powers[k]/totalpower)*n);       fracs[k] = (powers[k]/totalpower)*n-ns[k];    }    sum = n-sum;    while(sum>0)    {       for(k=0, maxfrac=0., imaxfrac=0; k<p; k++)          if(maxfrac<fracs[k])          {             maxfrac = fracs[k];             imaxfrac = k;          }       ns[imaxfrac]++;       sum––;       fracs[imaxfrac]=-1.0;    }    free(fracs);    for(k=0, sum=0; k<p; k++)       if(!ns[k])       {          ns[k]++;          sum++;       }    while(sum>0) {       for(k=0, max=0, imax=0; k<p; k++)          if(max<ns[k])          {             max = ns[k];             imax = k;          }       ns[imax]––;       sum––;    } }

The name of the third source file must be interval.h. This file contains serial routine code for processing lists of subintervals and is included in the program with an include directive. The contents of this file are as follows:

 struct Interval *CreateInterval(double *x, double *f,                                    double qprev,                                    struct Interval *next) {    struct Interval *i;    i = malloc(sizeof(struct Interval));    i->x = x;    i->f = f;    i->qprev = qprev;    i->next = next;    return i; } IntervalList BisectInterval(struct Interval *i,                                IntervalList ints,                                double qleft, double qright) {    struct Interval *left, *right, *cur;    int j;    right = CreateInterval(malloc(17*sizeof(double)),                                  malloc(17*sizeof (double)),                                  qright, NULL);    for(j=0; j<9; j++)    {       right->x[2*j] = i->x[j+8];       right->f[2*j] = i->f[j+8];    }    left = CreateInterval(i->x, i->f, qleft, right);    for(j=0; j>-9; j––)    {       left->x[2*j+16] = i->x[j+8];       left->f[2*j+16] = i->f[j+8];    }    if(ints == NULL)       return left;    for( cur = ints; cur->next != NULL; cur = cur->next); cur->next = left;    return ints; } void FreeList(IntervalList list) {    IntervalList l;    for(; list!=NULL; list=l)    {       l = list->next;       free(list);    } } void CentralQuanc8Calculation(struct Interval *i,                                double step,                               double *qleft,                                double *qright) {    static double w[] = {3956./14175., 23552./14175.,                               -3712./14175., 41984./14175.,                               -18160./14175.};    *qleft = (w[0]*(i->f[0]+i->f[8])+                    w[1]*(i->f[1]+i->f[7])+                    w[2]*(i->f[2]+i->f[6])+                    w[3]*(i->f[3]+i->f[5])+                    w[4]*i->f[4]                 )*step;    *qright = (w[0]*(i->f[8]+i->f[16])+                      w[1]*(i->f[9]+i->f[15])+                      w[2]*(i->f[10]+i->f[14])+                      w[3]*(i->f[11]+i->f[13])+                      w[4]*i->f[12]                   )*step; } void InitFirstInterval(double (*fun)(double), double a,                                    double b, double **x, double **f) {       int j;       *x = malloc(17*sizeof(double));       (*x)[0] = a;       (*x)[16] = b;       (*x)[8] = ((*x)[0]+(*x)[16])/2.;       (*x)[4] = ((*x)[0]+(*x)[8])/2.;       (*x)[12] = ((*x)[8]+(*x)[16])/2.;       (*x)[2] = ((*x)[0]+(*x)[4])/2.;       (*x)[6] = ((*x)[4]+(*x)[8])/2.;       (*x)[10] = ((*x)[8]+(*x)[12])/2.;       (*x)[14] = ((*x)[12]+(*x)[16])/2.;       *f = malloc(17*sizeof(double));       for(j=0; j<9; j++)          (*f)[2*j] = (*fun)((*x)[2*j]); }

The name of the fourth source file must be function.h. This file sup-plies a function fun for the integrand f(x) and is included in the program with an include directive. An example of the contents of this file that provides a function fun for the integrand is as follows:

 double infun(double x) {    return fabs(sin(16*x)); } double fun(double x) {    double y, errest, flag;    int nofun;    y = quanc8( infun, 0., x+0.5, 0., 1.e-10, &errest, &nofun, &flag);    return x/y; } 

The name of the fifth source file must be quanc8.h. This file contains the quanc8 C code presented in Section 9.3.3.