source: CIVL/examples/openacc/acc_c3a/acc_c3a.c

/* 
 *     Copyright (c) 2016, NVIDIA CORPORATION.  All rights reserved.
 *
 * NVIDIA CORPORATION and its licensors retain all intellectual property
 * and proprietary rights in and to this software, related documentation
 * and any modifications thereto.  Any use, reproduction, disclosure or
 * distribution of this software and related documentation without an express
 * license agreement from NVIDIA CORPORATION is strictly prohibited.
 *
 */


/*
 * Jacobi iteration example using OpenACC in C
 * Build with
 *   pgcc -acc -Minfo=accel -fast c3.c
 */

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <openacc.h>
#include <math.h>

#if defined(_WIN32) || defined(_WIN64)
#include <sys/timeb.h>
#define gettime(a) _ftime(a)
#define usec(t1,t2) ((((t2).time-(t1).time)*1000+((t2).millitm-(t1).millitm))*100)
typedef struct _timeb timestruct;
#else
#include <sys/time.h>
#define gettime(a) gettimeofday(a,NULL)
#define usec(t1,t2) (((t2).tv_sec-(t1).tv_sec)*1000000+((t2).tv_usec-(t1).tv_usec))
typedef struct timeval timestruct;
#endif

void
smooth( float*restrict a, float*restrict b, float w0, float w1, float w2, int n, int m, int niters )
{
    int i, j, iter;
    float* tmp;
    for( iter = 1; iter <= niters; ++iter ){
        #pragma acc kernels loop present(b[0:n*m],a[0:n*m]) independent
        for( i = 1; i < n-1; ++i ){
            for( j = 1; j < m-1; ++j ){
                a[i*m+j] = w0 * b[i*m+j] + 
                    w1*(b[(i-1)*m+j] + b[(i+1)*m+j] + b[i*m+j-1] + b[i*m+j+1]) +
                    w2*(b[(i-1)*m+j-1] + b[(i-1)*m+j+1] + b[(i+1)*m+j-1] + b[(i+1)*m+j+1]);
                }
        }
        tmp = a;  a = b;  b = tmp;
    }
}

void
smoothhost( float*restrict a, float*restrict b, float w0, float w1, float w2, int n, int m, int niters )
{
    int i, j, iter;
    float* tmp;
    for( iter = 1; iter <= niters; ++iter ){
        for( i = 1; i < n-1; ++i ){
            for( j = 1; j < m-1; ++j ){
                a[i*m+j] = w0 * b[i*m+j] + 
                    w1*(b[(i-1)*m+j] + b[(i+1)*m+j] + b[i*m+j-1] + b[i*m+j+1]) +
                    w2*(b[(i-1)*m+j-1] + b[(i-1)*m+j+1] + b[(i+1)*m+j-1] + b[(i+1)*m+j+1]);
                }
        }
        tmp = a;  a = b;  b = tmp;
    }
}

void
doprt( char* s, float*restrict a, float*restrict ah, int i, int j, int n, int m )
{
    printf( "%s[%d][%d] = %g  =  %g\n", s, i, j, a[i*m+j], ah[i*m+j] );
}

int
main( int argc, char* argv[] )
{
    float *aa, *bb, *aahost, *bbhost;
    int i,j;
    float w0, w1, w2;
    int n, m, aerrs, berrs, iters;
    float dif, rdif, tol;
    timestruct t1, t2, t3;
    long long cgpu, chost;

    n = 0;
    m = 0;
    iters = 0;

    if( argc > 1 ){
        n = atoi( argv[1] );
        if( argc > 2 ){
            m = atoi( argv[2] );
            if( argc > 3 ){
                iters = atoi( argv[3] );
                if( argc > 4 ){
                    if( !strcmp( argv[4], "host" ) ||
                        !strcmp( argv[4], "HOST" ) ){
                        acc_set_device( acc_device_host );
                        printf( "using host\n" );
                    }else
                    if( !strcmp( argv[4], "nvidia" ) ||
                        !strcmp( argv[4], "NVIDIA" ) ){
                        acc_set_device( acc_device_nvidia );
                        acc_init( acc_device_nvidia );
                        printf( "using nvidia\n" );
                    }else{
                        printf( "unknown device: %s\nUsing default\n", argv[4] );
                    }
                }
            }
        }
    }

    if( n <= 0 ) n = 100;
    if( m <= 0 ) m = n;
    if( iters <= 0 ) iters = 10;

    aa = (float*) malloc( sizeof(float) * n * m );
    aahost = (float*) malloc( sizeof(float) * n * m );
    bb = (float*)malloc( sizeof(float) * n * m );
    bbhost = (float*)malloc( sizeof(float) * n * m );
    for( i = 0; i < n; ++i ){
        for( j = 0; j < m; ++j ){
            aa[i*m+j] = 0;
            aahost[i*m+j] = 0;
            bb[i*m+j] = i*1000 + j;
            bbhost[i*m+j] = i*1000 + j;
        }
    }
    w0 = 0.5;
    w1 = 0.3;
    w2 = 0.2;
    gettime( &t1 );
    #pragma acc data copy(bb[0:n*m],aa[0:n*m])
    {
    smooth( aa, bb, w0, w1, w2, n, m, iters );
    }
    gettime( &t2 );
    smoothhost( aahost, bbhost, w0, w1, w2, n, m, iters );
    gettime( &t3 );

    cgpu = usec(t1,t2);
    chost = usec(t2,t3);

    printf( "matrix %d x %d, %d iterations\n", n, m, iters );
    printf( "%13ld microseconds optimized\n", cgpu );
    printf( "%13ld microseconds on host\n", chost );

    aerrs = berrs = 0;
    tol = 0.000005;
    for( i = 0; i < n; ++i ){
        for( j = 0; j < m; ++j ){
            rdif = dif = fabsf(aa[i*m+j] - aahost[i*m+j]);
            if( aahost[i*m+j] ) rdif = fabsf(dif / aahost[i*m+j]);
            if( rdif > tol ){
                ++aerrs;
                if( aerrs < 10 ){
                    printf( "aa[%d][%d] = %12.7e != %12.7e, dif=%12.7e\n", i, j, (double)aa[i*m+j], (double)aahost[i*m+j], (double)dif );
                }
            }
            rdif = dif = fabsf(bb[i*m+j] - bbhost[i*m+j]);
            if( bbhost[i*m+j] ) rdif = fabsf(dif / bbhost[i*m+j]);
            if( rdif > tol ){
                ++berrs;
                if( berrs < 10 ){
                    printf( "bb[%d][%d] = %12.7e != %12.7e, dif=%12.7e\n", i, j, (double)bb[i*m+j], (double)bbhost[i*m+j], (double)dif );
                }
            }
        }
    }
    if( aerrs == 0 && berrs == 0 ){
        printf( "Test PASSED\n" );
        return 0;
    }else{
        printf( "Test FAILED\n" );
        printf( "%d ERRORS found\n", aerrs + berrs );
        return 1;
    }
}