teigen1.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                                           *
 *************************************************************************************/

#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()
{
    Matrix<REAL> A(4,4,"5 2 3 1   2 8 1 1   3 1 8 5   1 1 5 10");
    Vector<REAL> v(4);    // Eigenvalues
    Matrix<REAL> * P = new Matrix<REAL> [4];
    for (int i=0; i<4; ++i) P[i].Set(4,4); // Eigenprojectors

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

    Syep(A,v,P);
    for (int i=0; i<4; ++i)
        cout << "P[" << i << "] =\n" << P[i] << endl;

    Matrix<REAL> Check(4,4);
    for (int i=0; i<4; ++i)
    for (int j=0; j<4; ++j)
        for (int k=0; k<4; ++k)
            Check(i,j) = Check(i,j) + v(k)*P[k](i,j);
    cout << "v0*P0 + v1*P1 + v2*P2 + v3*P3 =\n" << Check << endl;

    for (int i=0; i<4; ++i)
    {
        double trace=0.0;
        for (int j=0; j<4; ++j)
            trace += P[i](j,j);
        cout << "trace(P[" << i << "]) = " << trace << endl;
    }
    
    delete [] P;

    return 0;
}

---