1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

MMUL(K1, tmp1,tmp2);
//% K1 = (R1  dy^2/48*iR2)\(eye(4*K) + i*(32*sqrt(3))/12*dy*iR2);
#pragma omp parallel for private(ii,jj)
for (ii=1;ii<=4*K;ii++)
{
for (jj=1;jj<=4*K;jj++)
{
tmp2.ld(ii,jj, cmplx(0,(3.f+2.f*sqrt(3.f))/12.f*dy)* K1(ii,jj));
} //end
tmp2.ld(ii,ii, tmp2(ii,ii)+(1.f)); // %eye
} //end
//%K2 = iR2 *(eye(4*K) + i*(3+2*sqrt(3))/12*dy* K1);
//%K2 = iR2*tmp2;
MMUL(K2, iR2,tmp2);
#pragma omp parallel for private(ii,jj)
for (ii=1;ii<=4*K;ii++)
{
for (jj=1;jj<=4*K;jj++)
{
tmp2.ld(ii,jj, cmplx(0,dy/2)*(K1(ii,jj)+K2(ii,jj)));
} //end
tmp2.ld(ii,ii, tmp2(ii,ii)+(1.f));// %eye
} //end
//%Fnm = (eye(4*K) + i*dy/2*(K1+K2))*Fnm;
MMUL(tmp1, tmp2,Fnm);
Fnm=tmp1;
} //end
} //end
 