```
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
#define SIZE 50
#define REAL double
#define SQRT sqrt
#define ABS fabs
//complex number
typedef struct {
REAL R,I;
} COMPLEX;
typedef COMPLEX CMat[SIZE][SIZE];
int I,it,J,M,N;
REAL dta;
CMat A, VX;
REAL R[SIZE];
REAL Sqr(REAL a) {
return a*a;
}
//Absolute value of a complex number
REAL CABS(COMPLEX C) {
return (SQRT(Sqr(C.R)+Sqr(C.I)));
}
//Add two complex numbers
void CADD(COMPLEX c1, COMPLEX c2, COMPLEX *c3) {
c3->R=c1.R+c2.R; c3->I=c1.I+c2.I;
}
//Substract two complex numbers
void CDIF(COMPLEX c1, COMPLEX c2, COMPLEX *c3) {
c3->R=c1.R-c2.R; c3->I=c1.I-c2.I;
}
//Multiply two complex numbers
void CMUL(COMPLEX c1, COMPLEX c2, COMPLEX *c3) {
c3->R=c1.R*c2.R - c1.I*c2.I;
c3->I=c1.R*c2.I + c1.I*c2.R;
}
//Return conjugate of a complex number
void CONJ(COMPLEX C, COMPLEX *C1) {
C1->R=C.R; C1->I=-C.I;
}
//Multiply a complex number by a real number
void CPRO(REAL alpha, COMPLEX C, COMPLEX *C1) {
C1->R=alpha*C.R; C1->I=alpha*C.I;
}
void EPHJ(REAL dta, int M, int N, CMat A, int *it, REAL *R, CMat VX) {
int I,J,K,K0,L,L0;
REAL delta,s,s0,t0,t1,w0;
COMPLEX c0,c1,c2,c3,u0,u1,z0,z1;
z0.R=0.0; z0.I=0.0; z1.R=1.0; z1.I=0.0;
for (I=1; I<=N; I++)
for (J=1; J<=N; J++)
if (I==J) VX[I][J]=z1;
else VX[I][J]=z0;
*it=-1; L=1;
do {
s=0.0;
for (I=1; I<N; I++)
for (J=I+1; J<=N; J++) {
t0=CABS(A[I][J]);
if (t0 > s) {
s=t0; K0=I; L0=J;
}
}
if (s==0.0) *it=1;
else {
delta=Sqr(A[L0][L0].R-A[K0][K0].R)+4.0*Sqr(CABS(A[K0][L0]));
t0=A[L0][L0].R-A[K0][K0].R + SQRT(delta);
t1=A[L0][L0].R-A[K0][K0].R - SQRT(delta);
if (ABS(t0) >= ABS(t1)) w0=t0; else w0=t1;
s0=ABS(w0)/SQRT(Sqr(w0)+4.0*Sqr(CABS(A[K0][L0])));
t0=2.0*s0/w0; CPRO(t0,A[K0][L0],&c0);
CONJ(c0,&c1);
for (I=1; I<K0; I++) {
u0=A[I][K0];
CMUL(c0,u0,&c2); CPRO(s0,A[I][L0],&c3); CADD(c2,c3,&A[I][K0]);
CMUL(c1,A[I][L0],&c2); CPRO(s0,u0,&c3); CDIF(c2,c3,&A[I][L0]);
}
for (K=K0+1; K<L0; K++) {
u0=A[K0][K];
CMUL(c1,u0,&c2); CONJ(A[K][L0],&u1); CPRO(s0,u1,&c3);
CADD(c2,c3,&A[K0][K]);
CMUL(c1,A[K][L0],&c2); CONJ(u0,&u1); CPRO(s0,u1,&c3);
CDIF(c2,c3,&A[K][L0]);
}
for (J=L0+1; J<=N; J++) {
u0=A[K0][J];
CMUL(c1,u0,&c2); CPRO(s0,A[L0][J],&c3); CADD(c2,c3,&A[K0][J]);
CMUL(c0,A[L0][J],&c2); CPRO(s0,u0,&c3); CDIF(c2,c3,&A[L0][J]);
};
t0=A[K0][K0].R;
t1=4.0*Sqr(s0*CABS(A[K0][L0]))/w0;
A[K0][K0].R=Sqr(CABS(c0))*t0 + t1+Sqr(s0)*A[L0][L0].R;
A[L0][L0].R=Sqr(s0)*t0 - t1+Sqr(CABS(c0))*A[L0][L0].R;
A[K0][L0]=z0;
for (I=1; I<=N; I++) {
u0=VX[I][K0];
CMUL(c0,u0,&c2); CPRO(s0,VX[I][L0],&c3); CADD(c2,c3,&VX[I][K0]);
CMUL(c1,VX[I][L0],&c2); CPRO(s0,u0,&c3); CDIF(c2,c3,&VX[I][L0]);
}
t0=0.0;
for (I=1; I<N; I++)
for (J=I+1; J<=N; J++)
t0 += Sqr(CABS(A[I][J]));
s=2.0*t0;
if (s < dta) *it=1; else L++;
} //else
} while (L<=M && *it!=1);
if (*it==1)
for (I=1; I<=N; I++) R[I]=A[I][I].R;
}
void CDIV(double AR, double AI, double BR, double BI, double *ZR, double *ZI) {
double YR,YI,W;
YR=BR;
YI=BI;
if (ABS(YR) > ABS(YI)) {
W=YI/YR;
YR=W*YI+YR;
*ZR=(AR+W*AI)/YR;
*ZI=(AI-W*AR)/YR;
}
else {
W=YR/YI;
YI=W*YR+YI;
*ZR=(W*AR+AI)/YI;
*ZI=(W*AI-AR)/YI;
}
}
void NORMAL(int N, REAL *R, CMat Z) {
REAL Tr[SIZE], Ti[SIZE];
double ZM,Z1,Z2,VR;
COMPLEX V;
int IM,J,K;
// SORT SOLUTIONS IN ASCENDING ORDER
for (J=2; J<=N; J++) {
VR=R[J];
for (K=1; K<=N; K++) {
Tr[K]=Z[K][J].R;
Ti[K]=Z[K][J].I;
}
for (I=J-1; I>0; I--) {
if (ABS(R[I]) <= ABS(VR)) goto e5;
R[I+1]=R[I];
for (K=1; K<=N; K++) {
Z[K][I+1].R=Z[K][I].R;
Z[K][I+1].I=Z[K][I].I;
}
}
I=0;
e5: R[I+1]=VR;
for (K=1; K<=N; K++) {
Z[K][I+1].R=Tr[K];
Z[K][I+1].I=Ti[K];
}
}
// NORMALIZE WITH RESPECT TO BIGGEST ELEMENT
for (J=N; J>0; J--) {
ZM = 0.0;
for (I=1; I<=N; I++) {
V.R=Z[I][J].R;
V.I=Z[I][J].I;
Z1=CABS(V);
if (Z1 >= ZM) {
IM = I;
ZM = Z1;
}
}
Z1 = Z[IM][J].R;
Z2 = Z[IM][J].I;
for (I=1; I<=N; I++) {
CDIV(Z[I][J].R,Z[I][J].I,Z1,Z2,&V.R,&V.I);
Z[I][J].R = V.R;
Z[I][J].I = V.I;
}
}
} //NORMAL()
//Read from input file complex hermitian matrix
void LME() {
int i, j;
FILE *fp;
fp=fopen("matrix.dat","r");
fscanf(fp,"%d", &N);
for (i=1; i<=N; i++)
for (j=1; j<=N; j++)
fscanf(fp,"%lf %lf", &A[i][j].R, &A[i][j].I);
fclose(fp);
};
int main() {
printf("\n");
printf(" Precision: "); scanf("%lf", &dta);
printf(" Max. number of iterations: "); scanf("%d", &M);
printf("\n");
LME();
EPHJ(dta, M, N, A, &it, R, VX);
NORMAL(N, R, VX);
if (it==-1)
printf(" No convergence.\n");
else // display results
for (J=1; J<=N; J++) {
cout << R[J] << endl;
}
return 0;
}
```