C++

Forum

 
sorting eigenvalues problem

Katka01 (6)
Hello,
I have a program, which returns eigenvalues and eigenvectors in form d[i] and A[i][j]. After getting the values, I need to sort them by increasing value. I found a function and a piece of code online, which I modified a bit and it works great for some matrices. I'm acually studying eigenvalues of peturbed Hamiltonian. I won't go into details, but to obtain these matrices, whose eigensystems I'm looking for, I use several different methods. Therefore, I obtain different matrices, although their elements are of similar orders of magnitude. One has elements that go as low as 10^(-4) and - I'm not sure this is the reason - the sorting algorithm instead of sorted eigenvalues, returns exponentially falling, negative numbers, which have nothing to do with original numbers. Can anyone take a look at the sorting algorithm and tell me what is wrong with it?


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
int comp (const void *p1, const void *p2)
{
  int *a, *b;
  a = (int*) p1;
  b = (int*) p2;
  
  if (energies[*a] > energies[*b])
    return 1;
  return -1;
}

int index[NN];
    for (int u=0; u<NN;u++)
	{index[u]=u;}
	energies=d;
	qsort(index,NN,sizeof(index[0]), comp);
	for (int i=0; i<NN; i++)
		{ eigsorted[i]=d[index[i]];
	}
ne555 (10692)
¿why do you need the `index' array? ¿why don't just sort `d' directly?
http://www.cplusplus.com/reference/algorithm/sort/
std::sort( d, d+NN );


Try to post a minimal example that does reproduce your issue
Katka01 (6)
Here's almost all of it. I didn't copy the routines tqli and tred2, but if necessary, I can add them. Although I think the problem is in that sorting algorithm. If I print out the eigenvalues before they are sorted (vector d) for NN=30, I get
4 1.19293928e+000
5 1.03285686e+000
6 6.12672163e-001
7 1.20634833e+000
8 1.02091839e+000
9 5.37507933e-001
10 1.15730840e+000
11 8.81965745e-001
12 1.99320249e-001
13 1.04903576e+000
14 5.71159352e-001
15 1.13790663e+000
16 8.04462364e-001
17 1.19675733e+000
18 9.54570019e-001
19 3.26299214e-001
20 1.04993919e+000
21 5.52807343e-001
22 1.11347989e+000
23 7.18139262e-001
24 1.15979695e+000
25 8.45560415e-001
26 1.19704384e+000
27 9.41135224e-001
28 2.85617174e-001
29 1.01227511e+000

If I print them after the sorting algorithm (vector eigsorted[i]), they are
4 1.19293928e+000
5 1.03285686e+000
6 6.12672163e-001
7 4.82177119e-001
8 4.81451269e-001
9 4.75249808e-001
10 -3.96904312e-001
11 -1.72016870e+000
12 -3.42807829e+000
13 -5.52093282e+000
14 -8.00206466e+000
15 -1.08756008e+001
16 -1.41457677e+001
17 -1.78166413e+001
18 -2.18920586e+001
19 -2.63755943e+001
20 -3.12705644e+001
21 -3.65800418e+001
22 -4.23068750e+001
23 -4.84537080e+001
24 -5.50229997e+001
25 -6.20170409e+001
26 -6.94379708e+001
27 -7.72877906e+001
28 -8.55683768e+001
29 -9.42814923e+001


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
 
#include <iostream>
#include <cmath>
#include <string>
#include <stdlib.h>

void  tred2(double **a, int n, double *d, double *e);
void  tqli(double *d, double *e, int n, double **z);

// defintions of matrices (C/C++ style -- null determined)
template <class T> T** matrix(long nrow, long ncol);
template <class T> void free_matrix(T **m);
double *energies;
int comp (const void *p1, const void *p2);
// **************************** Glavni program ********************************
#include <iomanip>
#include <fstream>
#include <array>

