rheolef  6.3
tensor_eig_tst.cc
Go to the documentation of this file.
1 #include "rheolef/tensor.h"
2 #include "rheolef/tensor-eig3.h"
3 #include "rheolef/tensor-svd-algo.h"
4 #include "rheolef/iorheo.h"
5 using namespace rheolef;
6 using namespace std;
7
8 int eig_check (const tensor& a, size_t dim) {
9  tensor q;
10  point d = a.eig (q,dim);
11  tensor err_a = a - q*diag(d)*trans(q);
12  Float err = norm(err_a);
13  cout << setprecision(15) << matlab
14  << "a = " ; a.put(cout,dim); cout << endl
15  << "d = " ; diag(d).put(cout,dim); cout << endl
16  << "q = " ; q.put(cout,dim); cout << endl
17  << "err = norm(a-q*d*q')" << endl
18  << "%err =" << err << endl
19  ;
20  return (err < 1e-10) ? 0 : 1;
21 }
22 int eig2x2_tst () {
23  tensor a;
24  a(0,0) = 1;
25  a(1,0) = 2;
26  a(0,1) = 2;
27  a(1,1) = 3;
28  return eig_check (a, 2);
29 }
30 int eig3x3_tst () {
31  tensor a;
32  a(0,0) = 1;
33  a(1,0) = 2;
34  a(2,0) = 3;
35  a(0,1) = 2;
36  a(1,1) = 3;
37  a(2,1) = 4;
38  a(0,2) = 3;
39  a(1,2) = 4;
40  a(2,2) = 5;
41  return eig_check (a, 3);
42 }
43 int main(int argc, char**argv) {
44  int status = 0;
45  status |= eig2x2_tst() ;
46  status |= eig3x3_tst() ;
47  return status;
48 }
49