tmv.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/lawrap.h"
#include "linalg/matrix.h"
#include "linalg/vector.h"
#include "util/exception.h"

using namespace std;
using namespace LinAlg;

int main() try
{
    // =============================================================== Matrix x Vector

    cout << "\n============================================== Matrix x Vector\n";
    
    Matrix<REAL> A(5,5,"3 1 3 1 8   1 5 9 -2 4   2 6 5 0 -7   1 2 2 3 1   1 2 3 1 -1");
    Vector<REAL> X(5,"-1 -1  1 0 2");
    Vector<REAL> Y(5); Y.SetValues(1);

    cout << "A = \n" << A << endl;
    cout << "X = \n" << X << endl;
    cout << "Y = \n" << Y << endl;
    
    double a =  1.1;
    double b = -8.0;
    Gemv(a,A,X,b,Y);
    cout << "Y := " << a << "*A*X + " << b << "*Y = \n" << Y << endl;
    
    // =============================================================== Matrix x Matrix
    
    cout << "\n============================================== Matrix x Matrix\n";

    Matrix<REAL> B(5,5,"4 1 2 1 1   1 5 9 0 3   2 4 2 -1 2   0 0 1 -1 1   2 0 0 -1 1");
    Matrix<REAL> C(5,5,"0 1 -1 0 1   1 2 3 1 1  -1 0 1 1 1   2 -1 0 1 1  -2 -1 -1 -1 0");
    Matrix<REAL> E(5,5); E.SetValues(1);

    cout << "B = \n" << B << endl;
    cout << "C = \n" << C << endl;
    cout << "E = \n" << E << endl;

    a = -2.0;
    b =  3.0;
    Gemm(a,B,C,b,E);
    cout << "E := " << a << "*B*C + " << b << "*E = \n" << E << endl;

    // =========================================================== Matrix^(T) x Matrix
    
    cout << "\n========================================== Matrix^(T) x Matrix\n";

    E.SetValues(1);

    Gemtm(a,B,C,b,E);
    cout << "E := " << a << "*(B^T)*C + " << b << "*E = \n" << E << endl;

    // =========================================================== Matrix x Matrix^(T)
    
    cout << "\n========================================== Matrix x Matrix^(T)\n";

    E.SetValues(1);

    Gemmt(a,B,C,b,E);
    cout << "E := " << a << "*B*(C^T) + " << b << "*E = \n" << E << endl;

    // ======================================================= Matrix^(T) x Matrix^(T)
    
    cout << "\n====================================== Matrix^(T) x Matrix^(T)\n";

    E.SetValues(1);

    Gemtmt(a,B,C,b,E);
    cout << "E := " << a << "*(B^T)*(C^T) + " << b << "*E = \n" << E << endl;

    // ======================================================= Vector x Vector^(T)
    
    cout << "\n========================================== Vector x Vector^(T)\n";

    Matrix<REAL> F(5,9); F.SetValues(0);
    Vector<REAL> Z(9,"4 3 2 1 0 -1 2 3 4");
    X.Reset("0 1 2 3 4");

    cout << "F = \n" << F << endl;
    cout << "X = \n" << X << endl;
    cout << "Z = \n" << Z << endl;
    
    Ger(1,X,Z,F);
    cout << "F := F + " << 1 << "*X*Z^T = \n" << F << endl;

    return 0;
}
catch (Exception * e) //{{{ 
{
    e->Cout();
    if (e->IsFatal()) {delete e; exit(1);}
    delete e;
}
catch (char const * m)
{
    std::cout << "Fatal: " << m << std::endl;
    exit (1);
}
catch (...)
{
    std::cout << "Some exception (...) ocurred\n";
} //}}} 

// vim:fdm=marker

---