📄 matlib.cpp
字号:
alfa=-c[k]*c[k+1]/delta;
beta=d[k]*c[k+1]/delta;
gama=c[k+1]*beta;
c[k]=alfa;
I[k]=2;
c[k+1]=beta;
I[k+1]=1;
d[k+2]-=gama;
k++;
}
}
}
else
{
if(k==(mr-2))
{
dv[k]=c[k];
c[k]=0;
I[k]=2;
I[k+1]=0;
k++;
}
else
{
I[k]=0;
}
}
}
vector v;
v=b*q;
b[0]=v[0];
b[1]=v[1];
if(I[0]==1)
b[1]-=c[0]*b[0];
for(k=2;k<mr;k++)
{
b[k]=v[k];
if(I[k-1]==1)
b[k]-=c[k-1]*b[k-1];
if(I[k-2]==2)
b[k]-=c[k-2]*b[k-2];
}
for(k=0;k<mr;k++)
{
double data;
if(I[k]>1)
{
fmindir(b[k],b[k+1],d[k],d[k+1],dv[k]);
k++;
}
else
{
data=fabs(d[k]);
if(data>0.0000001)
b[k]/=-data;
else
b[k]=-b[k];
}
}
v[mr-1]=b[mr-1];
for(k=mr-2;k>=0;k--)
{
int n=k+I[k];
delta=c[k]*v[k+I[k]];
v[k]=b[k]-c[k]*v[k+I[k]];
}
b=q*v;
delete []I;
}
catch(char* str)
{
cout<<str<<endl;
}
catch(...)
{
cout<<"Solve MNewton directon other exception"<<endl;
}
return 1;
}
int SolveMNewtonDir1(matrix& a,vector& b)
{
int mr=mrow(a);
int* Li,*Ri,*I;
I=new int[mr];
Li=new int[mr];
Ri=new int[mr];
try
{
matrix L;
L=eye(mr,mr);
double eps=DBL_EPSILON;
double delta;
int S=0;
for(int k=0;k<mr;k++)
{
Li[k]=-1;
Ri[k]=-1;
}
for(k=0;k<mr;k++)
{
double dmax=fabs(a(k,k));
double data;
int j=k;
for(int i=k+1;i<mr;i++)
{
data=fabs(a(k,i));
if(dmax<data)
{
dmax=data;
j=i;
}
}
if(dmax<=eps)
continue;
if( k<(mr-2) )
{
if(j==k)
{
I[k]=1;
dmax=a(k,k);
for(i=k+1;i<mr;i++)
L(i,k)=a(i,k)/dmax;
for(i=k+1;i<mr;i++)
{
for(int t=i;t<mr;t++)
a(t,i)-=L(i,k)*a(k,t);
}
for(i=k+1;i<mr;i++)
{
for(int t=i+1;t<mr;t++)
a(i,t)=a(t,i);
}
}
else
{
dmax=fabs(a(k,k));
data=fabs(a(k,k)*a(j,j)-a(j,k)*a(j,k));
if(data<=dmax)
{
I[k]=1;
dmax=a(k,k);
for(i=k+1;i<mr;i++)
L(i,k)=a(i,k)/dmax;
for(i=k+1;i<mr;i++)
{
for(int t=i;t<mr;t++)
a(t,i)-=L(i,k)*a(k,t);
}
for(i=k+1;i<mr;i++)
{
for(int t=i+1;t<mr;t++)
a(i,t)=a(t,i);
}
}
else
{
if(j>(k+1))
{
Li[S]=k+1;
Ri[S]=j;
S++;
for(int t=0;t<k;t++)
{
data=L(k+1,t);
L(k+1,t)=L(j,t);
L(j,t)=data;
}
for(t=k;t<mr;t++)
{
data=a(k+1,t);
a(k+1,t)=a(j,t);
a(j,t)=data;
}
for(t=k;t<mr;t++)
{
data=a(t,k+1);
a(t,k+1)=a(t,j);
a(t,j)=data;
}
data=b[k+1];
b[k+1]=b[j];
b[j]=data;
}
I[k]=2;
I[k+1]=1;
dmax=a(k,k)*a(k+1,k+1)-a(k+1,k)*a(k+1,k);
for(i=k+2;i<mr;i++)
{
L(i,k)=(a(i,k)*a(k+1,k+1)-a(i,k+1)*a(k+1,k))/dmax;
L(i,k+1)=(a(i,k+1)*a(k,k)-a(i,k)*a(k+1,k))/dmax;
}
for(i=k+2;i<mr;i++)
{
for(j=i;j<mr;j++)
a(j,i)-=(L(i,k)*a(k,j)+L(i,k+1)*a(k+1,j));
}
for(i=k+2;i<mr;i++)
{
for(j=i+1;j<mr;j++)
a(i,j)=a(j,i);
}
k++;
}
}
}
else
{
if(k==(mr-2))
{
if(j==k)
{
I[k]=1;
dmax=a(k,k);
j=k+1;
L(j,k)=a(j,k)/dmax;
a(j,j)-=L(j,k)*a(k,j);
a(k,j)=a(j,k);
}
else
{
I[k]=2;
I[k+1]=0;
k++;
}
}
else
{
I[k]=0;
}
}
}
vector v(mr);
for(k=0;k<mr;k++)
{
delta=b[k];
for(int j=0;j<k;j++)
delta-=v[j]*L(k,j);
v[k]=delta;
}
for(k=0;k<mr;k++)
{
double data;
if(I[k]>1)
{
fmindir(v[k],v[k+1],a(k,k),a(k+1,k+1),a(k+1,k));
k++;
}
else
{
data=fabs(a(k,k));
if(data>0.0000001)
v[k]/=-data;
else
v[k]=-b[k];
}
}
for(k=mr-1;k>=0;k--)
{
delta=v[k];
for(int j=k+1;j<mr;j++)
delta-=b[j]*L(j,k);
b[k]=delta;
}
if(S>0)
for(k=S-1;k>=0;k--)
{
int j=Li[k];
int i=Ri[k];
delta=b[j];
b[j]=b[i];
b[i]=delta;
}
}
catch(char* str)
{
cout<<str<<endl;
}
catch(...)
{
cout<<"Solve MNewton directon other exception"<<endl;
}
delete []I;
delete []Li;
delete []Ri;
return 1;
}
static void ppp(matrix& a,vector& e,vector& s,matrix& v,int m,int n)
{
int i,j;
double d;
try
{
if (m>=n) i=n;
else i=m;
for (j=1; j<=i-1; j++){
a[j-1][j-1]=s[j-1];
a[j-1][j]=e[j-1];
}
a[i-1][i-1]=s[i-1];
if (m<n) a[i-1][i]=e[i-1];
for (i=1; i<=n-1; i++)
for (j=i+1; j<=n; j++){
d=v[i-1][j-1]; v[i-1][j-1]=v[j-1][i-1];
v[j-1][i-1]=d;
}
}
catch(char* str)
{
cout<<str<<endl;
}
catch(...)
{
cout<<"Solve MNewton directon other exception"<<endl;
}
}
static void sss(double fg[],double cs[])
{
double r,d;
try
{
if ((fabs(fg[0])+fabs(fg[1]))<0.000000001){
cs[0]=1.0; cs[1]=0.0; d=0.0;
}
else
{
d=sqrt(fg[0]*fg[0]+fg[1]*fg[1]);
if (fabs(fg[0])>fabs(fg[1])){
d=fabs(d);
if (fg[0]<0.0) d=-d;
}
if (fabs(fg[1])>=fabs(fg[0])){
d=fabs(d);
if (fg[1]<0.0) d=-d;
}
cs[0]=fg[0]/d; cs[1]=fg[1]/d;
}
r=1.0;
if (fabs(fg[0])>fabs(fg[1])) r=cs[1];
else
if (fabs(cs[0])>0.0000000001) r=1.0/cs[0];
fg[0]=d; fg[1]=r;
}
catch(char* str)
{
cout<<str<<endl;
}
catch(...)
{
cout<<"Solve MNewton directon other exception"<<endl;
}
}
void svd(matrix& a, matrix& u, matrix& v,double eps)
{
int m=a.mrow(),n=a.mcol();
int ka=MAX(m,n)+1;
int i,j,k,l,it,ll,kk,mm,nn,m1,ks;
double d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];
vector s(ka),e(ka),w(ka);
it=60; k=n;
if (m-1<n) k=m-1;
l=m;
if (n-2<m) l=n-2;
if (l<0) l=0;
ll=k;
if (l>k) ll=l;
if (ll>=1){
for (kk=1; kk<=ll; kk++){
if (kk<=k){
d=0.0;
for (i=kk; i<=m; i++)
d=d+a[i-1][kk-1]*a[i-1][kk-1];
s[kk-1]=sqrt(d);
if (s[kk-1]>0.000000001)
{
if (fabs(a[kk-1][kk-1])>0.000000001){
s[kk-1]=fabs(s[kk-1]);
if (a[kk-1][kk-1]<0.0)
s[kk-1]=-s[kk-1];
}
for (i=kk; i<=m; i++)
a[i-1][kk-1]/=s[kk-1];
a[kk-1][kk-1]+=1.0;
}
s[kk-1]=-s[kk-1];
}
if (n>=kk+1){
for (j=kk+1; j<=n; j++){
if ((kk<=k)&&(fabs(s[kk-1])>0.000000001)){
d=0.0;
for (i=kk; i<=m; i++)
d+=a[i-1][kk-1]*a[i-1][j-1];
d=-d/a[kk-1][kk-1];
for (i=kk; i<=m; i++)
a[i-1][j-1]+=d*a[i-1][kk-1];
}
e[j-1]=a[kk-1][j-1];
}
}
if (kk<=k)
for (i=kk; i<=m; i++)
u[i-1][kk-1]=a[i-1][kk-1];
if (kk<=l){
d=0.0;
for (i=kk+1; i<=n; i++)
d=d+e[i-1]*e[i-1];
e[kk-1]=sqrt(d);
if (e[kk-1]>0.000000001){
if (fabs(e[kk])>0.000000001){
e[kk-1]=fabs(e[kk-1]);
if (e[kk]<0.0)
e[kk-1]=-e[kk-1];
}
for (i=kk+1; i<=n; i++)
e[i-1]/=e[kk-1];
e[kk]=1.0+e[kk];
}
e[kk-1]=-e[kk-1];
if ((kk+1<=m)&&(e[kk-1]!=0.0)){
for (i=kk+1; i<=m; i++)
w[i-1]=0.0;
for (j=kk+1; j<=n; j++)
for (i=kk+1; i<=m; i++)
w[i-1]=w[i-1]+e[j-1]*a[i-1][j-1];
for (j=kk+1; j<=n; j++)
for (i=kk+1; i<=m; i++)
a[i-1][j-1]-=w[i-1]*e[j-1]/e[kk];
}
for (i=kk+1; i<=n; i++)
v[i-1][kk-1]=e[i-1];
}
}
}
mm=n;
if (m+1<n) mm=m+1;
if (k<n) s[k]=a[k][k];
if (m<mm) s[mm-1]=0.0;
if (l+1<mm) e[l]=a[l][mm-1];
e[mm-1]=0.0;
nn=m;
if (m>n) nn=n;
if (nn>=k+1)
for (j=k+1; j<=nn; j++){
for (i=1; i<=m; i++)
u[i-1][j-1]=0.0;
u[j-1][j-1]=1.0;
}
if (k>=1){
for (ll=1; ll<=k; ll++){
kk=k-ll+1;
if (s[kk-1]!=0.0){
if (nn>=kk+1)
for (j=kk+1; j<=nn; j++){
for (d=0.0,i=kk; i<=m; i++)
d+=u[i-1][kk-1]*u[i-1][j-1]/u[kk-1][kk-1];
d=-d;
for (i=kk; i<=m; i++)
u[i-1][j-1]+=d*u[i-1][kk-1];
}
for (i=kk; i<=m; i++)
u[i-1][kk-1]=-u[i-1][kk-1];
u[kk-1][kk-1]+=1.0;
if (kk-1>=1)
for (i=1; i<=kk-1; i++)
u[i-1][kk-1]=0.0;
}
else{
for (i=1; i<=m; i++)
u[i-1][kk-1]=0.0;
u[kk-1][kk-1]=1.0;
}
}
}
for (ll=1; ll<=n; ll++){
kk=n-ll+1;
if ((kk<=l)&&(fabs(e[kk-1])>0.000000001)){
for (j=kk+1; j<=n; j++){
for (d=0.0,i=kk+1; i<=n; i++)
d+=v[i-1][kk-1]*v[i-1][j-1]/v[kk][kk-1];d=-d;
for (i=kk+1; i<=n; i++)
v[i-1][j-1]=v[i-1][j-1]+d*v[i-1][kk-1];
}
}
for (i=1; i<=n; i++)
v[i-1][kk-1]=0.0;
v[kk-1][kk-1]=1.0;
}
for (i=1; i<=m; i++)
for (j=1; j<=n; j++)
a[i-1][j-1]=0.0;
m1=mm; it=60;
while (1){
if (mm==0){
ppp(a,e,s,v,m,n);
return;
}
if (it==0){
ppp(a,e,s,v,m,n);
// AfxMessageBox("叠代 60 步,尚未求得某个奇异值");
return;
}
kk=mm-1;
while ((kk!=0)&&(fabs(e[kk-1])!=0.0)){
d=fabs(s[kk-1])+fabs(s[kk]);
dd=fabs(e[kk-1]);
if (dd>eps*d) kk=kk-1;
else e[kk-1]=0.0;
}
if (kk==mm-1){
kk++;
if (s[kk-1]<0.0){
s[kk-1]=-s[kk-1];
for (i=1; i<=n; i++)
v[i-1][kk-1]=-v[i-1][kk-1];
}
while ((kk!=m1)&&(s[kk-1]<s[kk])){
d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d;
if (kk<n)
for (i=1; i<=n; i++){
d=v[i-1][kk-1];
v[i-1][kk-1]=v[i-1][kk];
v[i-1][kk-1]=d;
}
if (kk<m)
for (i=1; i<=m; i++){
d=u[i-1][kk-1];
u[i-1][kk-1]=u[i-1][kk];
u[i-1][kk]=d;
}
kk++;
}
it=60;
mm--;
}
else{
ks=mm;
while ((ks>kk)&&(fabs(s[ks-1])!=0.0)){
d=0.0;
if (ks!=mm) d=d+fabs(e[ks-1]);
if (ks!=kk+1) d=d+fabs(e[ks-2]);
dd=fabs(s[ks-1]);
if (dd>eps*d) ks=ks-1;
else s[ks-1]=0.0;
}
if (ks==kk){
kk=kk+1;
d=fabs(s[mm-1]);
t=fabs(s[mm-2]);
if (t>d) d=t;
t=fabs(e[mm-2]);
if (t>d) d=t;
t=fabs(s[kk-1]);
if (t>d) d=t;
t=fabs(e[kk-1]);
if (t>d) d=t;
if(d>0.000000001)
{
sm=s[mm-1]/d; sm1=s[mm-2]/d;
em1=e[mm-2]/d;
sk=s[kk-1]/d; ek=e[kk-1]/d;
}
b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0;
c=sm*em1; c=c*c; shh=0.0;
if ((fabs(b)>0.000000001)||(c>0.000000001)){
shh=sqrt(b*b+c);
if (b<0.0) shh=-shh;
shh=c/(b+shh);
}
fg[0]=(sk+sm)*(sk-sm)-shh;
fg[1]=sk*ek;
for (i=kk; i<mm; i++){
sss(fg,cs);
if (i!=kk) e[i-2]=fg[0];
fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];
e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];
fg[1]=cs[1]*s[i];
s[i]=cs[0]*s[i];
if ((cs[0]!=1.0)||(cs[1]!=0.0))
for (j=1; j<=n; j++){
d=cs[0]*v[j-1][i-1]+cs[1]*v[j-1][i];
v[j-1][i]=-cs[1]*v[j-1][i-1]+cs[0]*v[j-1][i];
v[j-1][i-1]=d;
}
sss(fg,cs);
s[i-1]=fg[0];
fg[0]=cs[0]*e[i-1]+cs[1]*s[i];
s[i]=-cs[1]*e[i-1]+cs[0]*s[i];
fg[1]=cs[1]*e[i];
e[i]=cs[0]*e[i];
if (i<m)
if ((cs[0]!=1.0)||(cs[1]!=0.0))
for (j=1; j<=m; j++){
d=cs[0]*u[j-1][i-1]+cs[1]*u[j-1][i];
u[j-1][i]=-cs[1]*u[j-1][i-1]+cs[0]*u[j-1][i];
u[j-1][i-1]=d;
}
}
e[mm-2]=fg[0];
it--;
}
else{
if (ks==mm){
kk++;
fg[1]=e[mm-2]; e[mm-2]=0.0;
for (ll=kk; ll<mm; ll++){
i=mm+kk-ll-1;
fg[0]=s[i-1];
sss(fg,cs);
s[i-1]=fg[0];
if (i!=kk){
fg[1]=-cs[1]*e[i-2];
e[i-2]=cs[0]*e[i-2];
}
if ((cs[0]!=1.0)||(cs[1]!=0.0))
for (j=1; j<=n; j++){
d=cs[0]*v[j-1][i-1]+cs[1]*v[j-1][mm-1];
v[j-1][mm-1]=-cs[1]*v[j-1][i-1]+cs[0]*v[j-1][mm-1];
v[j-1][i-1]=d;
}
}
}
else{
kk=ks+1;
fg[1]=e[kk-2];
e[kk-2]=0.0;
for (i=kk; i<=mm; i++){
fg[0]=s[i-1];
sss(fg,cs);
s[i-1]=fg[0];
fg[1]=-cs[1]*e[i-1];
e[i-1]=cs[0]*e[i-1];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -