rendered paste body/*- * Copyright (c) 2006 pluknet. All rights reserved. * Copyright (c) 2006 fenix. All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */#include <sys/time.h>#include <err.h>#include <mpi.h>#include <stdio.h>#include <stdlib.h>#include <string.h>#include <sysexits.h>#include <unistd.h>static void usage(int);intmain(int argc, char *argv[]){ int numprocs, myid, ch; double *a, *b, *c; int n = 1000000; // horyzon nsteps int icount = n; double time = 0.00001; int i, iter = 0; // counters#ifndef QUIETLY int j; // counters#endif double k = 0.0001; // heat conduction double h; // horyzon step int remain; // distribute remain points to procs double tau; int cflag = 0; int ncur; // # of points in cur proc#ifndef QUIETLY int ncurn = 0;#endif int n2send = 0; // # of points to send/recv for neighbor proc struct timeval before, after; // TIME struct { struct timeval x; struct timeval y; } balance[3]; // TIMING int bnum[3]; // distribution of points among neighbor procs#define SELF 1#define RIGHT 2#define LEFT 0 MPI_Status stat; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &numprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myid); while ((ch = getopt(argc, argv, "k:n:t:c:")) != -1) switch (ch) { case 'c': cflag = 1; icount = atoi(optarg); break; case 'k': k = atof(optarg); if (k <= 0) { if (myid == 0) warnx("invalid argument to option k: %s", optarg); usage(myid); } break; case 'n': n = atoi(optarg); if (n <= numprocs) { if (myid == 0) warnx("invalid argument to option n: %s", optarg); usage(myid); } break; case 't': time = atof(optarg); if (time <= 0) { if (myid == 0) warnx("invalid argument to option t: %s", optarg); usage(myid); } break; default: usage(myid); } argc -= optind; argv += optind; if (argc > 4) usage(myid); h = 1.0 / (n - 1); tau = 0.5 * h * h / k; if (!cflag) icount = time/tau; else if (icount > time/tau) { if (myid == 0) warnx("invalid argument to option c: %s, set at least: %d", optarg, (int) (time/tau)); usage(myid); } ncur = n / numprocs; remain = n - numprocs * ncur; ncur += 2; while (remain--) if (remain == myid) ncur++; a = (double *)calloc(ncur, sizeof(double)); b = (double *)calloc(ncur, sizeof(double)); bzero(a, ncur * sizeof(double)); bzero(b, ncur * sizeof(double)); bzero(bnum, 3 * sizeof(int)); for (i = 0; i < 3; i++) { balance[i].x.tv_sec = 0; balance[i].x.tv_usec = 0; balance[i].y.tv_sec = 0; balance[i].y.tv_usec = 0; } if (myid == numprocs - 1) { a[ncur - 2] = 1.0; b[ncur - 2] = 1.0; } if (myid == 0) {#ifdef TIME gettimeofday(&before, NULL);#endif#ifdef DEBUG printf("# of iterations: %f\n", (time / tau)); printf("Time = %f, tau = %.10f\n", time, tau);#endif } while ((time -= tau) > 0) {#ifdef TIMING iter++;#endif if (!(myid % 2)) { if (myid != numprocs - 1) { if (iter == icount) {#ifdef BALANCE MPI_Send(&balance[SELF].x, sizeof (struct timeval), MPI_CHAR, myid + 1, 0, MPI_COMM_WORLD); MPI_Send(&bnum[SELF], 1, MPI_INT, myid + 1, 0, MPI_COMM_WORLD);#ifdef DEBUG printf("%d->%d %ld:%ld\n", myid, myid + 1, balance[SELF].x.tv_sec, balance[SELF].x.tv_usec); printf("%d->%d %d\n", myid, myid + 1, bnum[SELF]);#endif#endif } MPI_Send(&a[ncur - 2], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD); if (iter == icount) {#ifdef BALANCE MPI_Recv(&balance[RIGHT].x, sizeof (struct timeval), MPI_CHAR, myid + 1, 0, MPI_COMM_WORLD, &stat); MPI_Recv(&bnum[RIGHT], 1, MPI_INT, myid + 1, 0, MPI_COMM_WORLD, &stat);#ifdef DEBUG printf("%d<-%d %ld:%ld\n", myid, myid + 1, balance[RIGHT].x.tv_sec, balance[RIGHT].x.tv_usec); printf("%d<-%d %d\n", myid, myid + 1, bnum[RIGHT]);#endif#endif } MPI_Recv(&a[ncur - 1], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD, &stat); } if (myid != 0) { if (iter == icount) {#ifdef BALANCE MPI_Send(&balance[SELF].x, sizeof (struct timeval), MPI_CHAR, myid - 1, 0, MPI_COMM_WORLD); MPI_Send(&bnum[SELF], 1, MPI_INT, myid - 1, 0, MPI_COMM_WORLD);#ifdef DEBUG printf("%d<-%d %ld:%ld\n", myid - 1, myid, balance[SELF].x.tv_sec, balance[SELF].x.tv_usec); printf("%d<-%d %d\n", myid - 1, myid, bnum[SELF]);#endif#endif } MPI_Send(&a[1], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD); if (iter == icount) {#ifdef BALANCE MPI_Recv(&balance[LEFT].x, sizeof (struct timeval), MPI_CHAR, myid - 1, 0, MPI_COMM_WORLD, &stat); MPI_Recv(&bnum[LEFT], 1, MPI_INT, myid - 1, 0, MPI_COMM_WORLD, &stat);#ifdef DEBUG printf("%d->%d %ld:%ld\n", myid - 1, myid, balance[LEFT].x.tv_sec, balance[LEFT].x.tv_usec); printf("%d->%d %d\n", myid - 1, myid, bnum[LEFT]);#endif#endif } MPI_Recv(&a[0], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD, &stat); } } else { if (myid != 0) { if (iter == icount) {#ifdef BALANCE MPI_Recv(&balance[LEFT].x, sizeof (struct timeval), MPI_CHAR, myid - 1, 0, MPI_COMM_WORLD, &stat); MPI_Recv(&bnum[LEFT], 1, MPI_INT, myid - 1, 0, MPI_COMM_WORLD, &stat);#ifdef DEBUG printf("%d->%d %ld:%ld\n", myid - 1, myid, balance[LEFT].x.tv_sec, balance[LEFT].x.tv_usec); printf("%d->%d %d\n", myid - 1, myid, bnum[LEFT]);#endif#endif } MPI_Recv(&a[0], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD, &stat); if (iter == icount) {#ifdef BALANCE MPI_Send(&balance[SELF].x, sizeof (struct timeval), MPI_CHAR, myid - 1, 0, MPI_COMM_WORLD); MPI_Send(&bnum[SELF], 1, MPI_INT, myid - 1, 0, MPI_COMM_WORLD);#ifdef DEBUG printf("%d<-%d %ld:%ld\n", myid - 1, myid, balance[SELF].x.tv_sec, balance[SELF].x.tv_usec); printf("%d<-%d %d\n", myid - 1, myid, bnum[SELF]);#endif#endif } MPI_Send(&a[1], 1, MPI_DOUBLE, myid - 1, 0, MPI_COMM_WORLD); } if (myid != numprocs - 1) { if (iter == icount) {#ifdef BALANCE MPI_Recv(&balance[RIGHT].x, sizeof (struct timeval), MPI_CHAR, myid + 1, 0, MPI_COMM_WORLD, &stat); MPI_Recv(&bnum[RIGHT], 1, MPI_INT, myid + 1, 0, MPI_COMM_WORLD, &stat);#ifdef DEBUG printf("%d<-%d %ld:%ld\n", myid, myid + 1, balance[RIGHT].x.tv_sec, balance[RIGHT].x.tv_usec); printf("%d<-%d %d\n", myid, myid + 1, bnum[RIGHT]);#endif#endif } MPI_Recv(&a[ncur - 1], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD, &stat); if (iter == icount) {#ifdef BALANCE MPI_Send(&balance[SELF].x, sizeof (struct timeval), MPI_CHAR, myid + 1, 0, MPI_COMM_WORLD); MPI_Send(&bnum[SELF], 1, MPI_INT, myid + 1, 0, MPI_COMM_WORLD);#ifdef DEBUG printf("%d->%d %ld:%ld\n", myid, myid + 1, balance[SELF].x.tv_sec, balance[SELF].x.tv_usec); printf("%d->%d %d\n", myid, myid + 1, bnum[SELF]);#endif#endif } MPI_Send(&a[ncur - 2], 1, MPI_DOUBLE, myid + 1, 0, MPI_COMM_WORLD); } }#ifdef BALANCE// n2send = 1/4 * (tC - tL) / (tC + tL) * min(#C, #L); #endif#ifdef BALANCE#ifdef TIMING if (iter == icount) { iter = 0; printf ("%f %d %ld %ld %ld\t%ld %ld %ld\n", time, myid, balance[SELF].x.tv_usec, balance[SELF].y.tv_usec, balance[SELF].y.tv_usec - balance[SELF].x.tv_usec, balance[LEFT].x.tv_usec, balance[SELF].x.tv_usec, balance[RIGHT].x.tv_usec); } gettimeofday(&balance[SELF].x, NULL);#endif#endif if (((ncur == 3) && (myid != 0) && (myid != numprocs - 1)) || ((ncur > 3) && (myid != 0))) b[1] = a[1] + k * tau / h / h * (a[2] + a[0] - 2 * a[1]); for (i = 2; i < ncur - 2; ++i) b[i] = a[i] + k * tau / h / h * (a[i + 1] + a[i - 1] - 2 * a[i]); if ((ncur > 3) && (myid != numprocs - 1)) b[ncur - 2] = a[ncur - 2] + k * tau / h / h * (a[ncur - 1] + a[ncur - 3] - 2 * a[ncur - 2]); c = a; a = b; b = c;#ifdef TIMING gettimeofday(&balance[SELF].y, NULL);#endif }#ifdef TIME if (myid == 0) gettimeofday (&after, NULL);#endif #ifndef QUIETLY if (myid == 0) { ncurn += (ncur - 2); MPI_Send (&ncurn, 1, MPI_INT, myid + 1, 0, MPI_COMM_WORLD); } for (i = 1; i < numprocs - 1; ++i) if (myid == i) { MPI_Recv (&ncurn, 1, MPI_INT, myid - 1, 0, MPI_COMM_WORLD, &stat); ncurn += (ncur - 2); MPI_Send (&ncurn, 1, MPI_INT, myid + 1, 0, MPI_COMM_WORLD); } if (myid == numprocs - 1) { MPI_Recv (&ncurn, 1, MPI_INT, myid - 1, 0, MPI_COMM_WORLD, &stat); ncurn += (ncur - 2); } if (myid == 0) { for (i = 1; i < ncur - 1; ++i) printf ("b[%d]=%f\n", i + ncurn - ncur + 2 - 1, b[i]); for (j = 1; j < numprocs; ++j) { MPI_Recv (&ncur, 1, MPI_INT, j, 0, MPI_COMM_WORLD, &stat); MPI_Recv ( b, ncur, MPI_DOUBLE, j, 0, MPI_COMM_WORLD, &stat); ncurn += (ncur - 2); for (i = 1; i < ncur - 1; ++i) printf ("b[%d]=%f\n", i + ncurn - ncur + 2 - 1, b[i]); } } else { MPI_Send (&ncur, 1, MPI_INT, 0, 0, MPI_COMM_WORLD); MPI_Send ( b, ncur, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD); }#endif#ifdef TIME if (myid == 0) printf ("%ld.%ld\n", after.tv_usec > before.tv_usec ? after.tv_sec - before.tv_sec : after.tv_sec - before.tv_sec - 1, after.tv_usec > before.tv_usec ? after.tv_usec - before.tv_usec : 1000000 - (before.tv_usec - after.tv_usec));#endif MPI_Finalize(); return EXIT_SUCCESS;}static voidusage (int myid){ if (myid == 0) (void)fprintf(stderr, "usage: stereo [-k rate] [-n number] [-t time] [-c icount]\n"); MPI_Finalize(); exit (EX_USAGE);}/* * #balance = 1/4 * (tC - tL) / (tC + tL) * min(#C, #L); */