tjacob.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 <cassert>

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

using namespace std;
using namespace LinAlg;

void FillEigenProjs(Matrix<REAL> const & Q, int N, Matrix<REAL> P[])
{
    assert(Q.Rows()==N);
    assert(Q.Cols()==N);
    for (int k=0; k<N; ++k)
        for (int i=0; i<N; ++i)
        for (int j=0; j<N; ++j)
            P[k](i,j) = Q(i,k) * Q(j,k);
}

int main()
{
    //Matrix<REAL> A(3,3,"2 0 1   0 3 0   1 0 4");
    Matrix<REAL> A(3,3,"4 5 6   5 8 7   6 7 9");
    Matrix<REAL> Q(3,3);
    Vector<REAL> v(3); // Eigenvalues
    Matrix<REAL> TMP(3,3);
    Matrix<REAL> * P = new Matrix<REAL> [3]; // Eigenprojectors
    for (int i=0; i<3; ++i) P[i].Set(3,3);

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

    cout << "\n=========================================== JacobiRot Function\n";

    Matrix<REAL> AA(3,3); AA=A;
    int itrs = Sym::JacobiRot(A,Q,v);
    cout << "Iterations = " << itrs << endl;
    cout << "A = \n"        << A    << endl;
    cout << "Q = \n"        << Q    << endl;
    cout << "v = \n"        << v    << endl;

    FillEigenProjs(Q,3,P);
    TMP.SetValues(0.0);
    for (int i=0; i<3; ++i)
    for (int j=0; j<3; ++j)
        TMP(i,j) = TMP(i,j) + v(0)*P[0](i,j)
                            + v(1)*P[1](i,j)
                            + v(2)*P[2](i,j);
    cout << "Sum{eigenvals * eigenprojs} = \n" << TMP << endl;

#ifdef USE_LAPACK
    cout << "\n=========================================== using LAPACK\n";

    A=AA;
    Syevr(A,v,Q);
    cout << "A = \n" << A << endl;
    cout << "v = \n" << v << endl;
    cout << "Q = \n" << Q << endl;

    FillEigenProjs(Q,3,P);
    TMP.SetValues(0.0);
    for (int i=0; i<3; ++i)
    for (int j=0; j<3; ++j)
        TMP(i,j) = TMP(i,j) + v(0)*P[0](i,j)
                            + v(1)*P[1](i,j)
                            + v(2)*P[2](i,j);
    cout << "Sum{eigenvals * eigenprojs} = \n" << TMP << endl;
#endif

    delete [] P;
    return 0;
}

---