tensor_eig_tst.cc

Go to the documentation of this file.

#include "rheolef/tensor.h"
#include "rheolef/tensor-eig3.h"
#include "rheolef/tensor-svd-algo.h"
#include "rheolef/iorheo.h"
using namespace rheolef;
using namespace std;

int eig_check (const tensor& a, size_t dim) {
    tensor q;
    point d = a.eig (q,dim);
    tensor err_a = a - q*diag(d)*trans(q);
    Float err = norm(err_a);
    cout << setprecision(15) << matlab
         << "a = " ; a.put(cout,dim); cout << endl
         << "d = " ; diag(d).put(cout,dim); cout << endl
         << "q = " ; q.put(cout,dim); cout << endl
         << "err = norm(a-q*d*q')" << endl
         << "%err =" << err << endl
        ;
    return (err < 1e-10) ? 0 : 1;
}
int eig2x2_tst () {
    tensor a;
    a(0,0) = 1;
    a(1,0) = 2;
    a(0,1) = 2;
    a(1,1) = 3;
    return eig_check (a, 2);
}
int eig3x3_tst () {
    tensor a;
    a(0,0) = 1;
    a(1,0) = 2;
    a(2,0) = 3;
    a(0,1) = 2;
    a(1,1) = 3;
    a(2,1) = 4;
    a(0,2) = 3;
    a(1,2) = 4;
    a(2,2) = 5;
    return eig_check (a, 3);
}
int main(int argc, char**argv) {
  int status = 0;
  status |= eig2x2_tst() ;
  status |= eig3x3_tst() ;
  return status;
}