i write codes of MRT-D3Q19 with c++,but find it can’t gain th correct solution, i check it for many times, but don’t find where is error. I change it to SRT-D3Q19,and gain good solution. by the way, what is formula of the square duct flow, how to understand.
thank you
fllowing:
#include
#include
#include
#include
#include
#include
#include
#include
using namespace std;
const int Q = 19;
const int NX = 11;
const int NY = 11;
const int NZ = 21;
const double U = 0.1;
double cr=25;
int e[Q][3] = {{0,0,0},{1,0,0},{-1,0,0},{0,1,0},{0,-1,0},{0,0,1},{0,0,-1},{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0},{1,0,1},{-1,0,1},{1,0,-1},{-1,0,-1},{0,1,1},{0,1,-1},{0,-1,1},{0,-1,-1}};
double w[Q] = {1.0/3,1.0/18,1.0/18,1.0/18,1.0/18,1.0/18,1.0/18,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36,1.0/36};
//
double rho[NX][NY][NZ],u[NX][NY][NZ][3],u0[NX][NX][NZ][3];
//
double f[NX][NY][NZ][Q],feqq[NX][NY][NZ][Q],ftemp[NX][NY][NZ][Q];
//
double fmom[NX][NY][NZ][Q],fmeq[NX][NY][NZ][Q],force[NX][NY][NZ][Q];
int opposxte[Q]={0,2,1,4,3,6,5,10,9,8,7,14,13,12,11,18,17,16,15};
int image[NX][NY][NZ];
int x,y,z,k,l,m,n,g;
double xy[NY];
double Re,rho0,tau,nu,error1,error2,omega,cx,cy,cz,suma;
double Fx,Fy,Fz;
double s[Q]; //D3Q19 relaxatxon ratxo matrxx
double tm[19][19] = {{1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0},
{-30.0,-11.0,-11.0,-11.0,-11.0,-11.0,-11.0,8.0,8.0,8.0,8.0,8.0,8.0,8.0,8.0,8.0,8.0,8.0,8.0},
{12.0,-4.0,-4.0,-4.0,-4.0,-4.0,-4.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0},
{0.0,1.0,-1.0,0.0,0.0,0.0,0.0,1.0,-1.0,1.0,-1.0,1.0,-1.0,1.0,-1.0,0,0,0,0},
{0.0,-4.0,4.0,0.0,0.0,0.0,0.0,1.0,-1.0,1.0,-1.0,1.0,-1.0,1.0,-1.0,0,0,0,0},
{0.0,0.0,0.0,1.0,-1.0,0.0,0.0,1.0,1.0,-1.0,-1.0,0,0,0,0,1.0,-1.0,1.0,-1.0},
{0.0,0.0,0.0,-4.0,4.0,0.0,0.0,1.0,1.0,-1.0,-1.0,0,0,0,0,1.0,-1.0,1.0,-1.0},
{0.0,0.0,0.0,0.0,0.0,1.0,-1.0,0.0,0.0,0.0,0.0,1.0,1.0,-1.0,-1.0,1.0,1.0,-1.0,-1.0},
{0.0,0.0,0.0,0.0,0.0,-4.0,4.0,0.0,0.0,0.0,0.0,1.0,1.0,-1.0,-1.0,1.0,1.0,-1.0,-1.0},
{0.0,2.0,2.0,-1.0,-1.0,-1.0,-1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,-2.0,-2.0,-2.0,-2.0},
{0.0,-4.0,-4.0,2.0,2.0,2.0,2.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,-2.0,-2.0,-2.0,-2.0},
{0.0,0.0,0.0,1.0,1.0,-1.0,-1.0,1.0,1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,-2.0,-2.0,2.0,2.0,1.0,1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,-1.0,-1.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0-1.0,-1.0,1.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,-1.0,-1.0,1.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,-1.0,1.0,-1.0,-1.0,1.0,-1.0,1.0,0.0,0.0,0.0,0.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,-1.0,-1.0,1.0,1.0,0.0,0.0,0.0,0.0,1.0,-1.0,1.0,-1.0},
{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0,1.0,1.0}};
double tminv[19][19],tme[19][19];//D3Q19 transver matrxx
//
double inrho = 1.1;
double outrho = 1.0;
double deltaP;
double umax;
double ua[NX][NY][NZ];
double jx,jy,jz;
//
void imagec();
void init();
void gaussy();
double feqc(int k,double rho,double u[3]);
void evolution();
void outputdata();
void Error();
int main()
{
using namespace std;
init();
imagec();
gaussy();
for(x=0; x<19; x++)
for(y=0; y<19; y++)
{
tme[x][y]=0;
for(z=0; z<19; z++)
{
tme[x][y]=tme[x][y]+tm[x][z]*tminv[z][y];
}
}
//
for(n=0; ;n++)
{
evolution();
if(n%1000 == 0)
{
Error();
cout<<“it times”<<n<<" error1="<<error1<<" error2="<<error2<<endl;
if(n >= 100)
{
if(n%1000== 0)
outputdata();
if (error1<1.0E-14)
break;
}
}
}
return 0;
}
//
//
void init()
{
rho0 = 1.0;
tau = 0.515;
omega=1.0/tau;
s[0] = 0.0; //free choice
s[1] = 1.0/tau;
s[2] = 1.0/tau;
s[3] = 0.0; // free choice
s[4] = 1.0/tau;
s[5] = 0.0; // free choice
s[6] =s[4]; //
s[7] =0.0; // free choice
s[8]= s[4];
s[9]= 1.0/tau;
s[10]= 1.0/tau;
s[11]=s[9];
s[12]=s[10];
s[13]=1.0/tau;
s[14] = s[13];
s[15] = s[13];
s[16] = 1.0/tau;
s[17] = s[16];
s[18] = s[16];
nu = 1.0/3.0*(tau-0.5);
cx=0.2NX;
cy=0.5NY;
std::cout<<“tau=”<<tau<<endl;
std::cout<<“nu=”<<nu<<endl;
for(x = 0; x <NX ; x++)
for(y = 0; y <NY ; y++)
for(z = 0; z <NZ ; z++)
{
u[x][y][z][0] = 0;
u[x][y][z][1] = 0;
u[x][y][z][2] = 0;
rho[x][y][z] = rho0;
// for(k = 0;k < Q;k++)
// {
// f[x][y][z][k]=feqc(k,rho[x][y][z],u[x][y][z]);
// }
f[x][y][z][0]=1.0/3.0;
f[x][y][z][1]=1.0/18.0;
f[x][y][z][2]=1.0/18.0;
f[x][y][z][3]=1.0/18.0;
f[x][y][z][4]=1.0/18.0;
f[x][y][z][5]=1.0/18.0;
f[x][y][z][6]=1.0/18.0;
f[x][y][z][7]=1.0/36.0;
f[x][y][z][8]=1.0/36.0;
f[x][y][z][9]=1.0/36.0;
f[x][y][z][10]=1.0/36.0;
f[x][y][z][11]=1.0/36.0;
f[x][y][z][12]=1.0/36.0;
f[x][y][z][13]=1.0/36.0;
f[x][y][z][14]=1.0/36.0;
f[x][y][z][15]=1.0/36.0;
f[x][y][z][16]=1.0/36.0;
f[x][y][z][17]=1.0/36.0;
f[x][y][z][18]=1.0/36.0;
}
//
deltaP =6.25000000E-05; //1.0/3.0*(inrho-outrho)/(NZ-1);
std::cout<<“deltaP=”<<deltaP<<endl;
// system(“pause”);
umax = deltaP/(16nu)(NY-2)*(NY-2);
std::cout<<“umax=”<<umax<<endl;
system(“pause”);
for(y=0;y<=(NY-1)/2;y++)
{
xy[y]=(NY-1)/2.0-y;
}
//
for(y=(NY-1)/2;y<NY;y++)
{
xy[y]=(NY-1)/2.0-y;
}
//
xy[0]=(NY-2)/2.0;
xy[NY-1]=-(NY-2)/2.0;
for(y=0;y<NY;y++)
{
std::cout<<“xy[y]=”<<xy[y]<<endl;
}
// system(“pause”);
for(x = 0; x <NX ; x++)
for(y = 0; y <NY ; y++)
for(z = 0; z <NZ ; z++)
{
if(image[x][y][z]==0)
{
ua[x][y][z] =deltaP/(4nu)((NY-2)*(NY-2)/4.0-xy[y]*xy[y]);
}
}
}
//
//
void imagec()
{
for(x = 0; x <NX ; x++)
for(y = 0; y <NY ; y++)
for(z = 0; z <NZ ; z++)
{
// xf (((x-cx)(x-cx)+(y-cy)(y-cy))<=cr*cr)
// xf(x==1&&x==NX)
// {
// xmage[x][y]=1;
// }
// else
// {
image[x][y][z]=0;
// }
}
// Bottom plate
for(x = 0; x < NX ; x++)
for(z = 0; z < NZ ; z++)
{
image[x][0][z]=1;
}
// Top plate
for(x = 0; x < NX ; x++)
for(z = 0; z < NZ ; z++)
{
image[x][NY-1][z]=1;
}
// Left plate
for(y = 0; y < NY ; y++)
for(z = 0; z < NZ ; z++)
{
image[0][y][z]=1;
}
// Rxgzt plate
for(y = 0; y < NY ; y++)
for(z = 0; z < NZ ; z++)
{
image[NX-1][y][z]=1;
}
// ofstream xmag;
// xmag.open
}
//
//
double feqc(int k,double rho,double u[3])
{
double eu,uv,feq;
eu = (e[k][0]*u[0] + e[k][1]*u[1] +e[k][2]u[2] );
uv = (u[0]u[0] + u[1]u[1] + u[2]u[2]);
feq = w[k]rho(1.0 + 3.0eu +4.5eueu - 1.5uv);
return feq;
}
//
void evolution()
{
int xe,yn,zt,zb,ys,xw;
double temp,suma1,suma2;
//
//
for(x = 0; x <NX; x++)
for(y = 0; y <NY; y++)
for(z = 0; z <NZ; z++)
{
u0[x][y][z][0] = u[x][y][z][0];
u0[x][y][z][1] = u[x][y][z][1];
u0[x][y][z][2] = u[x][y][z][2];
rho[x][y][z] = 0;
u[x][y][z][0] = 0;
u[x][y][z][1] = 0;
u[x][y][z][2] = 0;
if(image[x][y][z]==0)
{//
for(k = 0;k < Q;k++)
{
rho[x][y][z] +=f[x][y][z][k];
u[x][y][z][0] += e[k][0]*f[x][y][z][k];
u[x][y][z][1] += e[k][1]*f[x][y][z][k];
u[x][y][z][2] += e[k][2]*f[x][y][z][k];
}
u[x][y][z][0]= u[x][y][z][0]/rho[x][y][z];//+0.5*0.0/rho[x][y][z];
u[x][y][z][1]= u[x][y][z][1]/rho[x][y][z];//+0.5*0.0/rho[x][y][z];
u[x][y][z][2]= u[x][y][z][2]/rho[x][y][z];//+0.5*deltaP/rho[x][y][z];
}//
}
for(x = 0; x <NX ; x++)
for(y = 0; y <NY ; y++)
for(z = 0; z <NZ ; z++)
for(k = 0;k < Q;k++)
{
if(image[x][y][z]==0)
{
feqq[x][y][z][k]=feqc(k,rho[x][y][z],u[x][y][z]);
}
}
// computxng the velocxty moments
for(x = 0; x <NX ; x++)
for(y = 0; y <NY ; y++)
for(z = 0; z <NZ ; z++)
{
jx=u[x][y][z][0];
jy=u[x][y][z][1];
jz=u[x][y][z][2];
fmeq[x][y][z][0] = rho[x][y][z];
fmeq[x][y][z][1] = -11.0rho[x][y][z] + 19.0/rho0(jxjx+jyjy+jzjz);//e(eq)
fmeq[x][y][z][2] = 3.0rho[x][y][z] -11.0/(2.0rho0)(jxjx+jyjy+jzjz); //energy squareE(eq)
fmeq[x][y][z][3] = jx; //momentum xn x-dxr jx
fmeq[x][y][z][4] = -2.0/3.0jx; //energy flux xn x-dxr—/-jx qx
fmeq[x][y][z][5] = jy; //momentum xn y-dxr yy
fmeq[x][y][z][6] = -2.0/3.0jy; //energy flux xn y-dxr—/-jy qy
fmeq[x][y][z][7] = jz; //momentum xn z-dxr yz
fmeq[x][y][z][8] = -2.0/3.0jz; //energy flux xn z-dxr—/-jz qz
fmeq[x][y][z][9] = 1.0/rho0*(2.0jxjx-jyjy-jzjz); // pxx dxagonal comp of stress tensor
fmeq[x][y][z][10] = -1.0/(2.0rho0)(2.0jxjx-jyjy-jzjz); //paixx
fmeq[x][y][z][11] = 1.0/rho0*(jyjy-jzjz); //pww
fmeq[x][y][z][12] = -1.0/(2.0rho0)(jyjy-jzjz); //paiww
fmeq[x][y][z][13] = 1.0/rho0jxjy; //pxy
fmeq[x][y][z][14] = 1.0/rho0jyjz; //pyz
fmeq[x][y][z][15] = 1.0/rho0jxjz; //pxz
fmeq[x][y][z][16] = 0.0;
fmeq[x][y][z][17] = 0.0;
fmeq[x][y][z][18] = 0.0;
}
// computxng tze external body force
for(x = 0; x <NX ; x++)
for(y = 0; y <NY ; y++)
for(z = 0; z <NZ ; z++)
{
Fx = 0.0;
Fy = 0.0;
Fz = deltaP;
force[x][y][z][0] = 0;
force[x][y][z][1] = 38.0*(u[x][y][z][0]Fx+u[x][y][z][1]Fy+u[x][y][z][2]Fz);
force[x][y][z][2] = -11.0(u[x][y][z][0]Fx+u[x][y][z][1]Fy+u[x][y][z][2]Fz);
force[x][y][z][3] = Fx;
force[x][y][z][4] = -2.0/3.0Fx;
force[x][y][z][5] = Fy;
force[x][y][z][6] = -2.0/3.0Fy;
force[x][y][z][7]= Fz;
force[x][y][z][8]= -2.0/3.0Fz;
force[x][y][z][9]= 2.0(2.0u[x][y][z][0]*Fx-u[x][y][z][1]*Fy-u[x][y][z][2]Fz);
force[x][y][z][10]= -(2.0u[x][y][z][0]*Fx-u[x][y][z][1]*Fy-u[x][y][z][2]Fz);
force[x][y][z][11] = 2.0(u[x][y][z][1]*Fy-u[x][y][z][2]*Fz);
force[x][y][z][12] = -(u[x][y][z][1]*Fy-u[x][y][z][2]*Fz);
force[x][y][z][13] = u[x][y][z][0]*Fy+u[x][y][z][1]*Fx;
force[x][y][z][14] = u[x][y][z][1]*Fz+u[x][y][z][2]*Fy;
force[x][y][z][15] = u[x][y][z][2]*Fx+u[x][y][z][0]*Fz;
force[x][y][z][16] = 0.0;
force[x][y][z][17] = 0.0;
force[x][y][z][18] = 0.0;
}
//collxsxon start
//f->M
for(x=0; x<NX; x++)
for(y=0; y<NY; y++)
for(z=0; z<NZ; z++)
for(k=0;k<Q;k++)
{
suma1 = 0.0;
for(g=0;g<Q;g++)
{
suma1=suma1+tm[k][g]*f[x][y][z][g];
}
fmom[x][y][z][k]=suma1;
}
//
for(x=0;x<NX;x++)
for(y=0;y<NY;y++)
for(z=0;z<NZ;z++)
for(k=0;k<Q;k++)
{
if(image[x][y][z]==0)
{
// f[x][y][z][k] =f[x][y][z][k]+ omega*(feqq[x][y][z][k]-f[x][y][z][k])+3.0w[k](e[k][0]*0.0+e[k][1]0.0+e[k][2]deltaP);
fmom[x][y][z][k] =fmom[x][y][z][k]+ s[k](fmeq[x][y][z][k]-fmom[x][y][z][k]);//+(1-s[k]/2.0) force[x][y][z][k];
}
}
//M->f
for(x=0;x<NX;x++)
for(y=0;y<NY;y++)
for(z=0;z<NZ;z++)
for(k=0;k<Q;k++)
{
suma2=0;
for(g=0;g<Q;g++)
{
suma2=suma2+tminv[k][g]fmom[x][y][z][g];
}
f[x][y][z][k]=suma2+3.0w[k]*(e[k][0]*0.0+e[k][1]*0.0+e[k][2]*deltaP);
}
// Standard bounceback
for(x = 0; x<NX; x++)
for(y = 0; y<NY; y++)
for(z = 0; z<NZ; z++)
{
if(image[x][y][z]==1)
{
// Half-Way bounce back
temp = f[x][y][z][1];
f[x][y][z][1] = f[x][y][z][2];
f[x][y][z][2] = temp;
temp = f[x][y][z][3];
f[x][y][z][3] = f[x][y][z][4];
f[x][y][z][4] = temp;
temp = f[x][y][z][7];
f[x][y][z][7] = f[x][y][z][10];
f[x][y][z][10] = temp;
temp = f[x][y][z][9];
f[x][y][z][9] = f[x][y][z][8];
f[x][y][z][8] = temp;
temp = f[x][y][z][11];
f[x][y][z][11] = f[x][y][z][14];
f[x][y][z][14] = temp;
temp = f[x][y][z][13];
f[x][y][z][13] = f[x][y][z][12];
f[x][y][z][12] = temp;
temp = f[x][y][z][15];
f[x][y][z][15] = f[x][y][z][18];
f[x][y][z][18] = temp;
temp = f[x][y][z][16];
f[x][y][z][16] = f[x][y][z][17];
f[x][y][z][17] = temp;
temp = f[x][y][z][5];
f[x][y][z][5] = f[x][y][z][6];
f[x][y][z][6] = temp;
}
}
//collision end
//streaming start
for(x=0;x<NX;x++)
{
xe=(x<NX-1)?(x+1):(0 );
xw=(x>0)?(x-1):(NX-1);
for(y=0;y<NY;y++)
{
yn=(y<NY-1)?(y+1):(0 );
ys=(y>0)?(y-1):(NY-1);
for(z = 0; z<NZ;z++)
{
zt=(z<NZ-1)?(z+1):(0 );
zb=(z>0)?(z-1):(NZ-1);
ftemp[x][y][z][0]= f[x][y][z][0];
ftemp[xe][y][z][1]= f[x][y][z][1];
ftemp[xw][y][z][2]= f[x][y][z][2];
ftemp[x][yn][z][3]= f[x][y][z][3];
ftemp[x][ys][z][4]= f[x][y][z][4];
ftemp[x][y][zt][5]= f[x][y][z][5];
ftemp[x][y][zb][6]= f[x][y][z][6];
ftemp[xe][yn][z][7]= f[x][y][z][7];
ftemp[xe][ys][z][8]= f[x][y][z][8];
ftemp[xw][yn][z][9]= f[x][y][z][9];
ftemp[xw][ys][z][10]= f[x][y][z][10];
ftemp[xe][y][zt][11]= f[x][y][z][11];
ftemp[xw][y][zt][12]= f[x][y][z][12];
ftemp[xe][y][zb][13]= f[x][y][z][13];
ftemp[xw][y][zb][14]= f[x][y][z][14];
ftemp[x][yn][zt][15]= f[x][y][z][15];
ftemp[x][yn][zb][16]= f[x][y][z][16];
ftemp[x][ys][zt][17]= f[x][y][z][17];
ftemp[x][ys][zb][18]= f[x][y][z][18];
}
}
}
for(x=0;x<NX;x++)
for(y=0;y<NY;y++)
for(z=0;z<NZ;z++)
for(k=0;k<Q;k++)
{
f[x][y][z][k]= ftemp[x][y][z][k];
}
// streamxng end
}
//
void Error()
{
double temp1,temp2;
temp1 = 0;
temp2 = 0;
for(x = 0;x<NX;x++)
for(y = 0;y<NY;y++)
for(z = 0;z<NZ;z++)
{
temp1 += ((u[x][y][z][0] - u0[x][y][z][0])*(u[x][y][z][0] - u0[x][y][z][0]) +
(u[x][y][z][1] - u0[x][y][z][1])*(u[x][y][z][1] - u0[x][y][z][1])+
(u[x][y][z][2] - u0[x][y][z][2])*(u[x][y][z][2] - u0[x][y][z][2]));
temp2 += (u[x][y][z][0]*u[x][y][z][0] + u[x][y][z][1]*u[x][y][z][1]+u[x][y][z][2]*u[x][y][z][2]);
}
temp1 = sqrt(temp1);
temp2 = sqrt(temp2);
error1 = temp1/(temp2 + 1e-30);
//------------------------------------------------------------------
for(x = 0;x<NX;x++)
for(y = 0;y<NY;y++)
for(z = 0;z<NZ;z++)
{
temp1 += (u[x][y][z][2] - ua[x][y][z])*(u[x][y][z][2] - ua[x][y][z]);
temp2 += ua[x][y][z]*ua[x][y][z];
}
temp1 = sqrt(temp1);
temp2 = sqrt(temp2);
error2 = temp1/(temp2);
}
void outputdata()
{
if(n%100==0)
{
ofstream file;
char temp[10];
string str;
itoa(n,temp,10);
str = string(temp);
str += “u.dat”;
file.open(str.c_str());
file<<“Title=“LBM the Poiseillue Flow”\n”<<“VARxABLES=“Z”,“X”,“Y”,“u[0]”,“u[1]”,“u[2]”,“rho”\n”<<“ZONE T=“BOX”,z=”
<<NZ+1<<",x="<<NX+1<<",y="<<NY+1<<",F=POxNT"<<endl;
if (file.is_open())
{
for(z = 0; z<NZ; z++)
for(x = 0; x<NX; x++ )
for(y = 0; y<NY; y++ )
{
file<<setiosflags(ios::right)<<double(z)<<" "<<double(x)<<" "<<double(y)<<" "<<setiosflags(ios::scientific)<<setprecision(6)<<u[x][y][z][0]<<" "<<u[x][y][z][1]<<" "<<u[x][y][z][2]<<" "<<ua[x][y][z]<<" "<<rho[x][y][z]<<resetiosflags(ios::scientific)<<endl;
}
}
}
}
//
//采用部分主元法的高斯消去法求方阵A的逆矩阵B
//
//A 方阵 (xN)
//n 方阵的阶数 (xN)
//B 方阵 (OUT)
//返回true 说明正确计算出逆矩阵
//返回false说明矩阵A不存在逆矩阵
void gaussy()
{
int x,y,k;
double lMax,temp;
//临时矩阵存放A
double tt[19][19];
for(x=0;x<19;x++)
{
for(y=0;y<19;y++)
{
tt[x][y] = tm[x][y];
}
}
//初始化B为单位阵
for(x=0;x<19;x++)
{
for(y=0;y<19;y++)
{
if(x!=y)tminv[x][y] = 0;
else tminv[x][y] = 1;
}
}
for(x=0;x<19;x++)
{
//寻找主元
lMax = tt[x][x];
k = x;
for(y=x+1;y<19;y++) //扫描从x+1到n的各行
{
if( fabs(tt[y][x]) > fabs(lMax))
{
lMax = tt[y][x];
k = y;
}
}
//如果主元所在行不是第x行，进行行交换
if(k!=x)
{
for(y=0;y<19;y++)
{
temp = tt[x][y] ;
tt[x][y] = tt[k][y];
tt[k][y] = temp;
//B伴随计算
temp = tminv[x][y];
tminv[x][y] = tminv[k][y];
tminv[k][y] = temp;
}
}
//判断主元是否是0，如果是，则矩阵A不是满秩矩阵，不存在逆矩阵
if(tt[x][x] == 0) return;
//消去A的第x列除去x行以外的各行元素
temp = tt[x][x];
for(y=0;y<19;y++)
{
tt[x][y] = tt[x][y] / temp; //主对角线上元素变成1
tminv[x][y] = tminv[x][y] / temp; //伴随计算
}
for(y=0;y<19;y++) //行 0 -> n
{
if(y!=x) //不是第x行
{
temp = tt[y][x];
for(k=0;k<19;k++) // y行元素 - x行元素* y行x列元素
{
tt[y][k] = tt[y][k] - tt[x][k] * temp;
tminv[y][k] = tminv[y][k] - tminv[x][k] * temp;
}
}
}
}
return;
}