📄 emap4.c
字号:
complex *XVec,*AVec,**RIDat;
{
int II,JJ,KK;
for(II=0;II<Size;II++)
{
AVec[II]=COMplex_Null();
}
if(S==' ') /* Multiply matrix with vector */
{
for(II=0;II<Size;II++)
{
for(KK=0;KK<RIIndex[II];KK++)
{
JJ = RICol[II][KK];
AVec[II] = COMplex_Add(AVec[II],COMplex_Mul(RIDat[II][KK],XVec[JJ]));
}
}
}
else /* Multiply hermitian of matrix with vector */
/* Here the matrix is symmetric, so we don't need to transpose */
/* Only conjugate */
{
for(II=0;II<Size;II++)
{
for(KK=0;KK<RIIndex[II];KK++)
{
JJ = RICol[II][KK];
AVec[II] = COMplex_Add(AVec[II],COMplex_Mul(
COMplex_Conjg(RIDat[II][KK]),XVec[JJ]));
}
}
}
}
/****************************Main routine************************************/
void ConjugateSolver(MatrixSize,RowIIndex,RowICol,RowIDat,RHSVec,ITmax,TOL)
int ITmax,MatrixSize,**RowICol,*RowIIndex;
complex *RHSVec,**RowIDat;
double TOL;
{
int II,JJ,IterIndex;
double Residue,Vnrm;
printf("\n %d ",freq_step);
if(freq_step == 0) /* Only for the first frequency step */
{
X = CMPLX_Vector(MatrixSize);
P = CMPLX_Vector(MatrixSize);
PP = CMPLX_Vector(MatrixSize);
R = CMPLX_Vector(MatrixSize);
RR = CMPLX_Vector(MatrixSize);
A_P = CMPLX_Vector(MatrixSize);
AH_PP = CMPLX_Vector(MatrixSize);
}
for(II=0;II<MatrixSize;II++)
{
X[II] = COMplex_Null();
P[II] = COMplex_Null();
PP[II] = COMplex_Null();
R[II] = COMplex_Null();
RR[II] = COMplex_Null();
A_P[II] = COMplex_Null();
AH_PP[II] = COMplex_Null();
}
/* for debugging */
/* for (II=0;II<MatrixSize;II++) */
/* { */
/* for(JJ=0;JJ<RowIIndex[II];JJ++) */
/* { */
/*printf("%d %f %f \n",RowICol[II][JJ],RowIDat[II][JJ].x,RowIDat[II][JJ].y); */
/* } */
/* printf("\n %d %d\n",InnerEdgeStat[II], II); */
/* } */
/* Initialization */
MatrixVectorProd(' ',MatrixSize,X,R,RowIIndex,RowICol,RowIDat);
for(II=0;II<MatrixSize;II++)
{
R[II] = COMplex_Sub(RHSVec[II],R[II]);
P[II] = R[II];
RR[II] = COMplex_Conjg(R[II]);
PP[II] = RR[II];
}
/*for(II=0;II<MatrixSize;II++)
{
printf(" R = %f %f",R[II].x,R[II].y);
printf(" RR = %f %f \n ",RR[II].x,RR[II].y);
printf(" P = %f %f",P[II].x,P[II].y);
printf(" PP = %f %f \n",PP[II].x,PP[II].y);
}
*/
Vnrm = Vnorm(RHSVec,MatrixSize);
IterIndex = 0;
/* Actual Iteration */
for(;;)
{
MatrixVectorProd(' ',MatrixSize,P,A_P,RowIIndex,RowICol,RowIDat);
for(II=0;II<MatrixSize;II++)
{
AH_PP[II] = COMplex_Conjg(A_P[II]);
}
Alpha = Inner_Prod(RR,R,MatrixSize);
temp_var = Inner_Prod(PP,A_P,MatrixSize);
Alpha = COMplex_Div(Alpha,temp_var);
/* Alpha is the step length parameter */
for(II=0;II<MatrixSize;II++)
{
/* New Solution Estimate */
X[II] = COMplex_Add(X[II],COMplex_Mul(Alpha,P[II]));
R[II] = COMplex_Sub(R[II],COMplex_Mul(Alpha,A_P[II])); /* Residual */
RR[II] = COMplex_Sub(RR[II],COMplex_Mul(COMplex_Conjg(Alpha),AH_PP[II]));
/* Bi-residual */
}
Beta = Inner_Prod(AH_PP,R,MatrixSize);
Beta = COMplex_Div(Beta,temp_var);
Beta = Real_Mul(-1.0,Beta); /* Beta is the Bi-Conjugacy coefficient */
for(II=0;II<MatrixSize;II++)
{
P[II] = COMplex_Add(R[II],COMplex_Mul(Beta,P[II])); /* Direction */
PP[II] = COMplex_Add(RR[II],COMplex_Mul(COMplex_Conjg(Beta),PP[II]));
/* Bi-Direction */
}
IterIndex = IterIndex + 1; /* Next iteration */
Residue = sqrt(Vnorm(R,MatrixSize)/Vnrm);
/* printf("Iteration = %d Residue = %e TOL = %e\n",IterIndex,Residue,TOL); */
/* fflush(stdout); */
if(Residue<=TOL) /* Test termination condition */
{
printf("Convergence achieved\n");
printf("Iteration = %d Residue = %e TOL = %e\n",IterIndex,Residue,TOL);
fflush(stdout);
break;
}
else if(IterIndex==ITmax) /* No. of iterations has exceeded the limit */
{
printf("Iteration Exceeds\n");
fflush(stdout);
break;
}
else
{
continue ;
}
}/* End of for loop*/
for(II=0;II<MatrixSize;II++)
{
RHSVec[II] = X[II];
Old_Solution[II] = X[II]; /* Starting point for next iteration */
}
}
/***************************************************************************
* FUNCTION : Read_Input_Pass_1()
****************************************************************************/
void Read_Input_Pass_1()
{
int x1, y1, z1, x2, y2, z2, p, p1, p2, OutF_Count=0, output_flag=0;
char type[20], att1[20], att2[18], att3[18], att4[18];
float tmp1, tmp2;
double freq = 0;
char buffer[80];
Min_X = 10000; Min_Y=10000; Min_Z=10000;
Max_X = 0; Max_Y=0; Max_Z=0;
while (fgets(buffer, 80, InF))
{
if (buffer[0] == '#') continue; /* Comment statement */
if(!strncmp(buffer, "freqstep",8)) /* freqstep keyword */
{
sscanf(buffer, "%s%d%s",type, &NUM_OF_STEPS, att1);
FREQ_INC = atof(att1); /* Convert from string to float */
printf("N0. of steps%d inc%f\n",NUM_OF_STEPS,FREQ_INC);
continue;
}
if(!strncmp(buffer, "default_output",14))/*keyword for default output*/
{
sscanf(buffer, "%s%s", type, Out_FileName0);
OutF_0=fopen(Out_FileName0, "w");
if (OutF_0 == NULL)
{
fprintf(stderr,
" Error: Can't create default output file %s\n",
Out_FileName0);
exit(1);
}
output_flag=1;
continue;
}
if(!strncmp(buffer, "celldim", 7)) { /* celldim keyword */
continue;
}
if(sscanf(buffer, "%s%d%d%d%d%d%d%s%s%s%s",
type, &x1, &y1, &z1, &x2, &y2, &z2,
att1, att2, att3, att4)>2)
{
if(!strcmp(type,"boundary")) {
fprintf(stderr,
"\n\nERROR: boundary keyword is not allowed by EMAP4 \n");
exit(1);}
if(!strcmp(type, "dielectric"))/* dielectric block */
{
if (x1<Min_X) Min_X=x1;
if (x2<Min_X) Min_X=x2;
if (x1>Max_X) Max_X=x1;
if (x2>Max_X) Max_X=x2;
if (y1<Min_Y) Min_Y=y1;
if (y2<Min_Y) Min_Y=y2;
if (y1>Max_Y) Max_Y=y1;
if (y2>Max_Y) Max_Y=y2;
if (z1<Min_Z) Min_Z=z1;
if (z2<Min_Z) Min_Z=z2;
if (z1>Max_Z) Max_Z=z1;
if (z2>Max_Z) Max_Z=z2;
if((x1==x2)||(y1==y2)||(z1==z2))
{fprintf(stderr, "ERROR: diel thickness is not finite./n");
exit(1); } continue;
}
if(!strcmp(type, "PML")) /* Perfectly matched layer */
{
if (x1<Min_X) Min_X=x1;
if (x2<Min_X) Min_X=x2;
if (x1>Max_X) Max_X=x1;
if (x2>Max_X) Max_X=x2;
if (y1<Min_Y) Min_Y=y1;
if (y2<Min_Y) Min_Y=y2;
if (y1>Max_Y) Max_Y=y1;
if (y2>Max_Y) Max_Y=y2;
if (z1<Min_Z) Min_Z=z1;
if (z2<Min_Z) Min_Z=z2;
if (z1>Max_Z) Max_Z=z1;
if (z2>Max_Z) Max_Z=z2;
if((x1==x2)||(y1==y2)||(z1==z2))
{
fprintf(stderr, "ERROR: PML thickness is not finite./n");
exit(1);
} continue;
}
if(!strcmp(type, "box")) /* rectangular conducting box */
{
if((x1==x2)||(y1==y2)||(z1==z2)){
fprintf(stderr, "ERROR: box dimension is invalid.");
exit(1);}
if (x1<Min_X) Min_X=x1;
if (x2<Min_X) Min_X=x2;
if (x1>Max_X) Max_X=x1;
if (x2>Max_X) Max_X=x2;
if (y1<Min_Y) Min_Y=y1;
if (y2<Min_Y) Min_Y=y2;
if (y1>Max_Y) Max_Y=y1;
if (y2>Max_Y) Max_Y=y2;
if (z1<Min_Z) Min_Z=z1;
if (z2<Min_Z) Min_Z=z2;
if (z1>Max_Z) Max_Z=z1;
if (z2>Max_Z) Max_Z=z2;
continue;
}
if(!strcmp(type,"conductor")) /* PEC */
{
if (x1<Min_X) Min_X=x1;
if (x2<Min_X) Min_X=x2;
if (x1>Max_X) Max_X=x1;
if (x2>Max_X) Max_X=x2;
if (y1<Min_Y) Min_Y=y1;
if (y2<Min_Y) Min_Y=y2;
if (y1>Max_Y) Max_Y=y1;
if (y2>Max_Y) Max_Y=y2;
if (z1<Min_Z) Min_Z=z1;
if (z2<Min_Z) Min_Z=z2;
if (z1>Max_Z) Max_Z=z1;
if (z2>Max_Z) Max_Z=z2;
continue;
}
if(!strcmp(type,"resistor")) /* resistor */
{
if (x1<Min_X) Min_X=x1;
if (x2<Min_X) Min_X=x2;
if (x1>Max_X) Max_X=x1;
if (x2>Max_X) Max_X=x2;
if (y1<Min_Y) Min_Y=y1;
if (y2<Min_Y) Min_Y=y2;
if (y1>Max_Y) Max_Y=y1;
if (y2>Max_Y) Max_Y=y2;
if (z1<Min_Z) Min_Z=z1;
if (z2<Min_Z) Min_Z=z2;
if (z1>Max_Z) Max_Z=z1;
if (z2>Max_Z) Max_Z=z2;
continue;
}
if(!strcmp(type, "aperture")) /* aperture */
{
if((x1!=x2) & (y1!=y2) & (z1!=z2))
{
fprintf(stderr,"ERROR: aperture dimension is invalid.\n");
exit(1);
}
if(((x1==x2)&&(y1==y2))||((x1==x2)&&(z1==z2))
||((z1==z2)&&(y1==y2)))
{
fprintf(stderr, "ERROR: aperture dimension is invalid.\n");
exit(1); }
continue;
}
if(!strcmp(type, "seam")) continue; /* Not available */
if(!strcmp(type,"esource")) /* electric field source */
{
if (x1<Min_X) Min_X=x1;
if (x2<Min_X) Min_X=x2;
if (x1>Max_X) Max_X=x1;
if (x2>Max_X) Max_X=x2;
if (y1<Min_Y) Min_Y=y1;
if (y2<Min_Y) Min_Y=y2;
if (y1>Max_Y⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -