📄 matrix.cpp
字号:
{ t=value(j,k-1);
if (fabs(t)>fabs(d))
{ d=t; i=j;}
}
if (fabs(d)!=0.0)
{ if (i!=k)
{ for (j=k-1; j<rownum; j++)
{
t = value(i,j);
set(i,j,value(k,j));
set(k,j,t);
}
for (j=0; j<rownum; j++)
{
t = value(j,i);
set(j,i,value(j,k));
set(j,k,t);
}
}
for (i=k+1; i<rownum; i++)
{
t = value(i,k-1)/d;
set(i,k-1,0.0);
for (j=k; j<rownum; j++)
set(i,j,value(i,j)-t*value(k,j));
for (j=0; j<rownum; j++)
set(j,k,value(j,k)+t*value(j,i));
}
}
}
}
void matrix::qreigen(matrix & u1, matrix & v1, size_t jt, DOUBLE eps)
// 求一般实矩阵的所有特征根
// a和b均返回rownum行一列的列向量矩阵，返回所有特征根的实部和虚部
{
matrix a(*this);
a.hessenberg(); // 先算出赫申伯格矩阵
u1 = matrix(rownum);
v1 = matrix(rownum);
buffer * uu = getnewbuffer(rownum);
buffer * vv = getnewbuffer(rownum);
buffer &u = *uu;
buffer &v = *vv;
size_t m,it,i,j,k,l;
size_t iir,iic,jjr,jjc,kkr,kkc,llr,llc;
DOUBLE b,c,w,g,xy,p,q,r,x,s,e,f,z,y;
it=0; m=rownum;
while (m!=0)
{ l=m-1;
while ( l>0 && (fabs(a.value(l,l-1))>eps*
(fabs(a.value(l-1,l-1))+fabs(a.value(l,l))))) l--;
iir = m-1; iic = m-1;
jjr = m-1; jjc = m-2;
kkr = m-2; kkc = m-1;
llr = m-2; llc = m-2;
if (l==m-1)
{ u[m-1]=a.value(m-1,m-1); v[m-1]=0.0;
m--; it=0;
}
else if (l==m-2)
{ b=-(a.value(iir,iic)+a.value(llr,llc));
c=a.value(iir,iic)*a.value(llr,llc)-
a.value(jjr,jjc)*a.value(kkr,kkc);
w=b*b-4.0*c;
y=sqrt(fabs(w));
if (w>0.0)
{ xy=1.0;
if (b<0.0) xy=-1.0;
u[m-1]=(-b-xy*y)/2.0;
u[m-2]=c/u[m-1];
v[m-1]=0.0; v[m-2]=0.0;
}
else
{ u[m-1]=-b/2.0; u[m-2]=u[m-1];
v[m-1]=y/2.0; v[m-2]=-v[m-1];
}
m=m-2; it=0;
}
else
{
if (it>=jt) {
delete uu;
delete vv;
throw TMESSAGE("fail to calculate eigenvalue");
}
it++;
for (j=l+2; j<m; j++)
a.set(j,j-2,0.0);
for (j=l+3; j<m; j++)
a.set(j,j-3,0.0);
for (k=l; k+1<m; k++)
{ if (k!=l)
{ p=a.value(k,k-1); q=a.value(k+1,k-1);
r=0.0;
if (k!=m-2) r=a.value(k+2,k-1);
}
else
{
x=a.value(iir,iic)+a.value(llr,llc);
y=a.value(llr,llc)*a.value(iir,iic)-
a.value(kkr,kkc)*a.value(jjr,jjc);
iir = l; iic = l;
jjr = l; jjc = l+1;
kkr = l+1; kkc = l;
llr = l+1; llc = l+1;
p=a.value(iir,iic)*(a.value(iir,iic)-x)
+a.value(jjr,jjc)*a.value(kkr,kkc)+y;
q=a.value(kkr,kkc)*(a.value(iir,iic)+a.value(llr,llc)-x);
r=a.value(kkr,kkc)*a.value(l+2,l+1);
}
if ((fabs(p)+fabs(q)+fabs(r))!=0.0)
{ xy=1.0;
if (p<0.0) xy=-1.0;
s=xy*sqrt(p*p+q*q+r*r);
if (k!=l) a.set(k,k-1,-s);
e=-q/s; f=-r/s; x=-p/s;
y=-x-f*r/(p+s);
g=e*r/(p+s);
z=-x-e*q/(p+s);
for (j=k; j<m; j++)
{
iir = k; iic = j;
jjr = k+1; jjc = j;
p=x*a.value(iir,iic)+e*a.value(jjr,jjc);
q=e*a.value(iir,iic)+y*a.value(jjr,jjc);
r=f*a.value(iir,iic)+g*a.value(jjr,jjc);
if (k!=m-2)
{ kkr = k+2; kkc = j;
p=p+f*a.value(kkr,kkc);
q=q+g*a.value(kkr,kkc);
r=r+z*a.value(kkr,kkc);
a.set(kkr,kkc,r);
}
a.set(jjr,jjc,q);
a.set(iir,iic,p);
}
j=k+3;
if (j>=m-1) j=m-1;
for (i=l; i<=j; i++)
{ iir = i; iic = k;
jjr = i; jjc = k+1;
p=x*a.value(iir,iic)+e*a.value(jjr,jjc);
q=e*a.value(iir,iic)+y*a.value(jjr,jjc);
r=f*a.value(iir,iic)+g*a.value(jjr,jjc);
if (k!=m-2)
{ kkr = i; kkc = k+2;
p=p+f*a.value(kkr,kkc);
q=q+g*a.value(kkr,kkc);
r=r+z*a.value(kkr,kkc);
a.set(kkr,kkc,r);
}
a.set(jjr,jjc,q);
a.set(iir,iic,p);
}
}
}
}
}
for(i=0;i<rownum;i++) {
u1.set(i,u[i]);
v1.set(i,v[i]);
}
delete uu;
delete vv;
}
DOUBLE gassrand(int rr) // 返回一零均值单位方差的正态分布随机数
{
static DOUBLE r=3.0; // 种子
if(rr) r = rr;
int i,m;
DOUBLE s,w,v,t;
s=65536.0; w=2053.0; v=13849.0;
t=0.0;
for (i=1; i<=12; i++)
{ r=r*w+v; m=(int)(r/s);
r-=m*s; t+=r/s;
}
t-=6.0;
return(t);
}
gassvector::gassvector(matrix & r):a(r) //r必须是正定对称阵，为正态随机向量的协方差
{
if(r.rownum != r.colnum) throw TMESSAGE("must be a sqare matrix");//必顶为方阵
if(!r.issym) throw TMESSAGE("must be a symetic matrix");//必须为正定矩阵
matrix evalue;
a = a.eigen(evalue);
for(size_t i=0; i<a.colnum; i++)
{
DOUBLE e = sqrt(evalue(i));
for(size_t j=0; j<a.rownum; j++)
a.set(j,i,a.value(j,i)*e);
}
}
matrix gassvector::operator()(matrix & r) // 返回给定协方差矩阵的正态随机向量
{
a = r;
matrix evalue;
a = a.eigen(evalue);
size_t i;
for(i=0; i<a.colnum; i++) {
DOUBLE e = sqrt(evalue(i));
for(size_t j=0; j<a.rownum; j++)
a.set(j,i,a.value(j,i)*e);
}
matrix rr(a.rownum); // 产生列向量
for(i=0; i<a.rownum; i++)
rr.set(i,gassrand());
return a*rr;
}
matrix gassvector::operator()() // 返回已设定的协方差矩阵的正态随机向量
{
matrix rr(a.rownum);
for(size_t i=0; i<a.rownum; i++)
rr.set(i,gassrand());
return a*rr;//列向量
}
void gassvector::setdata(matrix & r) // 根据协方差矩阵设置增益矩阵
{
if(!r.issym) throw TMESSAGE("r must be symetic!");
a = r;
matrix evalue;
a = a.eigen(evalue);
for(size_t i=0; i<a.colnum; i++) {
if(evalue(i)<0.0) throw TMESSAGE("var matrix not positive!");
DOUBLE e = sqrt(evalue(i));
for(size_t j=0; j<a.rownum; j++)
a.set(j,i,a.value(j,i)*e);
}
}
int matrix::ispositive() // 判定矩阵是否对称非负定阵，如是返回1，否则返回0
{
if(!issym) return 0;
matrix ev;
eigen(ev);
for(size_t i=0; i<rownum; i++)
if(ev(i)<0) return 0;
return 1;
}
matrix matrix::cholesky(matrix& dd) // 用乔里斯基分解法求对称正定阵的线性
// 方程组ax=d返回方程组的解，本身为a，改变为分解阵u,d不变
{
if(!issym) throw TMESSAGE("not symetic!");
if(dd.rownum != colnum) throw TMESSAGE("dd's rownum not right!");
matrix md(dd);
size_t i,j,k,u,v;
if(value(0,0)<=0.0) throw TMESSAGE("not positive");
set(0,0,sqrt(value(0,0))); // a[0]=sqrt(a[0]);
buffer& a = (*buf);
buffer& d = (*(md.buf));
size_t n = rownum;
size_t m = dd.colnum;
for (j=1; j<n; j++) a[j]=a[j]/a[0];
for (i=1; i<n; i++)
{ u=i*n+i;
for (j=1; j<=i; j++)
{ v=(j-1)*n+i;
a[u]=a[u]-a[v]*a[v];
}
if (a[u]<=0.0) throw TMESSAGE("not positive");
a[u]=sqrt(a[u]);
if (i!=(n-1))
{ for (j=i+1; j<n; j++)
{ v=i*n+j;
for (k=1; k<=i; k++)
a[v]=a[v]-a[(k-1)*n+i]*a[(k-1)*n+j];
a[v]=a[v]/a[u];
}
}
}
for (j=0; j<m; j++)
{ d[j]=d[j]/a[0];
for (i=1; i<=n-1; i++)
{ u=i*n+i; v=i*m+j;
for (k=1; k<=i; k++)
d[v]=d[v]-a[(k-1)*n+i]*d[(k-1)*m+j];
d[v]=d[v]/a[u];
}
}
for (j=0; j<=m-1; j++)
{ u=(n-1)*m+j;
d[u]=d[u]/a[n*n-1];
for (k=n-1; k>=1; k--)
{ u=(k-1)*m+j;
for (i=k; i<=n-1; i++)
{ v=(k-1)*n+i;
d[u]=d[u]-a[v]*d[i*m+j];
}
v=(k-1)*n+k-1;
d[u]=d[u]/a[v];
}
}
if(n>1)
for(j=1; j<n; j++)
for(i=0; i<j; i++)
a[i+j*n]=0.0;
return md;
}
DOUBLE lineopt(matrix& aa,matrix& bb, matrix& cc, matrix & xx)
// 线性规划最优点寻找程序，aa为mXn不等式约束条件左端系数矩阵，bb为不等式约束
// 条件的右端项，为m维向量，cc为目标函数系数，n维向量，xx返回极小点，为n维向量
{
if(aa.rownum != bb.rownum || aa.colnum != cc.rownum ||
aa.colnum != xx.rownum) throw TMESSAGE("dimenstion not right!");
size_t n=aa.colnum, m=aa.rownum;
size_t i,mn,k,j;
matrix a(m,n+m);
for(i=0;i<m;i++) {
for(j=0;j<n;j++)
a.set(i,j,aa(i,j));
for(j=n;j<n+m;j++)
if(j-n == i) a.set(i,j,1.0);
else a.set(i,j,0.0);
}
matrix c(m+n);
c = 0.0;
for(i=0;i<m;i++)
c.set(i,cc(i));
lbuffer* jjs = getnewlbuffer(m);
lbuffer& js = (*jjs);
DOUBLE s,z,dd,y; //,*p,*d;
for (i=0; i<m; i++) js[i]=n+i;
matrix p(m,m);
matrix d;
mn=m+n; s=0.0;
matrix x(mn);
while (1)
{ for (i=0; i<m; i++)
for (j=0; j<m; j++)
p.set(i,j,a(i,(size_t)js[j]));
p.inv();
d = p*a;
x = 0.0;
for (i=0; i<m; i++)
{ s=0.0;
for (j=0; j<=m-1; j++)
s+=p(i,j)*bb(j);
x.set((size_t)js[i],s);
}
k=mn; dd=1.0e-35;
for (j=0; j<mn; j++)
{ z=0.0;
for (i=0; i<m; i++)
z+=c((size_t)js[i])*d(i,j);
z-=c(j);
if (z>dd) { dd=z; k=j;}
}
if (k==mn)
{ s=0.0;
for (j=0; j<n; j++) {
xx.set(j,x(j));
s+=c(j)*x(j);
}
delete jjs;
return s;
}
j=m;
dd=1.0e+20;
for (i=0; i<=m-1; i++)
if (d(i,k)>=1.0e-20) // [i*mn+k]>=1.0e-20)
{ y=x(size_t(js[i]))/d(i,k);
if (y<dd) { dd=y; j=i;}
}
if (j==m) { delete jjs;
throw TMESSAGE("lineopt failed!");
}
js[j]=k;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -