Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /* Copyright (c) 2012, NVIDIA CORPORATION. All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions
- * are met:
- * * Redistributions of source code must retain the above copyright
- * notice, this list of conditions and the following disclaimer.
- * * 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.
- * * Neither the name of NVIDIA CORPORATION nor the names of its
- * contributors may be used to endorse or promote products derived
- * from this software without specific prior written permission.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``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 COPYRIGHT OWNER 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 <math.h>
- #include <string.h>
- #include <stdlib.h>
- #include <omp.h>
- #include <stdio.h>
- double simplest_checksum(int imax, int jmax, double (*in)[jmax+2]) {
- double checksum = 0.0;
- int i;
- int j;
- for (i=1; i<(imax+1); i++){
- for (j=1; j<(jmax+1); j++){
- checksum+=in[i][j]*((double)((i)*jmax)+(j));
- }
- }
- return checksum;
- }
- void print_file(int imax, int jmax, double (*in)[jmax+2]) {
- char fname[10];
- sprintf(fname,"out.txt");
- FILE *f=fopen(fname,"w");
- int i;
- int j;
- for (i=0; i<(imax+2); i++){
- for (j=0; j<(jmax+2); j++){
- fprintf(f,"%g ",in[i][j]);
- }
- fprintf(f,"\n");
- }
- fclose(f);
- }
- void func(int imax, int jmax, double (*restrict A)[jmax+2], double (*restrict Anew)[jmax+2]) {
- #pragma acc parallel loop present(A,Anew) collapse(2)
- for( int i = 1; i < imax+1; i++ )
- {
- for( int j = 1; j < jmax+1; j++)
- {
- A[i][j] = Anew[i][j];
- }
- }
- }
- int main(int argc, char** argv)
- {
- //Size along y
- int jmax = 4094;
- //Size along x
- int imax = 4094;
- int iter_max = 1000;
- const double pi = 2.0 * asin(1.0);
- const double tol = 1.0e-5;
- double error = 1.0;
- double (*restrict A)[jmax+2];
- double (*restrict Anew)[jmax+2];
- double *restrict y0;
- A = (double (*)[jmax+2])malloc((imax+2) * (jmax+2) * sizeof(double));
- Anew = (double (*)[jmax+2])malloc((imax+2) * (jmax+2) * sizeof(double));
- y0 = (double *)malloc((imax+2) * sizeof(double));
- memset(A, 0, (imax+2) * (jmax+2) * sizeof(double));
- // set boundary conditions
- for (int j = 0; j < jmax+2; j++)
- A[0][j] = 0.0;
- for (int j = 0; j < jmax+2; j++)
- A[imax+1][j] = 0.0;
- for (int i = 0; i < imax+2; i++)
- {
- y0[i] = sin(pi * (i) / (imax+1));
- A[i][0] = y0[i];
- }
- for (int i = 0; i < imax+2; i++)
- {
- y0[i] = sin(pi * (i) / (imax+1));
- A[i][jmax+1] = y0[i]*exp(-pi);
- }
- printf("Jacobi relaxation Calculation: %d x %d mesh\n", imax+2, jmax+2);
- double t1 = omp_get_wtime();
- int iter = 0;
- for (int j = 1; j < jmax+2; j++)
- Anew[0][j] = 0.0;
- for (int j = 1; j < jmax+2; j++)
- Anew[imax+1][j] = 0.0;
- for (int i = 1; i < imax+2; i++)
- Anew[i][0] = y0[i];
- for (int i = 1; i < imax+2; i++)
- Anew[i][jmax+1] = y0[i]*expf(-pi);
- #pragma acc data copy(A[0:imax+2] [0:jmax+2], Anew[0:imax+2] [0:jmax+2])
- while ( error > tol && iter < iter_max )
- {
- error = 0.0;
- #pragma acc parallel loop reduction(max:error)
- for( int i = 1; i < imax+1; i++ )
- {
- for( int j = 1; j < jmax+1; j++)
- {
- Anew[i][j] = 0.25f * ( A[i][j+1] + A[i][j-1]
- + A[i-1][j] + A[i+1][j]);
- error = fmax( error, fabs(Anew[i][j]-A[i][j]));
- }
- }
- func(imax, jmax, A, Anew);
- if(iter % 100 == 0) printf("%5d, %0.8f\n", iter, error);
- iter++;
- }
- double checksum = simplest_checksum(imax,jmax,A);
- printf("Checksum: %4.8f\n",checksum);
- double runtime = omp_get_wtime()-t1;
- printf(" total: %f s\n", runtime);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement