📄 matrix.h
字号:
for(int i = 0; i < M; i++)
for(int j = 0; j < N; j++)
A[(i+1)*(N+1)+(j+1)] = matA(i, j);
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
COMPLEX_DOUBLE Q,R;
double B[100],C[100],T[100];
// double ETA=1.5E-8,TOL=1.0E-31;
double TOL=DBL_MIN, ETA=DBL_EPSILON;
int NP,N1,K,K1,I,J,L,KK,LL,L1;
double Z,W,CS,SN;
double F,X,Y,G,H;
double EPS;
NP=N+P;
N1=N+1;
//***************************************************************************
//****** HOUSEHOLDER REDUCTION *****************
//***************************************************************************
C[1]=0.0E0;
K=1;
LABEL10:
K1=K+1;
//***************************************************************************
//*************** ELIMINATION OF A[I,K] I=K+1,...,N *********
//***************************************************************************
Z=0.0E0;
for(I=K;I<=M;I++)
Z+=norm(A[I*(N+1)+K]);
B[K]=0.0E0;
if(Z<=TOL)
goto LABEL70;
Z=sqrt(Z);
B[K]=Z;
W=abs(A[K*(N+1)+K]);
Q=COMPLEX_DOUBLE (1.0E0,0.0E0);
if(W!=0.0E0)
Q=A[K*(N+1)+K]/W;
A[K*(N+1)+K]=Q*(Z+W);
if(K==NP)
goto LABEL70;
for(J=K1;J<=NP;J++)
{
Q=COMPLEX_DOUBLE (0.0E0,0.0E0);
for(I=K;I<=M;I++)
Q+=conjg(A[I*(N+1)+K])*A[I*(N+1)+J];
Q=Q/(Z*(Z+W));
for(I=K;I<=M;I++)
A[I*(N+1)+J]=A[I*(N+1)+J]-Q*A[I*(N+1)+K];
}
//***************************************************************************
//**************** PHASE TRANSFORMATION *************
//***************************************************************************
Q=-conjg(A[K*(N+1)+K])/abs(A[K*(N+1)+K]);
for(J=K1;J<=NP;J++)
A[K*(N+1)+J]=Q*A[K*(N+1)+J];
//***************************************************************************
//******************* ELIMINATION OF A[K,J] *************
//***************************************************************************
LABEL70:
if(K==N)
goto LABEL140;
Z=0.0E0;
for(J=K1;J<=N;J++)
Z+=norm(A[K*(N+1)+J]);
C[K1]=0.0E0;
if(Z<=TOL)
goto LABEL130;
Z=sqrt(Z);
C[K1]=Z;
W=abs(A[K*(N+1)+K1]);
Q=COMPLEX_DOUBLE (1.0E0,0.0E0);
if(W!=0.0E0)
Q=A[K*(N+1)+K1]/W;
A[K*(N+1)+K1]=Q*(Z+W);
for(I=K1;I<=M;I++)
{
Q=COMPLEX_DOUBLE (0.0E0,0.0E0);
for(J=K1;J<=N;J++)
Q+=conjg(A[K*(N+1)+J])*A[I*(N+1)+J];
Q=Q/(Z*(Z+W));
for(J=K1;J<=N;J++)
A[I*(N+1)+J]=A[I*(N+1)+J]-Q*A[K*(N+1)+J];
}
//***************************************************************************
//********** PHASE TRANSFORMATION *********
//***************************************************************************
Q=-conjg(A[K*(N+1)+K1])/abs(A[K*(N+1)+K1]);
for(I=K1;I<=M;I++)
A[I*(N+1)+K1]=A[I*(N+1)+K1]*Q;
LABEL130:
K=K1;
goto LABEL10;
//***************************************************************************
//********** TOLERANCE FOR NEGLIGIBLE ELEMENTS *********
//***************************************************************************
LABEL140:
EPS=0.0E0;
for(K=1;K<=N;K++)
{
S[K]=B[K];
T[K]=C[K];
EPS=__max(EPS,S[K]+T[K]);
}
EPS=EPS*ETA;
// EPS=ETA;
//***************************************************************************
//*********** INITIALIZATION OF U AND V *********
//***************************************************************************
if(NU==0)
goto LABEL180;
for(J=1;J<=NU;J++)
{
for(I=1;I<=M;I++)
U[I*(N+1)+J]=COMPLEX_DOUBLE (0.0E0,0.0E0);
U[J*(N+1)+J]=COMPLEX_DOUBLE (1.0E0,0.0E0);
}
LABEL180:
if(NV==0)
goto LABEL210;
for(J=1;J<=NV;J++)
{
for(I=1;I<=N;I++)
V[I*(N+1)+J]=COMPLEX_DOUBLE (0.0E0,0.0E0);
V[J*(N+1)+J]=COMPLEX_DOUBLE (1.0E0,0.0E0);
}
//***************************************************************************
//************** QR DIAGONALIZATION *************
//***************************************************************************
LABEL210:
for(KK=1;KK<=N;KK++)
{
K=N1-KK;
//***************************************************************************
//************* TEST FOR SPLIT *********
//***************************************************************************
LABEL220:
for(LL=1;LL<=K;LL++)
{
L=K+1-LL;
if(fabs(T[L])<=EPS)
{
goto LABEL290;
}
if(fabs(S[L-1])<=EPS)
{
goto LABEL240;
}
}
//***************************************************************************
//*************** CANCELLATION OF E[L] *************
//***************************************************************************
LABEL240:
CS=0.0E0;
SN=1.0E0;
L1=L-1;
for(I=L;I<=K;I++)
{
F=SN*T[I];
T[I]=CS*T[I];
if(fabs(F)<=EPS)
{
goto LABEL290;
}
H=S[I];
W=sqrt(F*F+H*H);
S[I]=W;
CS=H/W;
SN=-F/W;
if(NU==0)
goto LABEL260;
for(J=1;J<=N;J++)
{
X=real(U[J*(N+1)+L1]);
Y=real(U[J*(N+1)+I]);
U[J*(N+1)+L1]=COMPLEX_DOUBLE (X*CS+Y*SN,0.0E0);
U[J*(N+1)+I]=COMPLEX_DOUBLE (Y*CS-X*SN,0.0E0);
}
LABEL260:
if(NP==N)
// goto LABEL280;
continue;
for(J=N1;J<=NP;J++)
{
Q=A[L1*(N+1)+J];
R=A[I*(N+1)+J];
A[L1*(N+1)+J]=Q*CS+R*SN;
A[I*(N+1)+J]=R*CS-Q*SN;
}
}
//***************************************************************************
//************* TEST FOR CONVERGENCE *************
//***************************************************************************
LABEL290:
W=S[K];
if(L==K)
goto LABEL360;
//***************************************************************************
//***************** ORIGIN SHIFT *************
//***************************************************************************
X=S[L];
Y=S[K-1];
G=T[K-1];
H=T[K];
F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.0E0*H*Y);
G=sqrt(F*F+1.0E0);
if(F<0.0E0)
G=-G;
F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X;
//***************************************************************************
//************** QR STEP *************
//***************************************************************************
CS=1.0E0;
SN=1.0E0;
L1=L+1;
for(I=L1;I<=K;I++)
{
G=T[I];
Y=S[I];
H=SN*G;
G=CS*G;
W=sqrt(H*H+F*F);
T[I-1]=W;
CS=F/W;
SN=H/W;
F=X*CS+G*SN;
G=G*CS-X*SN;
H=Y*SN;
Y=Y*CS;
if(NV==0)
goto LABEL310;
for(J=1;J<=N;J++)
{
X=real(V[J*(N+1)+I-1]);
W=real(V[J*(N+1)+I]);
V[J*(N+1)+I-1]=COMPLEX_DOUBLE (X*CS+W*SN,0.0E0);
V[J*(N+1)+I]=COMPLEX_DOUBLE (W*CS-X*SN,0.0E0);
}
LABEL310:
W=sqrt(H*H+F*F);
S[I-1]=W;
CS=F/W;
SN=H/W;
F=CS*G+SN*Y;
X=CS*Y-SN*G;
if(NU==0)
goto LABEL330;
for(J=1;J<=N;J++)
{
Y=real(U[J*(N+1)+I-1]);
W=real(U[J*(N+1)+I]);
U[J*(N+1)+I-1]=COMPLEX_DOUBLE (Y*CS+W*SN,0.0E0);
U[J*(N+1)+I]=COMPLEX_DOUBLE (W*CS-Y*SN,0.0E0);
}
LABEL330:
if(N==NP)
// goto LABEL350;
continue;
for(J=N1;J<=NP;J++)
{
Q=A[(I-1)*(N+1)+J];
R=A[I*(N+1)+J];
A[(I-1)*(N+1)+J]=Q*CS+R*SN;
A[I*(N+1)+J]=R*CS-Q*SN;
}
}
//LABEL350:
T[L]=0.0E0;
T[K]=F;
S[K]=X;
goto LABEL220;
//***************************************************************************
//***************** CONVERGENCE *********************
//***************************************************************************
LABEL360:
if(W>=0.0E0)
// goto LABEL380;
continue;
S[K]=-W;
if(NV==0)
// goto LABEL380;
continue;
for(J=1;J<=N;J++)
V[J*(N+1)+K]=-V[J*(N+1)+K];
}
//***************************************************************************
//******************* SORT SINGULAR VALUES *********
//***************************************************************************
for(K=1;K<=N;K++)
{
G=-1.0E0;
J=K;
for(I=K;I<=N;I++)
{
if(S[I]<=G)
continue;
G=S[I];
J=I;
}
if(J==K)
// goto LABEL450;
continue;
S[J]=S[K];
S[K]=G;
if(NV==0)
goto LABEL410;
for(I=1;I<=N;I++)
{
Q=V[I*(N+1)+J];
V[I*(N+1)+J]=V[I*(N+1)+K];
V[I*(N+1)+K]=Q;
}
LABEL410:
if(NU==0)
goto LABEL430;
for(I=1;I<=N;I++)
{
Q=U[I*(N+1)+J];
U[I*(N+1)+J]=U[I*(N+1)+K];
U[I*(N+1)+K]=Q;
}
LABEL430:
if(N==NP)
// goto LABEL450;
continue;
for(I=N1;I<=NP;I++)
{
Q=A[J*(N+1)+I];
A[J*(N+1)+I]=A[K*(N+1)+I];
A[K*(N+1)+I]=Q;
}
}
//***************************************************************************
//*************** BACK TRANSFORMATION *************
//***************************************************************************
if(NU==0)
goto LABEL510;
for(KK=1;KK<=N;KK++)
{
K=N1-KK;
if(B[K]==0.0E0)
// goto LABEL500;
continue;
Q=-A[K*(N+1)+K]/abs(A[K*(N+1)+K]);
for(J=1;J<=NU;J++)
U[K*(N+1)+J]=Q*U[K*(N+1)+J];
for(J=1;J<=NU;J++)
{
Q=COMPLEX_DOUBLE (0.0E0,0.0E0);
for(I=K;I<=M;I++)
Q=Q+conjg(A[I*(N+1)+K])*U[I*(N+1)+J];
Q=Q/(abs(A[K*(N+1)+K])*B[K]);
for(I=K;I<=M;I++)
U[I*(N+1)+J]=U[I*(N+1)+J]-Q*A[I*(N+1)+K];
}
}
LABEL510:
if(NV==0)
return;
if(N<2)
return;
for(KK=2;KK<=N;KK++)
{
K=N1-KK;
K1=K+1;
if(C[K1]==0.0E0)
// goto LABEL560;
continue;
Q=-conjg(A[K*(N+1)+K1])/abs(A[K*(N+1)+K1]);
for(J=1;J<=NV;J++)
V[K1*(N+1)+J]=Q*V[K1*(N+1)+J];
for(J=1;J<=NV;J++)
{
Q=COMPLEX_DOUBLE (0.0E0,0.0E0);
for(I=K1;I<=N;I++)
Q=Q+A[K*(N+1)+I]*V[I*(N+1)+J];
Q=Q/(abs(A[K*(N+1)+K1])*C[K1]);
for(I=K1;I<=N;I++)
V[I*(N+1)+J]=V[I*(N+1)+J]-Q*conjg(A[K*(N+1)+I]);
}
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
for(int i = 0; i < M; i++)
for(int j = 0; j < N; j++)
matU(i, j) = U[(i+1)*(N+1)+(j+1)];
for (int i = 0; i < M; i++)
matU(i, 0) = -matU(i, 0);
for(int i = 0; i < M; i++)
for(int j = 0; j < N; j++)
matV(i, j) = V[(i+1)*(N+1)+(j+1)];
for (int i = 0; i < M; i++)
matV(i, 0) = -matV(i, 0);
delete[] A;
delete[] U;
delete[] V;
return;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -