Member Avatar for begyu

Hi all,

I have a .cpp code which reads a matrix from a .dat file and compute their eigenvalues. But I would like to read more matrices from the file not only one, and get the eigenvalues in the similar way.
Please help to solve this, I stuck at this point (i.e. how to read more than one matrix from the file).

The meaning of the content of the matrix.dat file is the following:

2 <-- this is the size of the NxN matrix, N=2
1 0 0 0 <-- (1 0) means the real and imaginary part of a matrix element 1,1;
                 (0 0) means the real and imaginary part of a matrix element 1,2;
0 0 1 0 <-- (0 0) means the real and imaginary part of a matrix element 2,1;
                 (1 0) means the real and imaginary part of a matrix element 2,2;
#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;
}
2                                                                  
1 0 0 0
0 0 1 0
2                                                                  
2 0 0 0
0 0 2 0
2
3 0 0 0 
0 0 3 0
2
4 0 0 0 
0 0 4 0
2
5 0 0 0 
0 0 5 0
  • 4 Contributors
  • 4 Replies
  • 273 Views
  • 11 Hours Discussion Span
  • Latest Post Latest Post by Lerner
Member Avatar for Freaky_Chris

Would it not just mean adding some form of a break in the text file?

For example enter an 'x' after the first matrix, then the program knowns it finished and can call again for the next.

Maybe i missed something?

Member Avatar for Sci@phy

You don't even need that, since you have a number before your matrix that tells you how big it is.
Simply add second number after that, do some sort of while/for/do-while/goto loop (ok, maybe only while or for :) ) and it should work.
If you want more help, post your attempt (that is, attempt to read more matrices, not just one)

Member Avatar for begyu

Could you post what you are thinkig about? I tried this thing but does not work.

Please, post here the code.

Member Avatar for Lerner

//algorhithm to read file containing matrix data where size of matrix (which will always be a square matrix) precedes matrix data. There may be data for 1 or more matrixes in the file.

int sizeofMatrix, i , j;
ifstream fin("myfile.dat");

//check if file open and continue only if it is

while(fin >> sizeofMatrix)
{
   //declare memory for matrix here

   //read in matrix here

   //calculate eigenvalue of matrix here
}
//fin is now in failed state either because it read in the entire file or because of an error in the file.  Find out why.  If file read successfully completed then move on to next part of program, else indicate unable to read file.
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.