teigen.cpp

/*************************************************************************************
 * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
 * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
 * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
 *                                                                                   *
 * This file is part of MechSys.                                                     *
 *                                                                                   *
 * MechSys is free software; you can redistribute it and/or modify it under the      *
 * terms of the GNU General Public License as published by the Free Software         *
 * Foundation; either version 2 of the License, or (at your option) any later        *
 * version.                                                                          *
 *                                                                                   *
 * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
 * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
 *                                                                                   *
 * You should have received a copy of the GNU General Public License along with      *
 * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
 * Fifth Floor, Boston, MA 02110-1301, USA                                           *
 *************************************************************************************/

#define _NMBFMT std::setw(12) <<

#include <iostream>
#include <cmath>
#include <iomanip>

#include "linalg/matrix.h"
#include "linalg/vector.h"
#include "linalg/lawrap.h"

using namespace std;
using namespace LinAlg;

int main()
{

    // ================================================================== Eigenproblem

    cout << "\n=========================================== EigenValues/Vectors\n";
    
    Matrix<REAL> P(5,5,"10 2 3 4 5  2 20 3 4 5  3 3 30 4 5  4 4 4 40 5  5 5 5 5 50");
    Vector<REAL> v(5); // Eigenvalues
    Matrix<REAL> Q(5,5); // Eigenvectors (columns of Q)

    Matrix<REAL> R(5,5); R = P;
    Syevr(R,v,Q);

    cout << "P = \n" << P << endl;
    cout << "R = \n" << R << endl;
    cout << "v = \n" << v << endl;
    cout << "Q = \n" << Q << endl;

    Vector<REAL> * V = new Vector<REAL> [5];
    for (int i=0; i<5; ++i) V[i].Set(5);
    Vector<REAL> S(5);

    for (int j=0; j<5; ++j)
        for (int i=0; i<5; ++i)
            V[j](i) = Q(i,j);

    for (int i=0; i<5; ++i)
    {
        Gemv(1,P,V[i],0,S);
        for (int j=0; j<5; ++j)
            S(j) = S(j) - v(i)*V[i](j);
        cout << "S = P*V" << i << " - v" << i << "*V" << i << "\n";
        cout << S << endl;
    }

    delete [] V;
    
    cout << "\n============================================== EigenProjectors\n";

    Matrix<REAL> F(4,4,"5 2 3 1   2 8 1 1   3 1 8 5   1 1 5 10");
    Vector<REAL> G(4);    // Eigenvalues
    Matrix<REAL> * H = new Matrix<REAL> [4];
    for (int i=0; i<4; ++i) H[i].Set(4,4); // Eigenprojectors

    cout << "F = \n" << F << endl;

    Syep(F,G,H);
    for (int i=0; i<4; ++i)
        cout << "H[" << i << "] =\n" << H[i] << endl;

    Matrix<REAL> I(4,4);
    for (int i=0; i<4; ++i)
    for (int j=0; j<4; ++j)
        for (int k=0; k<4; ++k)
            I(i,j) = I(i,j) + G(k)*H[k](i,j);
    cout << "I = v0*H0 + v1*H1 + v2*H2 + v3*H3 =\n" << I << endl;
    
    Matrix<REAL> T(4,4);
    for (int i=0; i<4; ++i)
    for (int j=0; j<4; ++j)
    {
        if (i==j)
        {
            T = H[i];
            Gemm(1,H[i],H[j],-1,T);
            cout << "H" << i << "*H" << j << " - H" << i << " =\n";
            cout << T << endl;
        }
        else
        {
            Gemm(1,H[i],H[j],0,T);
            cout << "H" << i << "*H" << j << " =\n" << T << endl;
        }
    }

    delete [] H;
    
    cout << "\n=========================================== Isotropic Functions\n";

    Matrix<REAL> IF(4,4,"5 2 3 1   2 8 1 1   3 1 8 5   1 1 5 10");

    Syif(F,IF,exp);
    cout << "exp(F) = \n" << IF << endl;

    F.Reset("5 2 3 1   2 8 1 1   3 1 8 5   1 1 5 10");
    cout << "F = \n" << F << endl;

    Syif(F,IF,sqrt);
    cout << "sqrt(F) = \n" << IF << endl;
    
    // Check
    Matrix<REAL> RR(4,4);
    Gemm(1,IF,IF,0,RR);
    cout << "sqrt(F)*sqrt(F) = \n" << RR << endl;
    
    cout << "\n=========================================== Isotropic Functions\n";

    int n = 7;
    Matrix<REAL> AA(7,7);
    for (int i=0; i<n-1; ++i) AA(i+1,i) = -1;
    for (int i=0; i<n-1; ++i) AA(i,i+1) = -1;
    for (int i=1; i<n-1; ++i) AA(i,i)   = 2;
    AA(0,0) = 1;
    AA(n-1,n-1) = 1;

    cout << "AA = \n" << AA << endl;
    
    Matrix<REAL> sqAA(7,7);
    Syif(AA,sqAA,sqrt);

    cout << "sqrt(AA) = \n" << sqAA << endl;

    // Check
    Matrix<REAL> RRR(7,7);
    Gemm(1,sqAA,sqAA,0,RRR);
    cout << "sqrt(AA)*sqrt(AA) = \n" << RRR << endl;

    return 0;
}

---