convo.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>

#define A(i,j) a[((i)*N+(j))]

#define OUT(i,j) out[((i)*N+(j))]


void convolve3(int M, int N, double *a, double scale1, double scale2, double *invr, double *out)
{
        int i, j;

        // Calculate the inner part of the array first, 
        // then the edges and corners later 
        for (i=1; i<M-1; i++) 
                for (j=2; j<N-1; j++) 
                        OUT(i,j) = scale1 * (
                                                  A(i-1,j) 
                                + A(i,j-1) - 4. * A(i,  j) + A(i,j+1)
                                                + A(i+1,j)   )
                        + scale2 * ( (-A(i,j-1) + A(i,j+1))*invr[j]);

        /* j=1: This is the r=0 row, so use L0 rather than L */
        for (i=1; i<M-1; i++) 
                OUT(i,1) = scale1 * (
                                                     A(i-1,1) 
                                + 2. * A(i,0) - 6. * A(i,  1) + 2. * A(i,2)
                                                   + A(i+1,1)   );

        // j=1, i=0 
        OUT(0,1) = scale1 * (
                                2. * A(0,0) - 6. * A(0,1) + 2. * A(0,2)
                                                 + A(1,1)   );

        // j=1, i=M-1 
        OUT(M-1,1) = scale1 * (
                                                     A(M-2,1) 
                                + 2. * A(M-1,0) - 6. * A(M-1,  1) + 2. * A(M-1,2)
                                                   );
        /* End of j=1 section */


        // i=0 
        for (j=1; j<N-1; j++) 
                OUT(0,j) = scale1 * (
                          A(0,j-1) - 4. * A(0,j) + A(0,j+1)
                                        + A(1,j)   )
                        + scale2 * ( (-A(0,j-1) + A(0,j+1))*invr[j] );

        // i=M-1 
        for (j=1; j<N-1; j++) 
                OUT(M-1,j) = scale1 * (
                                             A(M-2,j) 
                        + A(M-1,j-1) - 4. * A(M-1,j) + A(M-1,j+1)   )
                        + scale2 * ( (-A(M-1,j-1) + A(M-1,j+1))*invr[j] );

        // j=0 
        for (i=1; i<M-1; i++) 
                OUT(i,0) = scale1 * (
                               A(i-1,0) 
                        - 4. * A(i,  0) + A(i,1)
                             + A(i+1,0)   )
                        + scale2 * ( A(i,1)*invr[0] );

        // j=N-1 
        for (i=1; i<M-1; i++) 
                OUT(i,N-1) = scale1 * (
                                           A(i-1,N-1) 
                        + A(i,N-2) - 4. * A(i,  N-1) 
                                         + A(i+1,N-1)   )
                        + scale2 * ( -A(i,N-2)*invr[N-1] );

        // i=0, j=0 
        OUT(0,0) = scale1 * (
                        - 4. * A(0,0) + A(0,1)
                             + A(1,0)   )
                        + scale2 * ( A(0,1)*invr[0] );

        // i=0, j=N-1 
        OUT(0,N-1) = scale1 * (
                          A(0,N-2) - 4. * A(0,N-1) 
                                         + A(1,N-1)   )
                        + scale2 * ( -A(0,N-2)*invr[N-1] );

        // i=M-1, j=0 
        OUT(M-1,0) = scale1 * (
                               A(M-2,0) 
                        - 4. * A(M-1,0) + A(M-1,1)   )
                        + scale2 * ( A(M-1,1)*invr[0] );

        // i=M-1, j=N-1 
        OUT(M-1,N-1) = scale1 * (
                                              A(M-2,N-1) 
                        + A(M-1,N-2) - 4. * A(M-1,N-1)   )
                        + scale2 * ( -A(M-1,N-2)*invr[N-1] );
}