| 1 | #include <stdio.h>
|
|---|
| 2 | #include <assert.h>
|
|---|
| 3 |
|
|---|
| 4 | #ifdef _CIVL
|
|---|
| 5 | #include <civlc.cvh>
|
|---|
| 6 | #define n 10
|
|---|
| 7 | $input double u[n];
|
|---|
| 8 | #else
|
|---|
| 9 | #define n 10
|
|---|
| 10 | #endif
|
|---|
| 11 |
|
|---|
| 12 | int nEdges = n-1;
|
|---|
| 13 |
|
|---|
| 14 | void residualPrllel(double uin[n], double resout[n], int edges[nEdges][2], int colourIA[3]) {
|
|---|
| 15 | for(int c=0; c<2; c++) {
|
|---|
| 16 | #pragma omp parallel for default(none) shared(nEdges,edges,uin,resout,colourIA,c)
|
|---|
| 17 | for(int e=colourIA[c]; e<colourIA[c+1]; e++) {
|
|---|
| 18 | int i = edges[e][0];
|
|---|
| 19 | int j = edges[e][1];
|
|---|
| 20 | resout[i] += uin[i]*uin[j];
|
|---|
| 21 | resout[j] += 2*uin[i]+2*uin[j];
|
|---|
| 22 | }
|
|---|
| 23 | }
|
|---|
| 24 | }
|
|---|
| 25 |
|
|---|
| 26 | void residualSerial(double uin[n], double resout[n], int edges[nEdges][2]) {
|
|---|
| 27 | for(int e=0; e<nEdges; e++) {
|
|---|
| 28 | int i = edges[e][0];
|
|---|
| 29 | int j = edges[e][1];
|
|---|
| 30 | resout[i] = resout[i] + uin[i]*uin[j];
|
|---|
| 31 | resout[j] = resout[j] + 2*uin[i]+2*uin[j];
|
|---|
| 32 | }
|
|---|
| 33 | }
|
|---|
| 34 |
|
|---|
| 35 | int main(int argc, char** argv) {
|
|---|
| 36 | #ifndef _CIVL
|
|---|
| 37 | double u[n];
|
|---|
| 38 | for(int i=0; i<n; i++) {
|
|---|
| 39 | u[i] = 1;
|
|---|
| 40 | }
|
|---|
| 41 | #endif
|
|---|
| 42 | int edgesSerial[nEdges][2];
|
|---|
| 43 | for(int i=0; i<nEdges; i++) {
|
|---|
| 44 | edgesSerial[i][0] = i;
|
|---|
| 45 | edgesSerial[i][1] = i+1;
|
|---|
| 46 | printf("serial edge #%d=(%d, %d)\n", i, i, i+1);
|
|---|
| 47 | }
|
|---|
| 48 | int edgesPrllel[nEdges][2];
|
|---|
| 49 | for(int i=0; i<(nEdges+1)/2; i++) {
|
|---|
| 50 | edgesPrllel[i][0] = 2*i;
|
|---|
| 51 | edgesPrllel[i][1] = 2*i+1;
|
|---|
| 52 | printf("parallel edge #%d=(%d, %d) color A\n", i, 2*i, 2*i+1);
|
|---|
| 53 | }
|
|---|
| 54 | for(int i=1; i<(nEdges+1)/2; i++) {
|
|---|
| 55 | edgesPrllel[(nEdges+1)/2+i-1][0] = 2*i-1;
|
|---|
| 56 | edgesPrllel[(nEdges+1)/2+i-1][1] = 2*i;
|
|---|
| 57 | printf("parallel edge #%d=(%d, %d) color B\n", (nEdges+1)/2+i-1, 2*i-1, 2*i);
|
|---|
| 58 | }
|
|---|
| 59 | int colourIA[3] = {0, (nEdges+1)/2, nEdges};
|
|---|
| 60 | printf("colour markers at %d %d %d\n", colourIA[0], colourIA[1], colourIA[2]);
|
|---|
| 61 | double resPrllel[n];
|
|---|
| 62 | double resSerial[n];
|
|---|
| 63 | for(int i=0; i<n; i++) {
|
|---|
| 64 | resSerial[i] = 0;
|
|---|
| 65 | resPrllel[i] = 0;
|
|---|
| 66 | }
|
|---|
| 67 | residualSerial(u, resSerial, edgesSerial);
|
|---|
| 68 | residualPrllel(u, resPrllel, edgesPrllel, colourIA);
|
|---|
| 69 | for(int i=0; i<n; i++) {
|
|---|
| 70 | printf("residual(%d) = %e or %e\n", i, resPrllel[i], resSerial[i]);
|
|---|
| 71 | assert(resSerial[i] == resPrllel[i]);
|
|---|
| 72 | }
|
|---|
| 73 | }
|
|---|
| 74 |
|
|---|