Index
Subject
: LUG: Coding Horrors -- post your evils
From
: Alexander Ray <alexjray.ncsu@gmail.[redacted]>
Date
: Thu, 06 Aug 2009 02:23:39 -0400
Okay, so I've been wrestling w/ this for a good while, so I thought I'd post it.
It seems especially funny at 2AM, so here it is
I'm up for any suggestions on how to clean this up, but at this point
there are actually more than a dozen routines that look like this.
void Aniso_Pa_Field_Update(int z_min, int z_max, int x_min, int x_max,
COMPLEX** pxa, COMPLEX** pya, COMPLEX** pza,
COMPLEX** past_pxa, COMPLEX** past_pya, COMPLEX** past_pza,
COMPLEX** past_qx, COMPLEX** past_qy, COMPLEX** past_qz,
COMPLEX** rxx, COMPLEX** rxy, COMPLEX** rxz,
COMPLEX** ryx, COMPLEX** ryy, COMPLEX** ryz,
COMPLEX** rzx, COMPLEX** rzy, COMPLEX** rzz,
COMPLEX** cepx1, COMPLEX** cepx2, COMPLEX** cepx3,
COMPLEX** cepy1, COMPLEX** cepy2, COMPLEX** cepy3,
COMPLEX** cepz1, COMPLEX** cepz2, COMPLEX** cepz3,
double** deltaxx, double** deltaxy, double** deltaxz,
double** deltayx, double** deltayy, double** deltayz,
double** deltazx, double** deltazy, double** deltazz,
COMPLEX** cxi11, COMPLEX** cxi12, COMPLEX** cxi13,
COMPLEX** cxi21, COMPLEX** cxi22, COMPLEX** cxi23,
COMPLEX** cxi31, COMPLEX** cxi32, COMPLEX** cxi33,
double th, double S)
{
int z=0, x=0;
double beta;
beta=0.5/pow(cos(th),2);
#pragma omp parallel for schedule(static) private(x,z) default(shared)
for (z=z_min;z<z_max;z++)
{
for(x=x_min;x<x_max;x++)
{
//Pxa Update
if (x>0){
rxx[z][x].r=pxa[z][x].r
-cepx1[z-z_min][x-x_min].r*S*(past_qy[z+1][x].r-past_qy[z][x].r)
+0.5*cepx2[z-z_min][x-x_min].r*S*(past_qx[z+1][x].r+past_qx[z+1][x-1].r-past_qx[z][x].r-past_qx[z][x-1].r)
-0.5*cepx2[z-z_min][x-x_min].r*S*(past_qz[z][x+1].r-past_qz[z][x-1].r)
+0.25*cepx3[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].r+past_qy[z][x+1].r-past_qy[z+1][x-1].r-past_qy[z][x-1].r)
+0.5*sin(th)*(deltaxx[z-z_min][x]*cepx3[z-z_min][x].r+deltaxy[z-z_min][x]*cepy3[z-z_min][x].r+deltaxz[z-z_min][x]*cepz3[z-z_min][x].r)
*(past_qy[z][x].r+past_qy[z+1][x].r)
-sin(th)*(deltaxx[z-z_min][x]*cepx2[z-z_min][x].r+deltaxy[z-z_min][x]*cepy2[z-z_min][x].r+deltaxz[z-z_min][x]*cepz2[z-z_min][x].r)*past_qz[z][x].r
-deltaxx[z-z_min][x]*((1-2*beta)*past_pxa[z][x].r
+beta*pxa[z][x].r)
-0.5*deltaxy[z-z_min][x]*((1-2*beta)*(past_pya[z][x].r+past_pya[z][x-1].r)
+beta*(pya[z][x].r+pya[z][x-1].r))
-0.25*deltaxz[z-z_min][x]*((1-2*beta)*(past_pza[z][x].r+past_pza[z][x-1].r+past_pza[z+1][x].r+past_pza[z+1][x-1].r)
+beta*(pza[z][x].r+pza[z][x-1].r+pza[z+1][x].r+pza[z+1][x-1].r));
rxx[z][x].i=pxa[z][x].i
-cepx1[z-z_min][x-x_min].r*S*(past_qy[z+1][x].i-past_qy[z][x].i)
+0.5*cepx2[z-z_min][x-x_min].r*S*(past_qx[z+1][x].i+past_qx[z+1][x-1].i-past_qx[z][x].i-past_qx[z][x-1].i)
-0.5*cepx2[z-z_min][x-x_min].r*S*(past_qz[z][x+1].i-past_qz[z][x-1].i)
+0.25*cepx3[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].i+past_qy[z][x+1].i-past_qy[z+1][x-1].i-past_qy[z][x-1].i)
+0.5*sin(th)
*(deltaxx[z-z_min][x]*cepx3[z-z_min][x].r+deltaxy[z-z_min][x]*cepy3[z-z_min][x].r+deltaxz[z-z_min][x]*cepz3[z-z_min][x].r)
*(past_qy[z][x].i+past_qy[z+1][x].i)
-sin(th)
*(deltaxx[z-z_min][x]*cepx2[z-z_min][x].r+deltaxy[z-z_min][x]*cepy2[z-z_min][x].r+deltaxz[z-z_min][x]*cepz2[z-z_min][x].r)
*past_qz[z][x].i
-deltaxx[z-z_min][x]*((1-2*beta)*past_pxa[z][x].i
+beta*pxa[z][x].i)
-0.5*deltaxy[z-z_min][x]*((1-2*beta)*(past_pya[z][x].i+past_pya[z][x-1].i)+beta*(pya[z][x].i+pya[z][x-1].i))
-0.25*deltaxz[z-z_min][x]*((1-2*beta)*(past_pza[z][x].i+past_pza[z][x-1].i+past_pza[z+1][x].i+past_pza[z+1][x-1].i)
+beta*(pza[z][x].i+pza[z][x-1].i+pza[z+1][x].i+pza[z+1][x-1].i));
}
//Pya Update
ryy[z][x].r=pya[z][x].r
-0.5*cepy1[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].r+past_qy[z+1][x].r-past_qy[z][x+1].r-past_qy[z][x].r)
+cepy2[z-z_min][x-x_min].r*S*(past_qx[z+1][x].r-past_qx[z][x].r)
-cepy2[z-z_min][x-x_min].r*S*(past_qz[z][x+1].r-past_qz[z][x].r)
+0.5*cepy3[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].r+past_qy[z][x+1].r-past_qy[z+1][x].r-past_qy[z][x].r)
+0.25*sin(th)
*(deltayx[z-z_min][x]*cepx3[z-z_min][x].r+deltayy[z-z_min][x]*cepy3[z-z_min][x].r+deltayz[z-z_min][x]*cepz3[z-z_min][x].r)
*(past_qy[z][x].r+past_qy[z+1][x].r+past_qy[z][x+1].r+past_qy[z+1][x+1].r)
-0.5*sin(th)
*(deltayx[z-z_min][x]*cepx2[z-z_min][x].r+deltayy[z-z_min][x]*cepy2[z-z_min][x].r+deltayz[z-z_min][x]*cepz2[z-z_min][x].r)
*(past_qz[z][x].r+past_qz[z][x+1].r)
-0.5*deltayx[z-z_min][x]
*((1-2*beta)*(past_pxa[z][x].r+past_pxa[z][x+1].r)+beta*(pxa[z][x].r+pxa[z][x+1].r))
-deltayy[z-z_min][x]*((1-2*beta)*past_pya[z][x].r+beta*pya[z][x].r)
-0.5*deltayz[z-z_min][x]
*((1-2*beta)*(past_pza[z][x].r+past_pza[z+1][x].r)+beta*(pza[z][x].r+pza[z+1][x].r));
ryy[z][x].i=pya[z][x].i
-0.5*cepy1[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].i+past_qy[z+1][x].i-past_qy[z][x+1].i-past_qy[z][x].i)
+cepy2[z-z_min][x-x_min].r*S*(past_qx[z+1][x].i-past_qx[z][x].i)
-cepy2[z-z_min][x-x_min].r*S*(past_qz[z][x+1].i-past_qz[z][x].i)
+0.5*cepy3[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].i+past_qy[z][x+1].i-past_qy[z+1][x].i-past_qy[z][x].i)
+0.25*sin(th)
*(deltayx[z-z_min][x]*cepx3[z-z_min][x].r+deltayy[z-z_min][x]*cepy3[z-z_min][x].r+deltayz[z-z_min][x]*cepz3[z-z_min][x].r)
*(past_qy[z][x].i+past_qy[z+1][x].i+past_qy[z][x+1].i+past_qy[z+1][x+1].i)
-0.5*sin(th)
*(deltayx[z-z_min][x]*cepx2[z-z_min][x].r+deltayy[z-z_min][x]*cepy2[z-z_min][x].r+deltayz[z-z_min][x]*cepz2[z-z_min][x].r)
*(past_qz[z][x].i+past_qz[z][x+1].i)
-0.5*deltayx[z-z_min][x]
*((1-2*beta)*(past_pxa[z][x].i+past_pxa[z][x+1].i)+beta*(pxa[z][x].i+pxa[z][x+1].i))
-deltayy[z-z_min][x]*((1-2*beta)*past_pya[z][x].i+beta*pya[z][x].i)
-0.5*deltayz[z-z_min][x]
*((1-2*beta)*(past_pza[z][x].i+past_pza[z+1][x].i)+beta*(pza[z][x].i+pza[z+1][x].i));
//Pza Update
rzz[z][x].r=pza[z][x].r
-0.25*cepz1[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].r+past_qy[z+1][x].r-past_qy[z-1][x+1].r-past_qy[z-1][x].r)
+0.5*cepz2[z-z_min][x-x_min].r*S*(past_qx[z+1][x].r-past_qx[z-1][x].r)
-0.5*cepz2[z-z_min][x-x_min].r*S*(past_qz[z][x+1].r+past_qz[z-1][x+1].r-past_qz[z][x].r-past_qz[z-1][x].r)
+cepz3[z-z_min][x-x_min].r*S*(past_qy[z][x+1].r-past_qy[z][x].r)
+0.5*sin(th)
*(deltazx[z-z_min][x]*cepx3[z-z_min][x].r+deltazy[z-z_min][x]*cepy3[z-z_min][x].r+deltazz[z-z_min][x]*cepz3[z-z_min][x].r)
*(past_qy[z][x].r+past_qy[z][x+1].r)
-0.25*sin(th)
*(deltazx[z-z_min][x]*cepx2[z-z_min][x].r+deltazy[z-z_min][x]*cepy2[z-z_min][x].r+deltazz[z-z_min][x]*cepz2[z-z_min][x].r)
*(past_qz[z][x].r+past_qz[z-1][x].r+past_qz[z][x+1].r+past_qz[z-1][x+1].r)
-0.25*deltazx[z-z_min][x]
*((1-2*beta)*(past_pxa[z][x].r+past_pxa[z-1][x].r+past_pxa[z][x+1].r+past_pxa[z-1][x+1].r)
+beta*(pxa[z][x].r+pxa[z-1][x].r+pxa[z][x+1].r+pxa[z-1][x+1].r))
-0.5*deltazy[z-z_min][x]
*((1-2*beta)*(past_pya[z][x].r+past_pya[z-1][x].r)+beta*(pya[z][x].r+pya[z-1][x].r))
-deltazz[z-z_min][x]*((1-2*beta)*past_pza[z][x].r+beta*pza[z][x].r);
rzz[z][x].i=pza[z][x].i
-0.25*cepz1[z-z_min][x-x_min].r*S*(past_qy[z+1][x+1].i+past_qy[z+1][x].i-past_qy[z-1][x+1].i-past_qy[z-1][x].i)
+0.5*cepz2[z-z_min][x-x_min].r*S*(past_qx[z+1][x].i-past_qx[z-1][x].i)
-0.5*cepz2[z-z_min][x-x_min].r*S*(past_qz[z][x+1].i+past_qz[z-1][x+1].i-past_qz[z][x].i-past_qz[z-1][x].i)
+cepz3[z-z_min][x-x_min].r*S*(past_qy[z][x+1].i-past_qy[z][x].i)
+0.5*sin(th)
*(deltazx[z-z_min][x]*cepx3[z-z_min][x].r+deltazy[z-z_min][x]*cepy3[z-z_min][x].r+deltazz[z-z_min][x]*cepz3[z-z_min][x].r)
*(past_qy[z][x].i+past_qy[z][x+1].i)
-0.25*sin(th)
*(deltazx[z-z_min][x]*cepx2[z-z_min][x].r+deltazy[z-z_min][x]*cepy2[z-z_min][x].r+deltazz[z-z_min][x]*cepz2[z-z_min][x].r)
*(past_qz[z][x].i+past_qz[z-1][x].i+past_qz[z][x+1].i+past_qz[z-1][x+1].i)
-0.25*deltazx[z-z_min][x]
*((1-2*beta)*(past_pxa[z][x].i+past_pxa[z-1][x].i+past_pxa[z][x+1].i+past_pxa[z-1][x+1].i)
+beta*(pxa[z][x].i+pxa[z-1][x].i+pxa[z][x+1].i+pxa[z-1][x+1].i))
-0.5*deltazy[z-z_min][x]
*((1-2*beta)*(past_pya[z][x].i+past_pya[z-1][x].i)+beta*(pya[z][x].i+pya[z-1][x].i))
-deltazz[z-z_min][x]*((1-2*beta)*past_pza[z][x].i+beta*pza[z][x].i);
}
}
for (z=z_min;z<z_max;z++)
{
for(x=x_min;x<x_max;x++)
{
pxa[z][x].r=cxi11[z-z_min][x].r*rxx[z][x].r+cxi12[z-z_min][x].r*(ryy[z][x].r+ryy[z][x-1].r)/2.0+cxi13[z-z_min][x].r*(rzz[z][x].r+rzz[z+1][x].r+rzz[z][x-1].r+rzz[z+1][x-1].r)/4.0;
pxa[z][x].i=cxi11[z-z_min][x].r*rxx[z][x].i+cxi12[z-z_min][x].r*(ryy[z][x].i+ryy[z][x-1].i)/2.0+cxi13[z-z_min][x].r*(rzz[z][x].i+rzz[z+1][x].i+rzz[z][x-1].i+rzz[z+1][x-1].i)/4.0;
pya[z][x].r=cxi21[z-z_min][x].r*(rxx[z][x].r+rxx[z][x+1].r)/2.0+cxi22[z-z_min][x].r*ryy[z][x].r+cxi23[z-z_min][x].r*(rzz[z][x].r+rzz[z+1][x].r)/2.0;
pya[z][x].i=cxi21[z-z_min][x].r*(rxx[z][x].i+rxx[z][x+1].i)/2.0+cxi22[z-z_min][x].r*ryy[z][x].i+cxi23[z-z_min][x].r*(rzz[z][x].i+rzz[z+1][x].i)/2.0;
pza[z][x].r=cxi31[z-z_min][x].r*(rxx[z][x].r+rxx[z][x+1].r+rxx[z-1][x].r+rxx[z-1][x+1].r)/4.0+cxi32[z-z_min][x].r*(ryy[z][x].r+ryy[z-1][x].r)/2.0+cxi33[z-z_min][x].r*rzz[z][x].r;
pza[z][x].i=cxi31[z-z_min][x].r*(rxx[z][x].i+rxx[z][x+1].i+rxx[z-1][x].i+rxx[z-1][x+1].i)/4.0+cxi32[z-z_min][x].r*(ryy[z][x].i+ryy[z-1][x].i)/2.0+cxi33[z-z_min][x].r*rzz[z][x].i;
}
}
}
Replies
: