rheolef  6.5
tensor_svd_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 svd_check (const tensor& a, size_t dim) {
9  tensor u,v;
10  point s = a.svd (u,v,dim);
11  tensor err_a = a - u*diag(s)*trans(v);
12  Float err = norm(err_a);
13  cout << setprecision(15) << matlab
14  << "a = " ; a.put(cout,dim); cout << endl
15  << "s = " ; diag(s).put(cout,dim); cout << endl
16  << "u = " ; u.put(cout,dim); cout << endl
17  << "v = " ; v.put(cout,dim); cout << endl
18  << "err = norm(a-u*s*v')" << endl
19  << "%err =" << err << endl
20  ;
21  return (err < 1e-10) ? 0 : 1;
22 }
23 int svd2x2_tst () {
24  tensor a;
25  a(0,0) = 1;
26  a(1,0) = 2;
27  a(0,1) = 3;
28  a(1,1) = 4;
29  return svd_check (a, 2);
30 }
31 int svd3x3_tst () {
32  tensor a;
33  a(0,0) = 1;
34  a(1,0) = 2;
35  a(2,0) = 3;
36  a(0,1) = 4;
37  a(1,1) = 5;
38  a(2,1) = 6;
39  a(0,2) = 7;
40  a(1,2) = 8;
41  a(2,2) = 9;
42  return svd_check (a, 3);
43 }
44 int main(int argc, char**argv) {
45  int status = 0;
46  status |= svd2x2_tst() ;
47  status |= svd3x3_tst() ;
48  return status;
49 }
50