📄 main.cpp
字号:
while(stem[tt]!=icel&&tt<snu)
{
tt++;
}
if(tt==snu)/////查看是否是其他吸引域里面
{
long iflgf;
iflgf=flag[icel];
//////////已经确定的cell也不能加进去
long ddi;
ddi=fder[icel];
if(iflgf==i&&ddi!=i) ////
{
stem[snu]=icel;
snu++;
afn++;
}
}
}
}
}
///////////////// 此小批已经获得,放在stem[] , 有snu个 /////////////////////////////
///////// snu 个stem排序
long *sstem;
sstem=new long [snu];
for(long uu=0;uu<snu;uu++)
{
sstem[uu]=stem[uu];
}
/////已经胞排好了,放在sstem[snu]中
////开始计算
long threm[5000000];
long thrn=0;
long nhand=0;
while(nhand<snu)
{
thrn=0;
long cre;
long tt=snu-1;
if((snu-nhand)<200)
{
tt=0;
while(sstem[tt]==-1)///////////////从最后一个开始算
{
tt++;
}
}
else
{
tt=snu-1;
while(sstem[tt]==-1)///////////////从最后一个开始算
{
tt--;
}
}
cre=sstem[tt];
threm[thrn]=cre;
thrn++;
long an=0;
long bn=-1;
long psta=0;
while(an>bn)
{
bn=an;
for(long u=psta;u<thrn;u++)
{
long crm=threm[u];
psta++;
long Iz,loc,loc1;
gcm=fgcm[crm];
Iz=gcm.Iz;
loc=gcm.CpLoc;
for(long t=0;t<Iz;t++)
{
loc1=loc+t;
cp=fcpl[loc1];
long cie=cp.Cz;
long det;
det=fder[cie];
long flt;
flt=flag[cie];
if(det!=i&&flt==i)
{
//////保存了吗
long ttt=0;
while(threm[ttt]!=cie&&ttt<thrn)
{
ttt++;
}
if(ttt==thrn)
{
an++;
threm[thrn]=cie;
thrn++;
}
}
}
}
}
//////////////不稳定组, //////////
//////计算,已经放到threm[thrn]中
////////pj////////////////
handnum=handnum+thrn;
cout<<"total : "<<numb<<" handled : "<<handnum<<" equation : "<<thrn<<endl;
float *xr,*xr2;
xr=new float [thrn];/////////迭代用
xr2=new float [thrn];
float err=1.0;
for(k=0;k<thrn;k++)
{
xr[k]=0.0;
}
xr[0]=0.65; ////////////////初值
if(thrn<1000)
{
mat.Realloc(thrn,thrn);
mat=0.0;
}
xlink **head;
xlink *p;
head=new xlink *[thrn];
if(thrn>=1000)
{
for(int j=0;j<thrn;j++)
{
head[j]=new xlink [1];
if(j>=1)
head[j]->num=0;
}
}
cst.Realloc(thrn,1);
cst=0.0;
float *cs;
cs=new float [thrn];
if(thrn<1000)
{
for(k=0;k<thrn;k++)
{
if(thrn>1000&&thrn%500==0)
cout<<"get xishu "<<k<<" "<<thrn<<endl;
cs[k]=0.0;
long ce=threm[k];
long Iz;
gcm=fgcm[ce];
Iz=gcm.Iz;
long loc2=gcm.CpLoc;
mat.m_adValue[k][k]=1.0;///////// 别忘记了
for(long t=0;t<Iz;t++)
{
long loc=loc2+t;
cp=fcpl[loc];
long ceim=cp.Cz;////////
long tt=0;
while(ceim!=threm[tt]&&tt<thrn) /////// ????????tt<tnb
{
tt++;
}
if(tt==thrn) /////////象为已知胞, 也可能为未知胞
{ ////看 到
long ttfg;
ttfg=flag[ceim];
long dfdf;
dfdf=fder[ceim];
if(ttfg==i&&dfdf==i)
{
float pj;
pj=fres[ceim][i-1];
cs[k]=cs[k]+cp.Pz*pj; ////常数项
}
}
else
{
if(tt==k)////////与ce相同
{
////ce 到 ce 的 概率
mat.m_adValue[k][k]=mat.m_adValue[k][k]-cp.Pz;
}
else
{
mat.m_adValue[k][tt]=-cp.Pz;
}
}
}
}
}
else
{
for(k=0;k<thrn;k++)
{
if(thrn>1000&&thrn%500==0)
cout<<"get xishu "<<k<<" "<<thrn<<endl;
cs[k]=0.0;
long ce=threm[k];
long Iz;
gcm=fgcm[ce];
Iz=gcm.Iz;
long loc2=gcm.CpLoc;
p=head[k];
p->num=k;
p->x=1.0;
//mat.m_adValue[k][k]=1.0;///////// 别忘记了
for(long t=0;t<Iz;t++)
{
long loc=loc2+t;
cp=fcpl[loc];
long ceim=cp.Cz;////////
long tt=0;
while(ceim!=threm[tt]&&tt<thrn) /////// ????????tt<tnb
{
tt++;
}
if(tt==thrn) /////////象为已知胞, 也可能为未知胞
{ ////看 到
long ttfg;
ttfg=flag[ceim];
long dfdf;
dfdf=fder[ceim];
if(ttfg==i&&dfdf==i)
{
float pj;
pj=fres[ceim][i-1];
cs[k]=cs[k]+cp.Pz*pj; ////常数项
}
}
else
{
if(tt==k)////////与ce相同
{
////ce 到 ce 的 概率
//mat.m_adValue[k][k]=mat.m_adValue[k][k]-cp.Pz;
p=head[k];
p->x=p->x-cp.Pz;
}
else
{
//mat.m_adValue[k][tt]=-cp.Pz;
p=head[k];
while(p->next!=NULL)
{
p=p->next;
}
p->next=new xlink();
p=p->next;
p->x=-cp.Pz;
p->num=tt;
}
}
}
}
}
/////////上面获得系数
for(long yu=0;yu<thrn;yu++)
{
cst.m_adValue[yu][0]=cs[yu];
}
if(thrn<1000)
{
mat.MatInv();
cst=mat*cst;
///////////////////////
}
else
{
err=1.0;
while(err>0.01)
{
float tot=0.0;
for(long tr=0;tr<thrn;tr++)
{
if(tr%500==0)
cout<<tr<<" "<<thrn<<endl;
tot=0.0;
p=head[tr];
p=p->next;/////第二个胞了，非tr
while(p!=NULL)
{
tot=tot+(p->x)*xr[p->num];
p=p->next;
}
p=head[tr];
float njt=p->x;
xr2[tr]=(cs[tr]-tot)/njt;
}
err=0.0;
for(long jr=0;jr<thrn;jr++)
{
err=err+sqrt((xr2[jr]-xr[jr])*(xr2[jr]-xr[jr]));
}
for(jr=0;jr<thrn;jr++)
{
xr[jr]=xr2[jr];
}
cout<<err<<" "<<err<<endl;
}
////////////结果放在cst
for(long te=0;te<thrn;te++)
{
cst.m_adValue[te][0]=xr[te];
}
for(long mt=0;mt<thrn;mt++)
{
xlink *q1,*q2;
q1=head[mt]->next;
while(q1!=NULL)
{
q2=q1->next;
delete q1;
q1=q2;
}
}
}
////存入数据
for(long tr=0;tr<thrn;tr++)
{
long tre=threm[tr];
float pjr=cst.m_adValue[tr][0];
fres[tre][i-1]=pjr;
}
/////////////vj////////////////
delete [] cs;
delete [] xr;
delete [] xr2;
////////结果,改变sstem -1 以及ioder
for(long u=0;u<thrn;u++)
{
long cme=threm[u];
long det=i;
fder[cme]=det;
long tr=0;
while(sstem[tr]!=cme&&tr<snu)
{
tr++;
}
sstem[tr]=-1;
}
nhand=nhand+thrn;
//for(long loop=0;loop<thrn;loop++)//////////////////清理
delete [] head;
//delete [] head;
}
//delete [] stt;
//delete [] stn;
delete [] sstem;
///////////////////////////////////////
for(long u=0;u<snu;u++)
{
long cew;
cew=stem[u];
long yu=i;
fder[cew]=yu;
/////////////////dtem 中 -1
long tt=0;
while(dtem[tt]!=cew&&tt<nu)
{
tt++;
}
dtem[tt]=-1;/////////////////////
}
nt=nt+snu;
}
}
cet=ftemp[sta];
sta++;
///////////////此批计算完//////////////////////
}
}
////////////////////////
fstream ioresult; ////// 保存结果 gr,Vj,pj
ioresult.open("result.dat",ios::in|ios::out|ios::binary);
///long res=3*sizeof(float);/////////////////////////////////////////
for(long mt=0;mt<max;mt++)
{
for(long t=0;t<3;t++)
{
ioresult.write((char *)(&fres[mt][t]),sizeof(float));
}
}
delete [] fzpre;
delete [] flag;
delete [] fder;
delete [] fcpl;
delete [] fgcm;
delete [] ftemp;
}
////////////排序算法
void sort(long a[],long b[],long n)
{
long i,j,small;
long temp1,temp2;
for(i=0;i<n-1;i++)
{
small=i;
for(j=i+1;j<n;j++)
{
if(a[j]<a[small])
small=j;
}
if(small!=i)
{
temp1=a[i];
temp2=b[i];
a[i]=a[small];
a[small]=temp1;
b[i]=b[small];
b[small]=temp2;
}
}
}
long pai(long m)
{
long n=g_attr;//
long nt=1;
long mt=1;
if(m>n)
{
return -1;
}
else
{
for(long i=1;i<=m;i++)
{
nt=nt*n;
n=n-1;
}
if(m==0)
{
mt=1;
}
else
{
for(i=1;i<=m;i++)
{
mt=mt*i;
}
}
long re;
re=(long) (nt/mt);
return re;
}
}
void fenzu()
{
///// 根据 iodm 中的信息,改变ioresult 中的gr值
fstream iodm; ////// 保存 cell 的, domicle 的数量,做分组用
iodm.open("todomiledata.dat",ios::in|ios::out|ios::binary);
dom dm;
fstream ioresult; ////// 保存结果 gr,Vj,pj
ioresult.open("result.dat",ios::in|ios::out|ios::binary);
long res=sizeof(long)+2*g_attr*sizeof(float);//////////////////////////
long pai(long m);
for(long i=0;i<max;i++)
{
if(i%5000==0)
cout<<i<<endl;
long dn;
long total=0;
iodm.seekg(i*sizeof(dom));
iodm.read((char *)(&dm),sizeof(dom));
dn=dm.n;
total=dm.total;
long gpai=0;
if(dn==0)
{
ioresult.seekp(i*res);
ioresult.write((char *)(&total),sizeof(long));
}
}
///////////////////////////////////
}
void outputdata()
{
fstream iopersistent;
iopersistent.open("persistent.dat",ios::in|ios::out|ios::binary);
//g_attr=2;
long tr=0;
ofstream df("persistent.txt");
iopersistent.seekg(0,ios::beg);
long cur,pend;
cur=iopersistent.tellg();
iopersistent.seekg(0,ios::end);
pend=iopersistent.tellg();
iopersistent.seekg(cur);
while(pend!=iopersistent.tellg())
{
long cee,grg;
iopersistent.read((char *)(&cee),sizeof(long));
iopersistent.read((char *)(&grg),sizeof(long));
float x1,x2,x3,x4;
long z1,z2,z3,z4,m;
m=cee;
z4=(long)((m-1)/(step1*step2*step3))+1;
m=m-(z4-1)*step1*step2*step3;
z3=(long)((m-1)/(step1*step2))+1;
m=m-(z3-1)*step1*step2;
z2=(long) ((m-1)/step1)+1;
z1=m-(z2-1)*step1;
x1=z1*h1-h1+x01;
x2=z2*h2-h2+x02;
x3=z3*h3-h3+x03;
x4=z4*h4-h4+x04;
df<<x1<<" "<<x2<<" "<<x3<<" "<<x4<<endl;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -