AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

matrix_test.cpp

#include <animal/array.h>
#include <animal/vec2.h>
#include <animal/matrix.h>
#include <animal/linear.h>
#include <animal/io.h>

using std::cout;
using std::cin;
using std::endl;

const int nr = 2;
const int nc = 3;
const int nr2 = 3;
const int nc2 = 2;

// typedef double Vec[nc];
// typedef double Vec2[nr];
// typedef double Mat[nr][nc];
// typedef double Mat2[nr2][nc2];
// typedef double Mat3[nr][nc2];
// typedef double Mat4[nc][nr2];
// typedef double Mat5[nc][nc];   // square matrix

typedef animal::Array<nc,double> Vec;
//typedef animal::Array<animal::Array<double,nc>,nr> Mat;
//typedef animal::Array<animal::Array<double,nc2>,nr2> Mat2;
//typedef animal::Array<animal::Array<double,nc2>,nr> Mat3;
//typedef animal::Array<animal::Array<double,nr2>,nc> Mat4;
//typedef animal::Array<animal::Array<double,nc>,nc> Mat5;
typedef animal::Matrix<double> Mat;
typedef animal::Matrix<double> Mat2;
typedef animal::Matrix<double> Mat3;
typedef animal::Matrix<double> Mat4;
typedef animal::Matrix<double> Mat5;