int main(){
	using namespace std;

double lambda=0.5;
int NN=0;
     ofstream output("A_LV2_convergence_lambda_" + std::to_string(lambda) + ".txt");
     streambuf* saved_buffer=cout.rdbuf();
     cout.rdbuf(output.rdbuf());
     cout.setf(ios::scientific,ios::floatfield);
     cout << scientific << setprecision(8);

    	for (NN=4; NN<100; NN=NN+1)
    	{
    	long n = NN, i, j;
        double **A = matrix <double> (n,n), // 3x3
          *d = new double [n],
	           *e = new double [n],
	           *eigsorted = new double [n];
	           
        long double H0[NN][NN];
         for(int z3 = 0; z3<NN; z3++)
        {for (int w3 = 0; w3<NN; w3++){
        	double delta=0;
        	if(w3==z3){delta=1.0;}
        	H0[z3][w3]= (0.5+w3)*delta;
		}
   		}
    ////////////////////////////////////////////////////////////////////////////////
        long double qij[NN][NN];
 			for(int z1=0;z1<NN;z1++)
        	{for (int w1=0;w1<NN;w1++){
			qij[z1][w1]=0;
        	if(w1==z1+1){ if (z1==0){qij[z1][w1]=0;}
						  else {qij[z1][w1] = 0.5*sqrt(w1*2.0); }
					}
        	if(w1==z1-1){ qij[z1][w1] = 0.5*sqrt(2.0*w1+2.0); }
		}
    }
     /////////////////////////////////////////////////////////
     long double H[NN][NN];
     for(int i1=0; i1<NN; i1++)
     {for (int j1=0; j1<NN; j1++)
     {H[i1][j1] = H0[i1][j1] + lambda*pow(qij[i1][j1], 4.0);
     	 }
     }

   for(int i1=0; i1<NN; i1++)
     {for (int j1=0; j1<NN; j1++)
     {A[i1][j1] = H[i1][j1];}
 }
 
  tred2(A, n, d, e);
  tqli (d, e, n, A);
  
	int index[NN];
    for (int u=0; u<NN;u++)
	{index[u]=u;}
	energies=d;
	qsort(index,NN,sizeof(index[0]), comp);
	for (int i=0; i<NN; i++)
		{ eigsorted[i]=d[index[i]];
	}
  cout << NN << "      "<< setw(8) <<eigsorted[1] << endl;

  delete [] d;
  delete [] e;
  delete [] lastvr;
  free_matrix <double> (A);
}
  return EXIT_SUCCESS;
}
// ******************************* Rutines **********************************
template <class T> T** matrix (long nrow, long ncol) {

  T **m = new T* [nrow];
  m[0] = (T *) new char [nrow*ncol*sizeof(T)];
  for( long i = 1; i < nrow; i++) m[i] = m[i-1] + ncol;

  return m;
}

template <class T> void free_matrix(T **m) {
  delete [] (char*) m[0];
  delete [] m;
}

inline double sqr(double f){ return f*f; }

inline double pythag(double a, double b) {
  double absa = std::abs(a), absb = std::abs(b);
  if (absa > absb) return absa*sqrt(1.0 + sqr(absb/absa));
  return (absb == 0.0 ? 0.0 : absb*sqrt(1.0 + sqr(absa/absb)));
}

int comp (const void *p1, const void *p2)
{
  int *a, *b;
  a = (int*) p1;
  b = (int*) p2;
  if (energies[*a] > energies[*b])
    return 1;
  return -1;
}

void tred2(double **a, int n, double *d, double *e) {
...
}

#define sign_(a,b) ((b) >= 0.0 ? std::abs(a) : -std::abs(a))
void tqli(double *d, double *e, int n, double **z) {
...
}
#undef sign_
Last edited on
keskiverto (10402)
Those numbers look correctly sorted (descended).
  1.0
  0.6
  0.4
 -0.3
 -1.7
 -5.5
-10.8

Use std::fixed and precision to create more intuitive format.
ne555 (10692)
but the values are crap.
You may be accessing out of bounds in another place.

Post the missing functions, so we can compile and run your code.
Also, learn to indent.
Katka01 (6)
The eigenvalues, calculated by tqli and tred2 are correct. For a matrix 30x30 I posted them above. But, when I put them through sorting algorithm, they change completely and become negative. It has nothing to do with tqli or tred2. Besides, for two other matrices (let's call it B and C, and I'll say A to the one I posted above), obtained by a bit different functions in "for" loops for "qij" and "H", the sorting thingy works perfectly. I also tried with a simpler sorting code that I found online, 
1
2
 
	qsort(d,n,sizeof(double),
         (int(*)(const void *,const void *))comp);

where
1
2
3
4
5
6
7
8
9
10
11
12
 
/*
  The function   int comp()                  
  is a utility function for the library function qsort()
  to sort double numbers after increasing values.
*/       

