$title Test eigenvector utility (EIGVEC01,SEQ=410)

$onText

Eigenvector example.

Contributor: Erwin Kalvelagen, October 2008.

octave:1> a = [1 2 4 7 11; 2 3 5 8 12; 4 5 6 9 13; 7 8 9 10 14; 11 12 13 14 15]
a =

   1   2   4   7  11
   2   3   5   8  12
   4   5   6   9  13
   7   8   9  10  14
  11  12  13  14  15

octave:2> eig(a)
ans =

   -8.464425
   -1.116317
   -0.512109
   -0.027481
   45.120332

octave:3> [e1,e2] = eig(a)
e1 =

e2 =

   -8.46442    0.00000    0.00000    0.00000    0.00000
    0.00000   -1.11632    0.00000    0.00000    0.00000
    0.00000    0.00000   -0.51211    0.00000    0.00000
    0.00000    0.00000    0.00000   -0.02748    0.00000
    0.00000    0.00000    0.00000    0.00000   45.12033

$offText

set i /i1*i5/;
alias (i,j);

table a(i,j)
      i1   i2   i3   i4   i5
i1     1    2    4    7   11
i2     2    3    5    8   12
i3     4    5    6    9   13
i4     7    8    9   10   14
i5    11   12   13   14   15

;

Table expected(i,*)
           val          i1          i2          i3          i4          i5
i1   -8.464425   0.5550905  -0.2642556   0.2892854   0.6748602   0.2879604
i2   -1.116317   0.4820641  -0.2581518   0.2196341  -0.7349311   0.3355726
i3   -0.512109   0.2865066   0.2159261  -0.8437897   0.0411896   0.3970041
i4   -0.027481  -0.0992784   0.7711236   0.3943678   0.0055409   0.4898525
i5   45.120332  -0.6062562  -0.4714561  -0.0238286   0.0520829   0.6378888
;

parameter eval(i) 'eigenvalues';
parameter evec(i,j) 'eigenvectors';

execute_unload 'a.gdx', i, a;
executeTool.checkErrorLevel 'linalg.eigenvector i a eval evec -gdxin=a.gdx -gdxout=b.gdx';
execute_load 'b.gdx', eval, evec;
evec(i,j)$(evec('i1','i1')<0) = -evec(i,j);
abort$(sum(i$(abs(eval(i)-expected(i,'val'))>1e-4),1)) 'Wrong Eigenvalue', eval, expected;
abort$(sum((i,j)$(abs(evec(i,j)-expected(i,j))>1e-4),1)) 'Wrong Eigenvector', evec, expected;

option clear=eval, clear=evec;
executeTool.checkErrorLevel 'linalg.eigenvector i a eval evec';
evec(i,j)$(evec('i1','i1')<0) = -evec(i,j);
abort$(sum(i$(abs(eval(i)-expected(i,'val'))>1e-4),1)) 'Wrong Eigenvalue', eval, expected;
abort$(sum((i,j)$(abs(evec(i,j)-expected(i,j))>1e-4),1)) 'Wrong Eigenvector', evec, expected;
*
* only lower triangular part of A is used
*
table a2(i,j)
      i1   i2   i3   i4   i5
i1     1
i2     2    3
i3     4    5    6
i4     7    8    9   10
i5    11   12   13   14   15

;

execute_unload 'a.gdx', i, a2;
executeTool.checkErrorLevel 'linalg.eigenvector i a2 eval evec -gdxin=a.gdx -gdxout=b.gdx';
execute_load 'b.gdx', eval, evec;
evec(i,j)$(evec('i1','i1')<0) = -evec(i,j);
abort$(sum(i$(abs(eval(i)-expected(i,'val'))>1e-4),1)) 'Wrong Eigenvalue', eval, expected;
abort$(sum((i,j)$(abs(evec(i,j)-expected(i,j))>1e-4),1)) 'Wrong Eigenvector', evec, expected;

option clear=eval, clear=evec;
executeTool.checkErrorLevel 'linalg.eigenvector i a2 eval evec';
evec(i,j)$(evec('i1','i1')<0) = -evec(i,j);
abort$(sum(i$(abs(eval(i)-expected(i,'val'))>1e-4),1)) 'Wrong Eigenvalue', eval, expected;
abort$(sum((i,j)$(abs(evec(i,j)-expected(i,j))>1e-4),1)) 'Wrong Eigenvector', evec, expected;

*
* only upper triangular part of A is used
*
Parameter a3(i,j); a3(i,j) = a2(j,i)

execute_unload 'a.gdx', i, a3;
executeTool.checkErrorLevel 'linalg.eigenvector i a3 eval evec -gdxin=a.gdx -gdxout=b.gdx';
execute_load 'b.gdx', eval, evec;
evec(i,j)$(evec('i1','i1')<0) = -evec(i,j);
abort$(sum(i$(abs(eval(i)-expected(i,'val'))>1e-4),1)) 'Wrong Eigenvalue', eval, expected;
abort$(sum((i,j)$(abs(evec(i,j)-expected(i,j))>1e-4),1)) 'Wrong Eigenvector', evec, expected;

option clear=eval, clear=evec;
executeTool.checkErrorLevel 'linalg.eigenvector i a3 eval evec';
evec(i,j)$(evec('i1','i1')<0) = -evec(i,j);
abort$(sum(i$(abs(eval(i)-expected(i,'val'))>1e-4),1)) 'Wrong Eigenvalue', eval, expected;
abort$(sum((i,j)$(abs(evec(i,j)-expected(i,j))>1e-4),1)) 'Wrong Eigenvector', evec, expected;