int main()
{
    using namespace animal;

    // vector operations
    {
    cout << endl << "==== Vector operations ====" << endl;
    Vec v,u;
    cout<<"enter a vector v in dimension "<< size(v) << ": "; 
    cin >> v_input(v);
    v_eq(u,v);   
    cout << "v= " << v_output(u) << endl;
    cout<<"enter a vector u in dimension "<< size(u) << ": "; 
    cin >> v_input(u);
    cout << "u= " << v_output(u) << endl;

    v_peq(v,u);      cout << "v+=u = " << v_output(v) << endl;
    v_meq(v,u);      cout << "v-=u = " << v_output(v) << endl;
    v_eq_ab( v,2,u );    cout << "v=2*u = " << v_output(v) << endl;
    v_eq( v,v,2,u );    cout << "v=v+2*u = " << v_output(v) << endl;
    v_peq( v,2,u );    cout << "v+=2*u = " << v_output(v) << endl;
    v_teq( v,2 );      cout << "v*=2 = " << v_output(v) << endl;
    v_assign( v, 1.1); cout << "v assigned to 1.1: " << v_output(v) << endl;
    }


    // matrix operations
    {
    cout << endl << "==== Matrix operations ====" << endl;
    {
    Mat N(nr,nc),M(nr,nc);
    cout << "enter a matrix M " << nrows(M) << "x" << ncols(M) << ": " << endl;
    cin >> m_input(M); cout << "M = " << m_output(M) << endl;
    m_eq_A(N,M);       cout << "N = M = " << m_output(N) << endl;
    m_teq(M,2 );       cout << "M*=2 = " << m_output(M) << endl;
    m_eq_aA(N,2,M);    cout << "N = 2M = " << m_output(N) << endl;
    m_peq_A(N,M);      cout << "N += M = " << m_output(N) << endl;
    m_meq_A(N,M);      cout << "N -= M = " << m_output(N) << endl;
    m_peq_aA(N,-1,M);  cout << "N += (-1)*M = " << m_output(N) << endl;
    m_assign( M, 0);   cout << "M assigned to 0.: " << m_output(M) << endl;
    Mat P; P.resize(3,3); P.fill(2); cout<<"uniform matrix in dimension 3: "<< P << endl;
    P.setIdentity(); cout<<"identity in dimension 3: "<< P << endl;
    }
    {
    Mat N(nr,nc); 
    Mat2 M(nr2,nc2);
    cout << "enter a matrix M " << nrows(M) << "x" << ncols(M) << ": " << endl;
    cin >> m_input(M); cout << "M = " << m_output(M) << endl;
    m_eq_At(N,M);       cout << "N = Mt = " << m_output(N) << endl;
    m_eq_aAt(N,2,M);    cout << "N = 2Mt = " << m_output(N) << endl;
    m_peq_At(N,M);      cout << "N += Mt = " << m_output(N) << endl;
    m_meq_At(N,M);      cout << "N -= Mt = " << m_output(N) << endl;
    m_peq_aAt(N,-1,M);  cout << "N += (-1)*Mt = " << m_output(N) << endl;
    }
    
    
    // matrix-vector operations
    Mat A(nr,nc);
    cout<<"enter a matrix A "<< nr <<"x"<< nc << ": "; 
    cin >> m_input(A); cout << "A = " << m_output(A) << endl;
    Vec v;
    cout << "enter vector v in dimension "<< size(v) <<": ";
    cin >> v_input(v); cout << "v = " << v_output(v) << endl;

    Vec2 w;
    v_eq_Ab(w,A,v);  cout <<"w = Av = "<< v_output(w) << endl;
    v_eq_aAb(w,2,A,v);  cout <<"w = 2Av = "<< v_output(w) << endl;
    v_peq_Ab(w,A,v);  cout <<"w += Av = "<< v_output(w) << endl;
    v_meq_Ab(w,A,v);  cout <<"w -= Av = "<< v_output(w) << endl;
    v_peq_aAb(w,2,A,v);  cout <<"w += 2Av = "<< v_output(w) << endl;
    
    Vec2 u;
    cout << "enter vector u in dimension "<< size(u) <<": ";
    cin >> v_input(u); cout << "u = " << v_output(u) << endl;
    
    v_eq_Atb(v,A,u);  cout <<"v = Atu = "<< v_output(v) << endl;
    v_eq_aAtb(v,2,A,u);  cout <<"v = 2Atu = "<< v_output(v) << endl;
    v_peq_Atb(v,A,u);  cout <<"v += Atu = "<< v_output(v) << endl;
    v_meq_Atb(v,A,u);  cout <<"v -= Atu = "<< v_output(v) << endl;
    v_peq_aAtb(v,2,A,u);  cout <<"v += 2Atu = "<< v_output(v) << endl;


    // matrix-matrix operations
    Mat2 B(nr2,nc2);
    cout<<"enter matrix B "<< nc <<"x"<< nc2 << ": "; 
    cin >> m_input(B); 
    cout << "B = " << m_output(B) << endl;
    
    Mat3 M(nr,nc2);
    m_eq_AB(M,A,B);  cout<<"M=AB = "<< m_output(M) << endl;
    m_eq_aAB(M,2,A,B);  cout<<"M=2AB = "<< m_output(M) << endl;
    m_peq_AB(M,A,B);  cout<<"M+=AB = "<< m_output(M) << endl;
    m_meq_AB(M,A,B);  cout<<"M-=AB = "<< m_output(M) << endl;
    m_peq_aAB(M,2,A,B);  cout<<"M+=2AB = "<< m_output(M) << endl;

    m_eq_ABt(M,A,A);  cout<<"M=AA^t = "<< m_output(M) << endl;
    m_eq_aABt(M,2.,A,A);  cout<<"M=2AA^t = "<< m_output(M) << endl;
    m_peq_ABt(M,A,A);  cout<<"M+=AA^t = "<< m_output(M) << endl;
    m_meq_ABt(M,A,A);  cout<<"M-=AA^t = "<< m_output(M) << endl;
    m_peq_aABt(M,2.,A,A);  cout<<"M+=2AA^t = "<< m_output(M) << endl;
    
    m_eq_AtB(M,B,B),   cout<<"M=B^tB = "<< m_output(M) << endl;
    m_eq_aAtB(M,2,B,B),   cout<<"M=2B^tB = "<< m_output(M) << endl;
    m_peq_AtB(M,B,B),   cout<<"M+=B^tB = "<< m_output(M) << endl;
    m_meq_AtB(M,B,B),   cout<<"M-=B^tB = "<< m_output(M) << endl;
    m_peq_aAtB(M,2.,B,B);  cout<<"M+=2B^tB = "<< m_output(M) << endl;

    Mat4 P(nc,nr2);
    m_eq_AtBt(P,A,B);  cout<<"P=AtBt = "<< m_output(P) << endl;
    m_eq_aAtBt(P,2,A,B);  cout<<"P=2AtBt = "<< m_output(P) << endl;
    m_peq_AtBt(P,A,B);  cout<<"P+=AtBt = "<< m_output(P) << endl;
    m_meq_AtBt(P,A,B);  cout<<"P-=AtBt = "<< m_output(P) << endl;
    m_peq_aAtBt(P,2,A,B);  cout<<"P+=2AtBt = "<< m_output(P) << endl;
    }
    
    
    // Cholesky solution of equation systems
    {
    cout << endl << "===== Solution of linear equation system Ax=b  using Cholesky decomposition ===" << endl;
    Mat5 A(nc,nc),L(nc,nc);  
    cout << "enter a symetric positive definite matrix A of size "<< nrows(A) <<"x"<< ncols(A) <<": ";
    cin >> m_input(A);
    cout << endl << "matrix A = " << m_output(A) << endl;
    m_assign(L,0.);
    m_eq_modchol(L,A);
    cout << "L = modified cholesky of A: " << m_output(L) << endl;
    Vec b,x,ver;
    cout << "enter a vector b of size " << size(b) << ": ";
    cin >> v_input(b); cout << "b= " << v_output(b) << endl;
    v_eq_cholsolve(x,L,b);
    cout << "solution: x= " << v_output(x) << endl;
    v_eq_Ab(ver,A,x); 
    cout << "check solution: Ax = " << v_output(ver) << endl;
    Mat5 T(nc,nc),Inv(nc,nc),Verif(nc,nc);
    m_assign(T,0);
    m_eq_invsqrt(T,A);
    cout << "T = invsqrt(A) = " << m_output(T) << endl;
    m_eq_AtB( Inv, T, T );
    m_eq_AB( Verif,Inv,A );
    cout << "check invsqrt(A): TtTA should be equal to identity: " << m_output(Verif) << endl;
    m_eq_cholinv(T,A,L);
    cout << "T=inverse of A: " << m_output(T) << endl;
    m_eq_AB(L,T,A);
    cout << "check inverse: TA should be equal to identity: " << m_output(L) << endl;
    
    
//  Mat5 verif;  
//  m_eq_ABt(verif,L,L);
//  cout << "verif: LL^T= " << m_output(verif) << endl;
    }
     
    // LU decomposition, solution, inversion
    {
    cout << endl << "===== Solution of linear equation system Ax=b using LU decomposition ===" << endl;
    Mat5 A(nc,nc), LU(nc,nc);
    cout << "enter a square matrix A of size "<< nrows(A) <<"x"<< ncols(A) <<": ";
    cin >> m_input(A); cout << "A = " << m_output(A) << endl;
    // LU = LU decomposition in A
    m_eq_ludcmp( LU,A );
    cout << "LU decomposition of A: " << m_output(LU) << endl;
    
    // LU solution
    Vec b,x,verif;
    cout << "enter a vector of size " << size(b) << endl;
    cin >> v_input(b); cout << "b = " << v_output(b) << endl;
    v_eq_lusolve(x,LU,b);
    cout << "solution of Ax=b: " << v_output(x) << endl;
    v_eq_Ab(verif,A,x);
    cout << "verif: Ax = " << v_output(verif) << endl;
    
    // LU inversion
    Mat5 Inv(nc,nc), VerifInv(nc,nc);
    Vec aux;
    m_eq_luinv( Inv, LU, aux );
    cout << "inverse of A: " << m_output(Inv) << endl;
    m_eq_AB( VerifInv, A, Inv );
    cout << "verif: A*inv(A) = " << m_output(VerifInv) << endl;
    m_eq_AB( VerifInv, Inv, A );
    cout << "verif: inv(A)*A = " << m_output(VerifInv) << endl;
    }
    
    
    
return 0;
}

---