int comp(const double *val_1, const double *val_2)
{
  if((*val_1) <= (*val_2))       return -1;
  else  if((*val_1) > (*val_2))  return +1;
  else                     return  0; 
} // End: function comp()


It works fine with B and C, but not with A.
I'm posting tqli and tred2 below, although I really don't think they matter here.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
 
void tred2(double **a, int n, double *d, double *e) {

  int l, k, j, i;
  double scale, hh, h, g, f;
  for (i = n-1; i > 0; i--) {
    l = i-1;
    h = scale = 0.0;
    if (l > 0) {
      for (k = 0; k <= l; k++)
        scale += std::abs(a[i][k]);
	if (scale == 0.0)
          e[i] = a[i][l];
	else {
	  for (k = 0; k <= l;k++) {
	    a[i][k] /= scale;
	    h += a[i][k]*a[i][k];
	}
	f = a[i][l];
	g = (f >= 0.0 ? -sqrt(h) : sqrt(h));
	e[i] = scale*g;
	h -= f*g;
	a[i][l] = f-g;

	f = 0.0;
	for (j = 0; j <= l; j++) {
  	  a[j][i] = a[i][j]/h;
	  g = 0.0;

	  for (k = 0; k <= j; k++) g += a[j][k]*a[i][k];

	  for (k=j+1;k<=l;k++) g += a[k][j]*a[i][k];

	  e[j] = g/h;
	  f += e[j]*a[i][j];
	}

	hh = f/(h+h);
	for (j = 0; j <= l; j++) {
	  f = a[i][j];
	  e[j] = g = e[j] - hh*f;

	  for (k = 0; k <= j; k++) a[j][k] -= (f*e[k]+g*a[i][k]);
	}
      }
    } else e[i]=a[i][l];

    d[i]=h;
  }

  d[0]=0.0;
  e[0]=0.0;

  for (i = 0; i < n; i++) {
    l = i-1;
    if (d[i]) {
      for (j = 0; j <= l; j++) {
        g = 0.0;
        for (k = 0; k <= l; k++) g += a[i][k]*a[k][j];
        for (k = 0; k <= l; k++) a[k][j] -= g*a[k][i];
      }
    }
    d[i] = a[i][i];
    a[i][i] = 1.0;
    for (j = 0; j <= l; j++) a[j][i] = a[i][j] = 0.0;
  }
}

#define sign_(a,b) ((b) >= 0.0 ? std::abs(a) : -std::abs(a))

void tqli(double *d, double *e, int n, double **z) {

  int m, l, iter, i, k;

  double s, r, p, g, f, dd, c, b;

  for (i =1; i < n; i++) e[i-1] = e[i];

  e[n-1] = 0.0;

  for (l = 0; l < n; l++) {

    iter = 0;

    do {
      for (m = l; m < n - 1; m++) {
        dd = std::abs(d[m]) + std::abs(d[m+1]);
        if ((double)(std::abs(e[m]) + dd) == dd) break;
      }

      if (m != l) {

        if (iter++ == 30) {
          std::cerr << "Too many iterations in tqli";
	        exit(EXIT_FAILURE);
        }

        g = (d[l+1] - d[l])/(2.0*e[l]);
        r = pythag(g,1.0);
        g = d[m] - d[l] + e[l]/(g + sign_(r,g));
        s = c = 1.0;
        p = 0.0;

        for (i = m-1; i >= l; i--) {
          f = s*e[i];
          b = c*e[i];
      	  e[i+1] = (r = pythag(f,g));

      	  if (r == 0.0) {
      	    d[i+1] -= p;
      	    e[m] = 0.0;
      	    break;
      	  }
      	  s = f/r;
      	  c = g/r;
      	  g = d[i+1] - p;
      	  r = (d[i] - g)*s + 2.0*c*b;
      	  d[i+1] = g + (p = s*r);
      	  g = c*r - b;

          for (k = 0; k < n; k++) {
      	    f = z[k][i+1];
      	    z[k][i+1] = s*z[k][i] + c*f;
      	    z[k][i] = c*z[k][i] - s*f;
      	  }
        }
        if (r == 0.0 && i >= l) continue;

        d[l] -= p;
        e[l] = g;
        e[m] = 0.0;
      }
    } while (m != l);
  }
}

