/*        Generated by TAPENADE     (INRIA, Ecuador team)
    Tapenade 3.12 (r6213) - 13 Oct 2016 10:54
*/
/*
  Differentiation of asg_pstep2d in forward (tangent) mode:
   variations   of useful results: **p0 **p1
   with respect to varying inputs: **v0 **v1 **bulk **p0 **p1
   RW status of diff variables: **v0:in **v1:in **bulk:in **p0:in-out
                **p1:in-out
   Plus diff mem management of: v0:in *v0:in v1:in *v1:in bulk:in
                *bulk:in p0:in *p0:in p1:in *p1:in sdiv:in

void asg_pstep2d(float ** restrict bulk, 
		 float ** restrict p0, float ** restrict p1,
		 float ** restrict v0, float ** restrict v1,
		 float ** restrict ep, float ** restrict epp,
		 float * restrict sdiv,
		 int * gsc, int * gec, 
		 int * lbc, int * rbc,
		 int maxoff, float ** restrict c) {

 signature for tapenade*/

// hand mods - restrict variables, include file
#include "cstd.h"

void asg_pstep2d_d(float ** restrict bulk,
		   float ** restrict bulkd,
		   float ** restrict p0,
		   float ** restrict p0d,
		   float ** restrict p1,
		   float ** restrict p1d,
		   float ** restrict v0,
		   float ** restrict v0d,
		   float ** restrict v1,
		   float ** restrict v1d, 
		   float ** restrict ep,
		   float ** restrict epp,
		   float * restrict sdiv,
		   float * restrict sdivd,
		   int *gsc, int *gec,
		   int *lbc, int *rbc,
		   int maxoff, float ** restrict c) {
  
    int i0, i1;
    int ioff;
    for (i1 = gsc[1]; i1 < gec[1]+1; ++i1) {
        for (i0 = gsc[0]; i0 < gec[0]+1; ++i0) {
            sdivd[i0] = c[0][0]*(v0d[i1][i0]-v0d[i1][i0-1]) + c[1][0]*(v1d[i1]
                [i0]-v1d[i1-1][i0]);
            sdiv[i0] = c[0][0]*(v0[i1][i0]-v0[i1][i0-1]) + c[1][0]*(v1[i1][i0]
                -v1[i1-1][i0]);
        }
        for (ioff = 1; ioff < maxoff; ++ioff)
            for (i0 = gsc[0]; i0 < gec[0]+1; ++i0) {
                sdivd[i0] = sdivd[i0] + c[0][ioff]*(v0d[i1][i0+ioff]-v0d[i1][
                    i0-ioff-1]) + c[1][ioff]*(v1d[i1+ioff][i0]-v1d[i1-ioff-1][
                    i0]);
                sdiv[i0] += c[0][ioff]*(v0[i1][i0+ioff]-v0[i1][i0-ioff-1]) + c
                [1][ioff]*(v1[i1+ioff][i0]-v1[i1-ioff-1][i0]);
            }
        for (i0 = gsc[0]; i0 < gec[0]+1; ++i0) {
            //      fprintf(stderr,"i1=%d i0=%d bulk=%e\n",i1,i0,bulk[i1][i0]);
            p0d[i1][i0] = epp[0][i0]*(ep[0][i0]*p0d[i1][i0]-bulkd[i1][i0]*sdiv
                [i0]-bulk[i1][i0]*sdivd[i0]);
            p0[i1][i0] = (p0[i1][i0]*ep[0][i0]-bulk[i1][i0]*sdiv[i0])*epp[0][
                i0];
            p1d[i1][i0] = epp[1][i1]*(ep[1][i1]*p1d[i1][i0]-bulkd[i1][i0]*sdiv
                [i0]-bulk[i1][i0]*sdivd[i0]);
            p1[i1][i0] = (p1[i1][i0]*ep[1][i1]-bulk[i1][i0]*sdiv[i0])*epp[1][
                i1];
        }
    }
    /* boundary conditions - p is odd about index just before/after comp domain 
    */
    if (lbc[0])
        for (i1 = gsc[1]; i1 < gec[1]+1; ++i1) {
            p0d[i1][gsc[0] - 1] = 0.0;
            p0[i1][gsc[0] - 1] = 0.0f;
            for (ioff = 1; ioff < maxoff; ++ioff) {
                p0d[i1][gsc[0] - ioff - 1] = -p0d[i1][gsc[0]+ioff-1];
                p0[i1][gsc[0] - ioff - 1] = -p0[i1][gsc[0]+ioff-1];
            }
        }
    if (rbc[0])
        for (i1 = gsc[1]; i1 < gec[1]+1; ++i1) {
            p0d[i1][gec[0] + 1] = 0.0;
            p0[i1][gec[0] + 1] = 0.0f;
            for (ioff = 1; ioff < maxoff; ++ioff) {
                p0d[i1][gec[0] + ioff + 1] = -p0d[i1][gec[0]-ioff+1];
                p0[i1][gec[0] + ioff + 1] = -p0[i1][gec[0]-ioff+1];
            }
        }
    if (lbc[1])
        for (i0 = gsc[0]; i0 < gec[0]+1; ++i0) {
            p1d[gsc[1] - 1][i0] = 0.0;
            p1[gsc[1] - 1][i0] = 0.0f;
            for (ioff = 1; ioff < maxoff; ++ioff) {
                p1d[gsc[1] - ioff - 1][i0] = -p1d[gsc[1]+ioff-1][i0];
                p1[gsc[1] - ioff - 1][i0] = -p1[gsc[1]+ioff-1][i0];
            }
        }
    if (rbc[1])
        for (i0 = gsc[0]; i0 < gec[0]+1; ++i0) {
            p1d[gec[1] + 1][i0] = 0.0;
            p1[gec[1] + 1][i0] = 0.0f;
            for (ioff = 1; ioff < maxoff; ++ioff) {
                p1d[gec[1] + ioff + 1][i0] = -p1d[gec[1]-ioff+1][i0];
                p1[gec[1] + ioff + 1][i0] = -p1[gec[1]-ioff+1][i0];
            }
        }
}