#undef sign_
keskiverto (10402)
but the values are crap.

Doh! Missed that.

Ok, there is an array of integer values; continuous range of values [0..NN[
The array is reordered. It should still contain the same set of values.

I see no code that prints the d.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
int index[NN];
for ( int u=0; u<NN; u++ ) {
  index[u] = u;
}
std::cout << "Before sort\n";
for ( int i=0; i<NN; i++ ) {
  std::cout << i << " # " << d[i] << " # " << index[i] << " # " << d[ index[i] ] << '\n';
}
energies=d;
qsort(index,NN,sizeof(index[0]), comp);
std::cout << "After sort\n";
for ( int i=0; i<NN; i++ ) {
  std::cout << i << " # " << d[i] << " # " << index[i] << " # " << d[ index[i] ] << '\n';
}

That is trivial debugging.


There are other questions though.

Why the use of qsort() and preprocessor macro sign_? The same program does use templates and I/O streams. It is understandable that you have got part of the codebase from other sources, but you should thrive towards consistency and I'll wouch for std::sort.


Finally, the big one:

The code has this "beauty":
1
2
3
 
int NN=0;
// code
int index[NN];

That is not standard C++. Some compiler extensions do support VLAs (Variable Length Arrays), but the proposal was omitted from C++14. I would not rely on such feature. I would be against such feature, if I were in the standardization committee.


Luckily, this instance of the problem is easy to replace:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
#include <vector>    // std::vector
#include <numeric>   // std::iota
#include <algorithm> // std::sort


using std::begin;
using std::end;
std::vector<size_t> index( NN );
std::iota( begin(index), end(index), 0 );
std::sort( begin(index), end(index),
           [d](auto a, auto b) { return d[b] < d[a]; } );
std::cout << "After sort\n";
for ( int i=0; i<NN; ++i ) {
  std::cout << i << " # " << d[i] << " # " << index[i] << " # " << d[ index[i] ] << '\n';
}
Katka01 (6)
Thank you very much for your replies. I'm fairly new to programming and I'm learning from mistakes. I tried to use the sorting algorithm that keskiverto proposed, but I immediately ran into a problem. I'm using Windows environment (and Dev C++) and I think that C++11 is the newest standard here, which is causing a problem in 'auto' function. How can I delcare the variables a and b without using 'auto'?
Last edited on
keskiverto (10402)
Reference documentation:
http://www.cplusplus.com/reference/algorithm/sort/
comp: "two elements in the range as arguments"

The elements are elements of the index and index is a std::vector.
http://www.cplusplus.com/reference/vector/vector/
A vector of size_t elements.


However, Dev C++ is not the compiler. Dev C++ is an integrated development environment (IDE) that hides the real compiler under its hood. The latest Orwell Dev-C++ seems to include TDM-GCC 4.8.1 by default, but they have a "latest tested compilers" page:
http://sourceforge.net/projects/orwelldevcpp/files/Compilers/

Therefore, you should have access to C++14 awesomeness.
ne555 (10692)
> the sorting algorithm instead of sorted eigenvalues,
> returns exponentially falling, negative numbers,
> which have nothing to do with original numbers.
That's incorrect. The `d' array does have those negative numbers before sorting.

> The eigenvalues, calculated by tqli and tred2 are correct. For a matrix 30x30 I posted them above
This is what I'm getting for `d' when `NN=30'
$1 = {0.48849536172265384, 0.43616846105726848, 1.2605725444749634, 2.0262631756688285, 3.287341476478538,
	5.0979053270402677, -2.1367633773986303, 7.5473079990682672, 10.736061858016765, -6.2792684670322441,
	14.774117624854751, 19.785585011255154, -12.051740476170435, 25.914884237778789, -19.653377436806405,
	33.335055286150862, -29.388144614105325, 42.259599245705452, 52.960421939402977, -41.6995335525718,
	65.796869972943043, -57.256491125037364, 81.266424125527436, -77.150095589014455, 100.1021606846265,
	123.48693427420736, -103.4287962822827, 153.63194961599186, -141.22126782512993, 196.07136052357737}
as you can see, the negative numbers were already there.


Check that the matrix is assembled correctly
Review your eigen functions

Also, comment your code.


EDIT: used octave to compute the eigenvalues. They are all complex numbers of norm=1
Last edited on
Topic archived. No new replies allowed.