matvectest.cc File Reference

#include "mbconfig.h"
#include <cassert>
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <typeinfo>
#include "ac/f2c.h"
#include "clock_time.h"
#include "matvec.h"
#include "matvec3.h"
#include "Rot.hh"
#include "matvecass.h"

Include dependency graph for matvectest.cc:

Go to the source code of this file.

Classes
struct	testMatVecProductGradient_testData::test_b
struct	testMatVecProductGradient_testData::test_c
struct	testMatVecProductGradient_testData::test_d
struct	testMatVecProductGradient_testData::test_e
struct	testMatVecProductGradient_testData::test_f
struct	testMatVecProductGradient_testData::test_g
struct	testMatVecProductGradient_testData::test_i
struct	testMatVecProductGradient_testData::test_j
struct	testMatVecProductGradient_testData::test_l
struct	testMatVecProductGradient_testData::test_m
struct	testMatVecProductGradient_testData::test_r
struct	testMatVecProductGradient_testData::test_s
struct	testMatVecProductGradient_testData::test_t
struct	testMatVecProductGradient_testData::test_norm_g

Namespaces
	testMatVecProductGradient_testData
	gradVecAssTest

Macros
#define	MATVEC_DEBUG   1
#define	GRADIENT_DEBUG   1

Functions
void	tic ()
void	tic (doublereal &dTime)
doublereal	toc ()
doublereal	random1 ()
void	testScalarTypeTraits ()
template<typename T , index_type N>
void	func (const Vector< T, N > &a, const Vector< T, N > &b, Vector< T, N > &c)
template<typename T >
void	func (const T a[], const T b[], T c[], index_type N)
template<typename T , index_type N>
void	func2 (const Matrix< T, N, N > &A, const Vector< T, N > &b, const Vector< T, N > &c, Vector< T, N > &d, doublereal e, doublereal &dt)
template<typename T >
void	func2 (const T *A, const T *b, const T *c__, T *d__, const doublereal &e, doublereal &dt)
void	func2ad (const doublereal A[3][3], const doublereal b[3], const doublereal c[3], doublereal d[3], const doublereal &e)
void	func2ad_dv (const doublereal A[3][3], const doublereal Ad[3][3][nbdirsmax], const doublereal b[3], const doublereal bd[3][nbdirsmax], const doublereal c[3], const doublereal cd[3][nbdirsmax], doublereal d[3], doublereal dd[3][nbdirsmax], const doublereal &e, const integer &nbdirs)
void	func2ad (const Matrix< doublereal, 3, 3 > &A, const Vector< doublereal, 3 > &b, const Vector< doublereal, 3 > &c, Vector< doublereal, 3 > &d, const doublereal &e, LocalDofMap *, doublereal &dt)
template<index_type N_SIZE>
void	func2ad (const Matrix< Gradient< N_SIZE >, 3, 3 > &A, const Vector< Gradient< N_SIZE >, 3 > &b, const Vector< Gradient< N_SIZE >, 3 > &c, Vector< Gradient< N_SIZE >, 3 > &d, const doublereal &e, LocalDofMap *pDofMap, doublereal &dt)
template<typename T >
void	func3 (const Matrix< T, 3, 3 > &R1, const Matrix< T, 3, 3 > &R2, const Vector< T, 3 > &a, const Vector< T, 3 > &b, Vector< T, 3 > &c, doublereal e)
template void	func3< doublereal > (const Matrix< doublereal, 3, 3 > &R1, const Matrix< doublereal, 3, 3 > &R2, const Vector< doublereal, 3 > &a, const Vector< doublereal, 3 > &b, Vector< doublereal, 3 > &c, doublereal e)
template void	func3< Gradient< 0 > > (const Matrix< Gradient< 0 >, 3, 3 > &R1, const Matrix< Gradient< 0 >, 3, 3 > &R2, const Vector< Gradient< 0 >, 3 > &a, const Vector< Gradient< 0 >, 3 > &b, Vector< Gradient< 0 >, 3 > &c, doublereal e)
template<typename T >
bool	bCompare (const T &a, const T &b, doublereal dTolRel=0.)
template<index_type N_SIZE>
bool	bCompare (const Gradient< N_SIZE > &a, const Gradient< N_SIZE > &b, doublereal dTolRel=0.)
template<typename T , index_type N>
void	callFunc (LocalDofMap *pDofMap, const Vector< T, N > &a, const Vector< T, N > &b, Vector< T, N > &c, Vector< T, N > &c1)
template<typename T >
void	callFunc2 (LocalDofMap *pDofMap, const Matrix< T, 3, 3 > &A, const Vector< T, 3 > &b, const Vector< T, 3 > &c, Vector< T, 3 > &d, Vector< T, 3 > &d_C, Vector< T, 3 > &d_F)
template<index_type N>
void	testMatVecGradient (doublereal c_C[N], doublereal cd_C[N][N])
template<index_type N>
void	testMatVecDouble (doublereal c_C[N])
template<index_type N_SIZE>
void	testMatVecGradient2 ()
void	testMatVecDouble2 ()
void	testMatVecProduct ()
template<typename S , index_type N_rows, index_type N_SIZE>
void	testMatVecProductGradient_testData::testGradient (const S &ref, const Vector< Gradient< N_SIZE >, N_rows > &v, doublereal dTol)
template<typename S , index_type N_SIZE>
void	testMatVecProductGradient_testData::testGradient (const S &ref, const Gradient< N_SIZE > &g, doublereal dTol)
template<typename T , index_type N_rows>
T	Norm_1 (const Vector< T, N_rows > &u)
void	func2addad_dv (const doublereal x[], const doublereal xd[], const doublereal y[], const doublereal yd[], doublereal z[], doublereal zd[], const integer &n, const integer &nbdirs)
void	func2mulad_dv (const doublereal x[], const doublereal xd[], const doublereal y[], const doublereal yd[], doublereal z[], doublereal zd[], const integer &n, const integer &nbdirs)
template<typename T , index_type iRowCount>
void	doVecAdd (const Vector< T, iRowCount > &x, const Vector< T, iRowCount > &y, Vector< T, iRowCount > &z)
template<typename T , index_type iRowCount>
void	doVecMul (const Vector< T, iRowCount > &x, const Vector< T, iRowCount > &y, Vector< T, iRowCount > &z)
template<index_type iRowCount, index_type iMaxDeriv, typename Function , typename FunctionF77 >
void	testVecOp (const int M, const int N, Function f, FunctionF77 f77, const char *function)
void	testMatVecProductGradient ()
template<index_type N_SIZE>
void	testMatVecProductGradient2 (index_type iNumDeriv, int N)
template<index_type N_rows>
void	testMatVecCopy ()
template<typename T >
Matrix< T, 3, 3 > &	gradVecAssTest::Euler123ToMatR (const Vector< T, 3 > &v, Matrix< T, 3, 3 > &R)
void	testMatVec3 ()
void	testSubVecAss ()
void	testSubVecAssMatVec ()
void	testInv ()
void	testSolve ()
template<index_type N>
void	testMatVec ()
template<int iNumDofMax>
void	cppad_benchmark1 (const int N)
template<int iNumDofMax>
void	cppad_benchmark2 (const int N)
template<int iNumDofMax>
void	cppad_benchmark3 (const int N)
void	Mat3xN_test (int N, int M)
void	MatNxN_test (int N, int M)
void	MatDynamic_test (index_type iNumRows, index_type iNumCols, index_type iNumLoops)
template<index_type N_SIZE>
void	Mat3xN_test_grad (int iNumDeriv, int N, int M)
template<index_type N_SIZE>
void	Mat3xNT_test_grad (int iNumDeriv, int N, int M)
template<index_type N_SIZE>
void	MatNxNT_test_grad (int iNumDeriv, int N, int M)
template<index_type N_SIZE>
void	MatNxN_test_grad (int iNumDeriv, int N, int M)
template<index_type N_SIZE>
void	MatDynamic_test_grad (index_type iNumDeriv, index_type iNumRows, index_type iNumCols, int iNumLoops)
template<index_type N_ROWS, index_type N_COLS>
void	MatDynamicT_test (index_type iNumRows, index_type iNumCols, int iNumLoops)
template<index_type N_DERIV, index_type N_ROWS, index_type N_COLS>
void	MatDynamicT_test_grad (index_type iNumDeriv, index_type iNumRows, index_type iNumCols, int iNumLoops)
void	MatManip_test (int NLoops)
int	main (int argc, char *argv[])

Variables
int	NLoops = 1
int	NLoopsAss = 1
const integer	nbdirsmax = 12
static const struct testMatVecProductGradient_testData::test_b	testMatVecProductGradient_testData::oct_b
static const struct testMatVecProductGradient_testData::test_c	testMatVecProductGradient_testData::oct_c
static const struct testMatVecProductGradient_testData::test_d	testMatVecProductGradient_testData::oct_d
static const struct testMatVecProductGradient_testData::test_e	testMatVecProductGradient_testData::oct_e
static const struct testMatVecProductGradient_testData::test_f	testMatVecProductGradient_testData::oct_f
static const struct testMatVecProductGradient_testData::test_g	testMatVecProductGradient_testData::oct_g
static const struct testMatVecProductGradient_testData::test_i	testMatVecProductGradient_testData::oct_i
static const struct testMatVecProductGradient_testData::test_j	testMatVecProductGradient_testData::oct_j
static const struct testMatVecProductGradient_testData::test_l	testMatVecProductGradient_testData::oct_l
static const struct testMatVecProductGradient_testData::test_m	testMatVecProductGradient_testData::oct_m
static const struct testMatVecProductGradient_testData::test_r	testMatVecProductGradient_testData::oct_r
static const struct testMatVecProductGradient_testData::test_s	testMatVecProductGradient_testData::oct_s
static const struct testMatVecProductGradient_testData::test_t	testMatVecProductGradient_testData::oct_t
static const struct testMatVecProductGradient_testData::test_norm_g	testMatVecProductGradient_testData::oct_norm_g
const int	gradVecAssTest::I1 = 1
const int	gradVecAssTest::I2 = 2
const int	gradVecAssTest::I3 = 3
doublereal	dStartTime

Macro Definition Documentation

#define GRADIENT_DEBUG   1

Definition at line 70 of file matvectest.cc.

Referenced by main().

#define MATVEC_DEBUG   1

Definition at line 66 of file matvectest.cc.

Referenced by main().

Function Documentation

template<typename T >

bool bCompare	(	const T &	a,
		const T &	b,
		doublereal	dTolRel = 0.
	)

Definition at line 367 of file matvectest.cc.

Referenced by callFunc(), callFunc2(), Mat3xN_test(), Mat3xN_test_grad(), Mat3xNT_test_grad(), MatDynamic_test(), MatDynamic_test_grad(), MatDynamicT_test(), MatDynamicT_test_grad(), MatManip_test(), MatNxN_test(), MatNxN_test_grad(), MatNxNT_test_grad(), testMatVecProductGradient_testData::testGradient(), testInv(), testMatVec(), testMatVecDouble(), testMatVecGradient(), testMatVecGradient2(), testMatVecProduct(), and testMatVecProductGradient().

                                                               {
    assert(!std::isnan(a));
    assert(!std::isnan(b));
        doublereal dTolAbs = std::max<T>(1., std::max<T>(std::abs(a), std::abs(b))) * dTolRel;
        return std::abs(a - b) <= dTolAbs;
}

template<index_type N_SIZE>

bool bCompare	(	const Gradient< N_SIZE > &	a,
		const Gradient< N_SIZE > &	b,
		doublereal	dTolRel = 0.
	)

Definition at line 375 of file matvectest.cc.

References grad::Gradient< N_SIZE >::dGetDerivativeLocal(), grad::Gradient< N_SIZE >::dGetValue(), grad::Gradient< N_SIZE >::iGetEndIndexLocal(), and grad::Gradient< N_SIZE >::iGetStartIndexLocal().

                                                                                             {
    assert(!std::isnan(a.dGetValue()));
    assert(!std::isnan(b.dGetValue()));
        doublereal dTolAbs = std::max(1., std::max(std::abs(a.dGetValue()), std::abs(b.dGetValue()))) * dTolRel;

        if (std::abs(a.dGetValue() - b.dGetValue()) > dTolAbs) {
                std::cerr << "a = " << a.dGetValue() << " != b = " << b.dGetValue() << std::endl;
                return false;
        }

        index_type iStart = std::min(a.iGetStartIndexLocal(), b.iGetStartIndexLocal());
        index_type iEnd = std::max(a.iGetEndIndexLocal(), b.iGetEndIndexLocal());

        for (index_type i = iStart; i < iEnd; ++i) {
                doublereal dTolAbs = std::max<scalar_deriv_type>(1., std::max<scalar_deriv_type>(std::abs(a.dGetDerivativeLocal(i)), std::abs(b.dGetDerivativeLocal(i)))) * dTolRel;
        assert(!std::isnan(a.dGetDerivativeLocal(i)));
        assert(!std::isnan(b.dGetDerivativeLocal(i)));
                if (std::abs(a.dGetDerivativeLocal(i) - b.dGetDerivativeLocal(i)) > dTolAbs) {
                        std::cerr << "ad(" << i << ") = " << a.dGetDerivativeLocal(i) << " != bd(" << i << ") = " << b.dGetDerivativeLocal(i) << std::endl;
                        return false;
                }
        }

        return true;
}

Here is the call graph for this function:

template<typename T , index_type N>

void callFunc	(	LocalDofMap *	pDofMap,
		const Vector< T, N > &	a,
		const Vector< T, N > &	b,
		Vector< T, N > &	c,
		Vector< T, N > &	c1
	)

Definition at line 402 of file matvectest.cc.

References bCompare(), c, func(), NLoops, grad::Vector< T, N_rows >::pGetVec(), grad::sqrt(), tic(), and toc().

Referenced by testMatVecDouble(), and testMatVecGradient().

                                                                                                                     {
        srand(0);
    tic();
    for (int i = 0; i < NLoops; ++i) {
        func(a, b, c);
    }

    std::cerr << "Vector (" << typeid(T).name() << "): " << toc() << "s" << std::endl;

    srand(0);
    tic();
    for (int i = 0; i < NLoops; ++i) {
        func(a.pGetVec(), b.pGetVec(), c1.pGetVec(), N);
    }

    std::cerr << "Array (" << typeid(T).name() << "): " << toc() << "s" << std::endl;

    std::cout << "a=" << std::endl << a << std::endl;
    std::cout << "b=" << std::endl << b << std::endl;
    std::cout << "c=" << std::endl << c << std::endl;
    std::cout << "c1=" << std::endl << c1 << std::endl;

    const doublereal dTol = sqrt(std::numeric_limits<doublereal>::epsilon());

    for (int i = 1; i <= N; ++i) {
        assert(bCompare(c(i), c1(i), dTol));
    }
}

Here is the call graph for this function:

template<typename T >

void callFunc2	(	LocalDofMap *	pDofMap,
		const Matrix< T, 3, 3 > &	A,
		const Vector< T, 3 > &	b,
		const Vector< T, 3 > &	c,
		Vector< T, 3 > &	d,
		Vector< T, 3 > &	d_C,
		Vector< T, 3 > &	d_F
	)

Definition at line 432 of file matvectest.cc.

References bCompare(), func2(), func2ad(), NLoops, grad::Matrix< T, N_rows, N_cols >::pGetMat(), grad::Vector< T, N_rows >::pGetVec(), and random1().

Referenced by testMatVecDouble2(), and testMatVecGradient2().

                                                                                                                                                                    {
        srand(0);
    doublereal dtMatVec = 0., dtC = 0., dtF = 0.;
    for (int i = 0; i < NLoops; ++i) {
        func2(A, b, c, d, random1(), dtMatVec);
    }

    std::cerr << "matvec (" << typeid(T).name() << "): " << dtMatVec << "s" << std::endl;

    srand(0);
    for (int i = 0; i < NLoops; ++i) {
        func2(A.pGetMat(), b.pGetVec(), c.pGetVec(), d_C.pGetVec(), random1(), dtC);
    }

    std::cerr << "C (" << typeid(T).name() << "): " << dtC << "s" << std::endl;

    srand(0);
    for (int i = 0; i < NLoops; ++i) {
        func2ad(A, b, c, d_F, random1(), pDofMap, dtF);
    }

    std::cerr << "F77 (" << typeid(T).name() << "): " << dtF << "s" << std::endl;

    std::cerr << "overhead matvec:" << dtMatVec / std::max(std::numeric_limits<double>::epsilon(), dtF) << std::endl;
    std::cerr << "overhead C:" << dtC / std::max(std::numeric_limits<double>::epsilon(), dtF) << std::endl;

    std::cout << "A=" << std::endl << A << std::endl;
    std::cout << "b=" << std::endl << b << std::endl;
    std::cout << "c=" << std::endl << c << std::endl;
    std::cout << "d=" << std::endl << d << std::endl;
    std::cout << "d_C=" << std::endl << d_C << std::endl;
    std::cout << "d_F=" << std::endl << d_F << std::endl;

    const doublereal dTol = 10 * std::numeric_limits<scalar_deriv_type>::epsilon();

    for (int i = 1; i <= 3; ++i) {
        assert(bCompare(d(i), d_C(i), dTol));
        assert(bCompare(d(i), d_F(i), dTol));
    }
}

Here is the call graph for this function:

template<int iNumDofMax>

void cppad_benchmark1

(

const int

N

)

Definition at line 2919 of file matvectest.cc.

References dof, gradVecAssTest::Euler123ToMatR(), grad::MapVectorBase::GLOBAL, mbdyn_clock_time(), and R.

                                   {
    const int iNumDof = 6;
    
    assert(iNumDofMax == 0 || iNumDofMax >= iNumDof);
    
    Matrix<Gradient<iNumDofMax>, 3, 3> R;
    Vector<Gradient<iNumDofMax>, 3> X, Y, Phi;
    Vector<doublereal, 3> o;
    LocalDofMap dof;
        
    o(1) = 1.;
    o(2) = 2.;
    o(3) = 3.;
    
    Matrix<doublereal, 3, iNumDof> jac;
    
    const double start = mbdyn_clock_time();
    double calc = 0.;
    
    for (int loop = 0; loop < N; ++loop) {            
        for (int i = 0; i < 3; ++i) {
            X(i + 1).SetValuePreserve(0.);
            X(i + 1).DerivativeResizeReset(&dof, i, MapVectorBase::GLOBAL, 1.);
            Phi(i + 1).SetValuePreserve(0.);
            Phi(i + 1).DerivativeResizeReset(&dof, i + 3, MapVectorBase::GLOBAL, 1.);
        }
        
        const double start_calc = mbdyn_clock_time();
        
        gradVecAssTest::Euler123ToMatR(Phi, R);
        
        Y = X + R * o;
        
        calc += mbdyn_clock_time() - start_calc;
        
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < iNumDof; ++j) {
                jac(i + 1, j + 1) = Y(i + 1).dGetDerivativeGlobal(j);
            }
        }
    }
    
    std::cout << "calculation time: " << calc/N << "s\n";
    std::cout << "elapsed time: " << (mbdyn_clock_time() - start)/N << "s\n";
    
    std::cout << "x=" << X << std::endl;
    std::cout << "y=" << Y << std::endl;
    std::cout << "jac=\n";

    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < iNumDof; ++j) {
            std::cout << jac(i + 1, j + 1) << '\t';
        }
        std::cout << std::endl;
    }
}

Here is the call graph for this function:

template<int iNumDofMax>

void cppad_benchmark2

(

const int

N

)

Definition at line 2977 of file matvectest.cc.

References dof, gradVecAssTest::Euler123ToMatR(), grad::MapVectorBase::GLOBAL, and mbdyn_clock_time().

                                   {
    const int iNumDof = 12;
    
    assert(iNumDofMax == 0 || iNumDofMax >= iNumDof);
    
    Matrix<Gradient<iNumDofMax>, 3, 3> R1, R2;
    Vector<Gradient<iNumDofMax>, 3> X1, X2, Y, Phi1, Phi2;
    Vector<doublereal, 3> o1, o2;
    LocalDofMap dof;
        
    o1(1) = 1.;
    o1(2) = 2.;
    o1(3) = 3.;

    o2(1) = 1.;
    o2(2) = 2.;
    o2(3) = 3.;
    
    Matrix<doublereal, 3, iNumDof> jac;
    
    const double start = mbdyn_clock_time();
    double calc = 0.;
    
    for (int loop = 0; loop < N; ++loop) {            
        for (int i = 0; i < 3; ++i) {
            X1(i + 1).SetValuePreserve(0.);
            X1(i + 1).DerivativeResizeReset(&dof, i, MapVectorBase::GLOBAL, 1.);
            Phi1(i + 1).SetValuePreserve(0.);
            Phi1(i + 1).DerivativeResizeReset(&dof, i + 3, MapVectorBase::GLOBAL, 1.);
            X2(i + 1).SetValuePreserve(0.);
            X2(i + 1).DerivativeResizeReset(&dof, i + 6, MapVectorBase::GLOBAL, 1.);
            Phi2(i + 1).SetValuePreserve(0.);
            Phi2(i + 1).DerivativeResizeReset(&dof, i + 9, MapVectorBase::GLOBAL, 1.);            
        }
        
        const double start_calc = mbdyn_clock_time();
        
        gradVecAssTest::Euler123ToMatR(Phi1, R1);
        gradVecAssTest::Euler123ToMatR(Phi2, R2);
        
        Y = X1 + R1 * o1 - X2 - R2 * o2;
        
        calc += mbdyn_clock_time() - start_calc;
        
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < iNumDof; ++j) {
                jac(i + 1, j + 1) = Y(i + 1).dGetDerivativeGlobal(j);
            }
        }
    }
    
    std::cout << "calculation time: " << calc/N << "s\n";
    std::cout << "elapsed time: " << (mbdyn_clock_time() - start)/N << "s\n";
    
    std::cout << "X1=" << X1 << std::endl;
    std::cout << "X2=" << X2 << std::endl;
    std::cout << "Y=" << Y << std::endl;
    std::cout << "jac=\n";

    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < iNumDof; ++j) {
            std::cout << jac(i + 1, j + 1) << '\t';
        }
        std::cout << std::endl;
    }
}

Here is the call graph for this function:

template<int iNumDofMax>

void cppad_benchmark3

(

const int

N

)

Definition at line 3045 of file matvectest.cc.

References dof, gradVecAssTest::Euler123ToMatR(), grad::MapVectorBase::GLOBAL, mbdyn_clock_time(), and grad::Transpose().

                                   {
    const int iNumDof = 12;
    
    assert(iNumDofMax == 0 || iNumDofMax >= iNumDof);
    
    typedef Matrix<Gradient<iNumDofMax>, 3, 3> Mat3x3;
    typedef Vector<Gradient<iNumDofMax>, 3> Vec3;
    typedef Vector<doublereal, 3> CVec3;
    
    Mat3x3 R1, R2;
    Vec3 X1, X2, Y, Phi1, Phi2;
    CVec3 o1, o2;
    LocalDofMap dof;
        
    o1(1) = 1.;
    o1(2) = 2.;
    o1(3) = 3.;

    o2(1) = 1.;
    o2(2) = 2.;
    o2(3) = 3.;
    
    Matrix<doublereal, 3, iNumDof> jac;
    
    const double start = mbdyn_clock_time();
    double calc = 0.;
    
    for (int loop = 0; loop < N; ++loop) {            
        for (int i = 0; i < 3; ++i) {
            X1(i + 1).SetValuePreserve(0.);
            X1(i + 1).DerivativeResizeReset(&dof, i, MapVectorBase::GLOBAL, 1.);
            Phi1(i + 1).SetValuePreserve(0.);
            Phi1(i + 1).DerivativeResizeReset(&dof, i + 3, MapVectorBase::GLOBAL, 1.);
            X2(i + 1).SetValuePreserve(0.);
            X2(i + 1).DerivativeResizeReset(&dof, i + 6, MapVectorBase::GLOBAL, 1.);
            Phi2(i + 1).SetValuePreserve(0.);
            Phi2(i + 1).DerivativeResizeReset(&dof, i + 9, MapVectorBase::GLOBAL, 1.);            
        }
        
        const double start_calc = mbdyn_clock_time();
        
        gradVecAssTest::Euler123ToMatR(Phi1, R1);
        gradVecAssTest::Euler123ToMatR(Phi2, R2);
        
        Y = Transpose(R2) * Vec3(X1 + R1 * o1 - X2) - o2;
        
        calc += mbdyn_clock_time() - start_calc;
        
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < iNumDof; ++j) {
                jac(i + 1, j + 1) = Y(i + 1).dGetDerivativeGlobal(j);
            }
        }
    }
    
    std::cout << "calculation time: " << calc/N << "s\n";
    std::cout << "elapsed time: " << (mbdyn_clock_time() - start)/N << "s\n";
    
    std::cout << "X1=" << X1 << std::endl;
    std::cout << "X2=" << X2 << std::endl;
    std::cout << "Y=" << Y << std::endl;
    std::cout << "jac=\n";

    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < iNumDof; ++j) {
            std::cout << jac(i + 1, j + 1) << '\t';
        }
        std::cout << std::endl;
    }
}

Here is the call graph for this function:

template<typename T , index_type iRowCount>

void doVecAdd	(	const Vector< T, iRowCount > &	x,
		const Vector< T, iRowCount > &	y,
		Vector< T, iRowCount > &	z
	)

Definition at line 1297 of file matvectest.cc.

{
        z = x + y;
}

template<typename T , index_type iRowCount>

void doVecMul	(	const Vector< T, iRowCount > &	x,
		const Vector< T, iRowCount > &	y,
		Vector< T, iRowCount > &	z
	)

Definition at line 1303 of file matvectest.cc.

References grad::Vector< T, N_rows >::iGetNumRows().

{
        for (index_type i = 1; i <= z.iGetNumRows(); ++i)
        {
                z(i) = x(i) * y(i);
        }
}

Here is the call graph for this function:

template<typename T , index_type N>

void func	(	const Vector< T, N > &	a,
		const Vector< T, N > &	b,
		Vector< T, N > &	c
	)

Definition at line 200 of file matvectest.cc.

References grad::Dot(), and random1().

Referenced by callFunc().

                                                                         {
        doublereal r1, r2, r3;
        r1 = random1();
        r2 = random1();
        r3 = random1();
        const T d1 = Dot(a, b) * (2. / 1.5);
        c = (a * (r1 / 2) + b * (r2 / 2.)) - (a * (r3 * r1) - b * r1) + b * d1;
}

Here is the call graph for this function:

template<typename T >

void func	(	const T	a[],
		const T	b[],
		T	c[],
		index_type	N
	)

Definition at line 210 of file matvectest.cc.

References random1().

                                                         {
        doublereal r1, r2, r3;
        r1 = random1();
        r2 = random1();
        r3 = random1();

        T d1 = T(0.);

        for (int i = 0; i < N; ++i) {
                d1 += a[i] * b[i] * (2. / 1.5);
        }

        for (int i = 0; i < N; ++i) {
                c[i] = (a[i] * r1 + b[i] * r2) / 2. - (a[i] * r3 - b[i]) * r1 + b[i] * d1;
        }
}

Here is the call graph for this function:

template<typename T , index_type N>

void func2	(	const Matrix< T, N, N > &	A,
		const Vector< T, N > &	b,
		const Vector< T, N > &	c,
		Vector< T, N > &	d,
		doublereal	e,
		doublereal &	dt
	)

Definition at line 228 of file matvectest.cc.

References c, tic(), and grad::Transpose().

Referenced by callFunc2().

                                                                                                                                  {
        doublereal start, stop;
        tic(start);
        d = (Transpose(A) * Vector<T, N>(A * Vector<T, N>((b - c) * e)));
        tic(stop);
        dt += stop - start;
}

Here is the call graph for this function:

template<typename T >

void func2	(	const T *	A,
		const T *	b,
		const T *	c__,
		T *	d__,
		const doublereal &	e,
		doublereal &	dt
	)

Definition at line 237 of file matvectest.cc.

References tic().

                                                                                              {
        doublereal start, stop;
        tic(start);
    /* Parameter adjustments */
    --d__;
    --c__;
    --b;
    A -= 4;
    const T tmp1 = b[3] - c__[3];
    const T tmp2 = b[2] - c__[2];
    const T tmp3 = b[1] - c__[1];
    const T tmp4 = tmp3 * A[7];
    const T tmp5 = tmp1 * A[9];
    const T tmp6 = tmp1 * A[12];
    const T tmp7 = tmp1 * A[6];
    const T tmp8 = tmp2 * A[11];
    const T tmp10 = tmp3 * A[10];
    const T tmp11 = tmp6 + tmp8 + tmp10;
    const T tmp13 = tmp5 + tmp2 * A[8] + tmp4;
    const T tmp14 = tmp7 + tmp2 * A[5] + tmp3 * A[4];
    /* Function Body */
    d__[1] = (A[10] * tmp11 + A[7] * tmp13 + A[4] * tmp14) * e;
    d__[2] = (A[11] * tmp11 + A[8] * tmp13 + A[5] * tmp14) * e;
    d__[3] = (A[12] * tmp11 + A[9] * tmp13 + A[6] * tmp14) * e;

    tic(stop);
    dt += stop - start;
} /* func2_ */

Here is the call graph for this function:

void func2ad	(	const doublereal	A[3][3],
		const doublereal	b[3],
		const doublereal	c[3],
		doublereal	d[3],
		const doublereal &	e
	)

Referenced by callFunc2().

void func2ad ( const Matrix< doublereal, 3, 3 > &  A, const Vector< doublereal, 3 > &  b, const Vector< doublereal, 3 > &  c, Vector< doublereal, 3 > &  d, const doublereal &  e, LocalDofMap *  , doublereal &  dt  )

inline

Definition at line 297 of file matvectest.cc.

References grad::Vector< T, N_rows >::pGetVec(), and tic().

                                                                                                                                                                                                    {

        doublereal A_F[3][3];

        for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                        A_F[j][i] = A(i + 1, j + 1);
                }
        }

        doublereal start, stop;

        tic(start);
        func2ad_(A_F, b.pGetVec(), c.pGetVec(), d.pGetVec(), e);
        tic(stop);

        dt += stop - start;
}

Here is the call graph for this function:

template<index_type N_SIZE>

void func2ad ( const Matrix< Gradient< N_SIZE >, 3, 3 > &  A, const Vector< Gradient< N_SIZE >, 3 > &  b, const Vector< Gradient< N_SIZE >, 3 > &  c, Vector< Gradient< N_SIZE >, 3 > &  d, const doublereal &  e, LocalDofMap *  pDofMap, doublereal &  dt  )

inline

Definition at line 317 of file matvectest.cc.

References c, grad::Gradient< N_SIZE >::dGetDerivativeGlobal(), grad::Gradient< N_SIZE >::dGetValue(), grad::MapVectorBase::GLOBAL, nbdirsmax, and tic().

                                                                                                                                                                                                                                    {
        typedef typename MaxSizeCheck<N_SIZE <= index_type(nbdirsmax)>::CheckType check_nbdirsmax;

        doublereal A_F[3][3], Ad_F[3][3][nbdirsmax], b_F[3], bd_F[3][nbdirsmax], c_F[3], cd_F[3][nbdirsmax], d_F[3], dd_F[3][nbdirsmax];

        for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                        const Gradient<N_SIZE>& A_ij = A(i + 1, j + 1);
                        A_F[j][i] = A_ij.dGetValue();
                        for (index_type k = 0; k < N_SIZE; ++k) {
                                Ad_F[j][i][k] = A_ij.dGetDerivativeGlobal(k);
                        }
                }

                b_F[i] = b(i + 1).dGetValue();
                c_F[i] = c(i + 1).dGetValue();

                for (index_type k = 0; k < N_SIZE; ++k) {
                        bd_F[i][k] = b(i + 1).dGetDerivativeGlobal(k);
                        cd_F[i][k] = c(i + 1).dGetDerivativeGlobal(k);
                }
        }

        doublereal start, stop;
        tic(start);
        func2ad_dv_(A_F, Ad_F, b_F, bd_F, c_F, cd_F, d_F, dd_F, e, N_SIZE);
        tic(stop);
        dt += stop - start;

        for (int i = 0; i < 3; ++i) {
                d(i + 1).SetValuePreserve(d_F[i]);
                d(i + 1).DerivativeResizeReset(pDofMap, 0L, N_SIZE, MapVectorBase::GLOBAL, 0.);
                for (index_type k = 0; k < N_SIZE; ++k) {
                        d(i + 1).SetDerivativeGlobal(k, dd_F[i][k]);
                }
        }
}

Here is the call graph for this function:

void func2ad_dv	(	const doublereal	A[3][3],
		const doublereal	Ad[3][3][nbdirsmax],
		const doublereal	b[3],
		const doublereal	bd[3][nbdirsmax],
		const doublereal	c[3],
		const doublereal	cd[3][nbdirsmax],
		doublereal	d[3],
		doublereal	dd[3][nbdirsmax],
		const doublereal &	e,
		const integer &	nbdirs
	)

void func2addad_dv	(	const doublereal	x[],
		const doublereal	xd[],
		const doublereal	y[],
		const doublereal	yd[],
		doublereal	z[],
		doublereal	zd[],
		const integer &	n,
		const integer &	nbdirs
	)

Referenced by main().

void func2mulad_dv	(	const doublereal	x[],
		const doublereal	xd[],
		const doublereal	y[],
		const doublereal	yd[],
		doublereal	z[],
		doublereal	zd[],
		const integer &	n,
		const integer &	nbdirs
	)

Referenced by main().

template<typename T >

void func3	(	const Matrix< T, 3, 3 > &	R1,
		const Matrix< T, 3, 3 > &	R2,
		const Vector< T, 3 > &	a,
		const Vector< T, 3 > &	b,
		Vector< T, 3 > &	c,
		doublereal	e
	)

Definition at line 356 of file matvectest.cc.

References grad::Cross().

                                                                                                                                              {
        c = (Cross(R1 * a, R2 * b) + Cross(R2 * a, R1 * b)) * e;
}

Here is the call graph for this function:

template void func3< doublereal >	(	const Matrix< doublereal, 3, 3 > &	R1,
		const Matrix< doublereal, 3, 3 > &	R2,
		const Vector< doublereal, 3 > &	a,
		const Vector< doublereal, 3 > &	b,
		Vector< doublereal, 3 > &	c,
		doublereal	e
	)

template void func3< Gradient< 0 > >	(	const Matrix< Gradient< 0 >, 3, 3 > &	R1,
		const Matrix< Gradient< 0 >, 3, 3 > &	R2,
		const Vector< Gradient< 0 >, 3 > &	a,
		const Vector< Gradient< 0 >, 3 > &	b,
		Vector< Gradient< 0 >, 3 > &	c,
		doublereal	e
	)

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 3746 of file matvectest.cc.

References func2addad_dv(), func2mulad_dv(), GRADIENT_DEBUG, Mat3xN_test(), MatDynamic_test(), MatManip_test(), MatNxN_test(), MATVEC_DEBUG, NLoops, testInv(), testMatVec3(), testMatVecDouble2(), testMatVecProduct(), testMatVecProductGradient(), testScalarTypeTraits(), testSolve(), testSubVecAss(), and testSubVecAssMatVec().

                                 {
#ifdef HAVE_FEENABLEEXCEPT
        feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
#endif
        if (argc > 1) {
                NLoops = atol(argv[1]);
        }

    if (argc > 2) {
        NLoopsAss = atol(argv[2]);
    }

    if (NLoops < 1) {
        NLoops = 1;
    }

    if (NLoopsAss < 1) {
        NLoopsAss = 1;
    }

    std::cerr << "MatManip_test()\n";
    
    MatManip_test(NLoops);
    
        std::cerr << "---------------------------\ntestScalarTypeTraits()\n";
        testScalarTypeTraits();

        testMatVec<3>();
        testMatVec<6>();
        testMatVec<8>();
        testMatVec<12>();
        testMatVec<24>();

    std::cerr << "---------------------------\ntestMatVecGradient2<1>()\n";
        testMatVecGradient2<1>();

    std::cerr << "---------------------------\ntestMatVecGradient2<2>()\n";
        testMatVecGradient2<2>();

    std::cerr << "---------------------------\ntestMatVecGradient2<3>()\n";
        testMatVecGradient2<3>();

    std::cerr << "---------------------------\ntestMatVecGradient2<4>()\n";
        testMatVecGradient2<4>();

    std::cerr << "---------------------------\ntestMatVecGradient2<5>()\n";
        testMatVecGradient2<5>();

    std::cerr << "---------------------------\ntestMatVecGradient2<6>()\n";
        testMatVecGradient2<6>();

    std::cerr << "---------------------------\ntestMatVecGradient2<8>()\n";
        testMatVecGradient2<8>();

    std::cerr << "---------------------------\ntestMatVecGradient2<10>()\n";
        testMatVecGradient2<10>();

    std::cerr << "---------------------------\ntestMatVecGradient2<12>()\n";
        testMatVecGradient2<12>();

        std::cerr << "---------------------------\ntestMatVecDouble2()\n";
        testMatVecDouble2();

    std::cerr << "---------------------------\ntestMatVecProduct()\n";
    testMatVecProduct();
    std::cerr << "---------------------------\ntestMatVecProductGradient()\n";
    testMatVecProductGradient();

    std::cerr << "---------------------------\ntestMatVecProductGradient2()\n";
    testMatVecProductGradient2<1>(1, NLoops);
    testMatVecProductGradient2<4>(4, NLoops);
    testMatVecProductGradient2<8>(8, NLoops);
    testMatVecProductGradient2<12>(12, NLoops);
    testMatVecProductGradient2<32>(32, NLoops);
    testMatVecProductGradient2<64>(64, NLoops);
    testMatVecProductGradient2<0>(256, NLoops);
    
    std::cerr << "---------------------------\ntestMatVecCopy()\n";
    testMatVecCopy<1>();
    testMatVecCopy<3>();
    testMatVecCopy<5>();
    testMatVecCopy<9>();

#ifdef HAVE_BLITZ
    std::cerr << "---------------------------\ntestAssembly()\n";
    gradVecAssTest::testAssembly();
#endif

    std::cerr << "---------------------------\ntestMatVec3()\n";
    testMatVec3();

    std::cerr << "---------------------------\ntestSubVecAss()\n";
    testSubVecAss();

    std::cerr << "---------------------------\ntestSubVecAssMatVec()\n";
    testSubVecAssMatVec();

    std::cerr << "---------------------------\ntestInv()\n";
    testInv();

    std::cerr << "----------------------------\ntestSolve()\n";
    testSolve();

      {  
            const int N = argc > 3 ? atoi(argv[3]) : 1;
            const int nr = 4, nc = 4;
            Matrix<doublereal, nr, nc> A, B, C;

            for (int i = 1; i <= nr; ++i)
            {
                for (int j = 1; j <= nc; ++j)
                {
                    B(i, j) = rand();
                    C(i, j) = rand();
                }
            }

            clock_t start = clock();

            for (int i = 0; i < N; ++i)
            {
                A = B * C;
            }

            std::cerr << "time: " << double(clock()-start)/N/CLOCKS_PER_SEC << std::endl;
      }
      
    std::cerr << "----------------------------\ncppad_benchmark1<0>()\n";
    cppad_benchmark1<0>(NLoops);

    std::cerr << "----------------------------\ncppad_benchmark1<6>()\n";
    cppad_benchmark1<6>(NLoops);

    std::cerr << "----------------------------\ncppad_benchmark2<0>()\n";
    cppad_benchmark2<0>(NLoops);

    std::cerr << "----------------------------\ncppad_benchmark2<12>()\n";
    cppad_benchmark2<12>(NLoops);

    std::cerr << "----------------------------\ncppad_benchmark3<0>()\n";
    cppad_benchmark3<0>(NLoops);

    std::cerr << "----------------------------\ncppad_benchmark3<12>()\n";
    cppad_benchmark3<12>(NLoops);
      

    testVecOp<3, 2>(NLoops, 2, doVecAdd<Gradient<2>, 3>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<3, 4>(NLoops, 4, doVecAdd<Gradient<4>, 3>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<3, 8>(NLoops, 8, doVecAdd<Gradient<8>, 3>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<3, 16>(NLoops, 16, doVecAdd<Gradient<16>, 3>, __FC_DECL__(func2addad_dv), "add");

    testVecOp<12, 2>(NLoops, 2, doVecAdd<Gradient<2>, 12>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<12, 4>(NLoops, 4, doVecAdd<Gradient<4>, 12>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<12, 8>(NLoops, 8, doVecAdd<Gradient<8>, 12>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<12, 16>(NLoops, 16, doVecAdd<Gradient<16>, 12>, __FC_DECL__(func2addad_dv), "add");

    testVecOp<3, 2>(NLoops, 2, doVecMul<Gradient<2>, 3>, __FC_DECL__(func2mulad_dv), "mul");
        testVecOp<3, 4>(NLoops, 4, doVecMul<Gradient<4>, 3>, __FC_DECL__(func2mulad_dv), "mul");
        testVecOp<3, 8>(NLoops, 8, doVecMul<Gradient<8>, 3>, __FC_DECL__(func2mulad_dv), "mul");
        testVecOp<3, 16>(NLoops, 16, doVecMul<Gradient<16>, 3>, __FC_DECL__(func2mulad_dv), "mul");

    testVecOp<12, 2>(NLoops, 2, doVecAdd<Gradient<2>, 12>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<12, 4>(NLoops, 4, doVecAdd<Gradient<4>, 12>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<12, 8>(NLoops, 8, doVecAdd<Gradient<8>, 12>, __FC_DECL__(func2addad_dv), "add");
        testVecOp<12, 16>(NLoops, 16, doVecAdd<Gradient<16>, 12>, __FC_DECL__(func2addad_dv), "add");

    testVecOp<3, 0>(NLoops, 2, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<3, 0>(NLoops, 4, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<3, 0>(NLoops, 8, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<3, 0>(NLoops, 16, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<3, 0>(NLoops, 256, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<3, 0>(NLoops, 512, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<3, 0>(NLoops, 1024, doVecAdd<Gradient<0>, 3>, __FC_DECL__(func2addad_dv), "add");

    testVecOp<12, 0>(NLoops, 2, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<12, 0>(NLoops, 4, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<12, 0>(NLoops, 8, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<12, 0>(NLoops, 16, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<12, 0>(NLoops, 256, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<12, 0>(NLoops, 512, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");
    testVecOp<12, 0>(NLoops, 1024, doVecAdd<Gradient<0>, 12>, __FC_DECL__(func2addad_dv), "add");

    testVecOp<3, 0>(NLoops, 2, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<3, 0>(NLoops, 4, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<3, 0>(NLoops, 8, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<3, 0>(NLoops, 16, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<3, 0>(NLoops, 256, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<3, 0>(NLoops, 512, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<3, 0>(NLoops, 1024, doVecMul<Gradient<0>, 3>, __FC_DECL__(func2mulad_dv), "mul");

    testVecOp<12, 0>(NLoops, 2, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<12, 0>(NLoops, 4, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<12, 0>(NLoops, 8, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<12, 0>(NLoops, 16, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<12, 0>(NLoops, 256, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<12, 0>(NLoops, 512, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");
    testVecOp<12, 0>(NLoops, 1024, doVecMul<Gradient<0>, 12>, __FC_DECL__(func2mulad_dv), "mul");

    Mat3xN_test(3, NLoops);
    Mat3xN_test(10, NLoops);
    Mat3xN_test(100, NLoops);

    MatNxN_test(3, NLoops);
    MatNxN_test(10, NLoops);
    MatNxN_test(100, NLoops);

    Mat3xN_test_grad<4>(4, 3, NLoops);
    Mat3xN_test_grad<5>(5, 10, NLoops);
    Mat3xN_test_grad<0>(20, 100, NLoops);

    Mat3xNT_test_grad<4>(4, 3, NLoops);
    Mat3xNT_test_grad<5>(5, 10, NLoops);
    Mat3xNT_test_grad<0>(20, 100, NLoops);
    
    MatNxN_test_grad<4>(4, 3, NLoops);
    MatNxN_test_grad<5>(5, 10, NLoops);
    MatNxN_test_grad<0>(20, 100, NLoops);

    MatNxNT_test_grad<4>(4, 3, NLoops);
    MatNxNT_test_grad<5>(5, 10, NLoops);
    MatNxNT_test_grad<0>(20, 100, NLoops);
    
    MatDynamic_test(3, 4, NLoops);
    MatDynamic_test(10, 20, NLoops);
    
    MatDynamic_test_grad<3>(3, 4, 5, NLoops);
    MatDynamic_test_grad<5>(5, 10, 20, NLoops);
    MatDynamic_test_grad<0>(15, 20, 30, NLoops);

    MatDynamicT_test<10, 20>(10, 20, NLoops);
    MatDynamicT_test<DYNAMIC_SIZE, 20>(10, 20, NLoops);
    MatDynamicT_test<10, DYNAMIC_SIZE>(10, 20, NLoops);
    MatDynamicT_test<DYNAMIC_SIZE, DYNAMIC_SIZE>(10, 20, NLoops);

    MatDynamicT_test_grad<5, 10, 20>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<5, DYNAMIC_SIZE, 20>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<5, 10, DYNAMIC_SIZE>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<5, DYNAMIC_SIZE, DYNAMIC_SIZE>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<0, 10, 20>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<0, DYNAMIC_SIZE, 20>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<0, 10, DYNAMIC_SIZE>(5, 10, 20, NLoops);
    MatDynamicT_test_grad<0, DYNAMIC_SIZE, DYNAMIC_SIZE>(5, 10, 20, NLoops);
    
#ifdef NDEBUG
    std::cerr << "\nNo tests have been done" << std::endl;
#else
    std::cerr << "\nAll tests passed" << std::endl;
#endif

    std::cerr << "MATVEC_DEBUG=" << MATVEC_DEBUG << std::endl;
    std::cerr << "GRADIENT_DEBUG=" << GRADIENT_DEBUG << std::endl;

    return 0;
}

Here is the call graph for this function:

void Mat3xN_test	(	int	N,
		int	M
	)

Definition at line 3116 of file matvectest.cc.

References bCompare(), Mat3xN::iGetNumCols(), Mat3xN::iGetNumRows(), mbdyn_clock_time(), and grad::sqrt().

Referenced by main().

                               {
    Mat3xN A(N, 0.);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<doublereal, DYNAMIC_SIZE> x1(A.iGetNumCols()), x2(A.iGetNumCols()), x(A.iGetNumCols());

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x1(i) = i;
        x2(i) = 10 * i;
    }
    
    Vector<doublereal, 3> b = A * x;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < M; ++i) {
        x = x1 * 3. + x2 * 5.;
        b = A * x;
    }

    std::cerr << "Mat3xN: " << (mbdyn_clock_time() - start) / M << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        doublereal b_i = 0.;
        
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            b_i += A(i, j) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_SIZE>

void Mat3xN_test_grad	(	int	iNumDeriv,
		int	N,
		int	M
	)

Definition at line 3238 of file matvectest.cc.

References bCompare(), dof, grad::DYNAMIC_SIZE, Mat3xN::iGetNumCols(), Mat3xN::iGetNumRows(), grad::MapVectorBase::LOCAL, mbdyn_clock_time(), and grad::sqrt().

                                                   {
    assert((N_SIZE == 0) || (N_SIZE >= iNumDeriv));
    
    LocalDofMap dof;
    
    Mat3xN A(N, 0.);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> x1(A.iGetNumCols()), x2(A.iGetNumCols()), x(A.iGetNumCols());

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x1(i).SetValuePreserve(i);
        x1(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x1(i).SetDerivativeLocal(k, -1. * k);
        }
        
        x2(i).SetValuePreserve(10 * i);
        x2(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x2(i).SetDerivativeLocal(k, 10. * k);
        }
    }
    
    Vector<Gradient<N_SIZE>, 3> b = A * x;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < M; ++i) {
        x = x1 * 3. + x2 * 5.;
        b = A * x;
    }

    std::cerr << "Mat3xN * Gradient: " << (mbdyn_clock_time() - start) / M << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        Gradient<N_SIZE> b_i;
        
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            b_i += A(i, j) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_SIZE>

void Mat3xNT_test_grad	(	int	iNumDeriv,
		int	N,
		int	M
	)

Definition at line 3294 of file matvectest.cc.

References bCompare(), dof, grad::DYNAMIC_SIZE, Mat3xN::iGetNumCols(), Mat3xN::iGetNumRows(), grad::MapVectorBase::LOCAL, mbdyn_clock_time(), grad::sqrt(), and grad::Transpose().

                                                    {
    assert((N_SIZE == 0) || (N_SIZE >= iNumDeriv));
    
    LocalDofMap dof;
    
    Mat3xN A(N, 0.);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<Gradient<N_SIZE>, 3> x1(A.iGetNumRows()), x2(A.iGetNumRows());

    for (index_type i = 1; i <= x1.iGetNumRows(); ++i) {
        x1(i).SetValuePreserve(i);
        x1(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x1(i).SetDerivativeLocal(k, -1. * k);
        }
        
        x2(i).SetValuePreserve(10 * i);
        x2(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x2(i).SetDerivativeLocal(k, 10. * k);
        }
    }
    
    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> b = Transpose(A) * x1;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < M; ++i) {
        b = Transpose(A) * (x1 * 3. + x2 * 5.);
    }

    std::cerr << "Transpose(Mat3xN) * Gradient: " << (mbdyn_clock_time() - start) / M << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumCols(); ++i) {
        Gradient<N_SIZE> b_i;
        
        for (index_type j = 1; j <= A.iGetNumRows(); ++j) {
            b_i += A(j, i) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

void MatDynamic_test	(	index_type	iNumRows,
		index_type	iNumCols,
		index_type	iNumLoops
	)

Definition at line 3196 of file matvectest.cc.

References bCompare(), grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), mbdyn_clock_time(), and grad::sqrt().

Referenced by main().

                                                                                     {
    Matrix<doublereal, DYNAMIC_SIZE, DYNAMIC_SIZE> A(iNumRows, iNumCols);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<doublereal, DYNAMIC_SIZE> x1(A.iGetNumCols()), x2(A.iGetNumCols()), x(A.iGetNumCols());

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x1(i) = i;
        x2(i) = 10 * i;
    }
    
    Vector<doublereal, DYNAMIC_SIZE> b = A * x;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < iNumLoops; ++i) {
        x = x1 * 3. + x2 * 5.;
        b = A * x;
    }

    std::cerr << "Matrix<DYNAMIC_SIZE, DYNAMIC_SIZE>: " << (mbdyn_clock_time() - start) / iNumLoops << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        doublereal b_i = 0.;
        
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            b_i += A(i, j) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_SIZE>

void MatDynamic_test_grad	(	index_type	iNumDeriv,
		index_type	iNumRows,
		index_type	iNumCols,
		int	iNumLoops
	)

Definition at line 3461 of file matvectest.cc.

References bCompare(), dof, grad::DYNAMIC_SIZE, grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), grad::MapVectorBase::LOCAL, mbdyn_clock_time(), and grad::sqrt().

{
    assert((N_SIZE == 0) || (N_SIZE >= iNumDeriv));
    
    LocalDofMap dof;
    
    Matrix<Gradient<N_SIZE>, DYNAMIC_SIZE, DYNAMIC_SIZE> A(iNumRows, iNumCols);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j).SetValuePreserve(100 * i + j);
            A(i, j).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
            
            for (index_type k = 0; k < iNumDeriv; ++k) {
                A(i, j).SetDerivativeLocal(k, 1000. * i + 100. * j + k + 1);
            }
        }
    }

    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> x1(A.iGetNumCols()), x2(A.iGetNumCols()), x(A.iGetNumCols());

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x1(i).SetValuePreserve(i);
        x1(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x1(i).SetDerivativeLocal(k, -1. * k);
        }
        
        x2(i).SetValuePreserve(10 * i);
        x2(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x2(i).SetDerivativeLocal(k, 10. * k);
        }
    }
    
    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> b = A * x;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < iNumLoops; ++i) {
        x = x1 * 3. + x2 * 5.;
        b = A * x;
    }

    std::cerr << "Matrix<Gradient,DYNAMIC_SIZE,DYNAMIC_SIZE> * Gradient: " << (mbdyn_clock_time() - start) / iNumLoops << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        Gradient<N_SIZE> b_i;
        
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            b_i += A(i, j) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_ROWS, index_type N_COLS>

void MatDynamicT_test	(	index_type	iNumRows,
		index_type	iNumCols,
		int	iNumLoops
	)

Definition at line 3523 of file matvectest.cc.

References bCompare(), grad::DYNAMIC_SIZE, grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), mbdyn_clock_time(), grad::sqrt(), and grad::Transpose().

{
    assert((N_ROWS == iNumRows) || (N_ROWS == DYNAMIC_SIZE));
    assert((N_COLS == iNumCols) || (N_COLS == DYNAMIC_SIZE));
     
    Matrix<doublereal, N_ROWS, N_COLS> A(iNumRows, iNumCols);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<doublereal, N_ROWS> x1(A.iGetNumRows()), x2(A.iGetNumRows());

    for (index_type i = 1; i <= x1.iGetNumRows(); ++i) {
        x1(i) = i;
        x2(i) = 10 * i;
    }
    
    Vector<doublereal, N_COLS> b = Transpose(A) * x1;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < iNumLoops; ++i) {
        b = Transpose(A) * (x1 * 3. + x2 * 5.);
    }

    std::cerr << "Transpose(Matrix<doublereal>,"
              << iNumRows << "(" << N_ROWS << ")," << iNumCols << "(" << N_COLS << ")"
              << ">) * Vector<doublereal>,"
              << iNumRows << "(" << N_ROWS << ")>: " << (mbdyn_clock_time() - start) / iNumLoops << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumCols(); ++i) {
        doublereal b_i = 0.;
        
        for (index_type j = 1; j <= A.iGetNumRows(); ++j) {
            b_i += A(j, i) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_DERIV, index_type N_ROWS, index_type N_COLS>

void MatDynamicT_test_grad	(	index_type	iNumDeriv,
		index_type	iNumRows,
		index_type	iNumCols,
		int	iNumLoops
	)

Definition at line 3570 of file matvectest.cc.

References bCompare(), dof, grad::DYNAMIC_SIZE, grad::Matrix< T, N_rows, N_cols >::GetCol(), grad::Matrix< T, N_rows, N_cols >::GetRow(), grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), grad::MapVectorBase::LOCAL, mbdyn_clock_time(), grad::sqrt(), and grad::Transpose().

{
    assert((N_DERIV == 0) || (N_DERIV >= iNumDeriv));
    assert((N_ROWS == iNumRows) || (N_ROWS == DYNAMIC_SIZE));
    assert((N_COLS == iNumCols) || (N_COLS == DYNAMIC_SIZE));
    
    LocalDofMap dof;
    
    Matrix<Gradient<N_DERIV>, N_ROWS, N_COLS> A(iNumRows, iNumCols);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j).SetValuePreserve(100 * i + j);
            A(i, j).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
            
            for (index_type k = 0; k < iNumDeriv; ++k) {
                A(i, j).SetDerivativeLocal(k, 1000. * i + 100. * j + k + 1);
            }
        }
    }

    Vector<Gradient<N_DERIV>, N_ROWS> x1(A.iGetNumRows()), x2(A.iGetNumRows());

    for (index_type i = 1; i <= x1.iGetNumRows(); ++i) {
        x1(i).SetValuePreserve(i);
        x1(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x1(i).SetDerivativeLocal(k, -1. * k);
        }
        
        x2(i).SetValuePreserve(10 * i);
        x2(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x2(i).SetDerivativeLocal(k, 10. * k);
        }
    }
    
    Vector<Gradient<N_DERIV>, N_COLS> b = Transpose(A) * x1;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < iNumLoops; ++i) {
        b = Transpose(A) * (x1 * 3. + x2 * 5.);
    }

    std::cerr << "Transpose(Matrix<Gradient<" << iNumDeriv << "(" << N_DERIV << ")"
              << ">," << iNumRows << "(" << N_ROWS << ")," << iNumCols << "(" << N_COLS << ")"
              << ">) * Vector<Gradient<" << iNumDeriv << "(" << N_DERIV << ")"
              << ">," << iNumRows << "(" << N_ROWS << ")>: " << (mbdyn_clock_time() - start) / iNumLoops << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumCols(); ++i) {
        Gradient<N_DERIV> b_i;
        
        for (index_type j = 1; j <= A.iGetNumRows(); ++j) {
            b_i += A(j, i) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            assert(bCompare(A(i, j), A.GetRow(i)(j), 0.));
            assert(bCompare(A(i, j), A.GetCol(j)(i), 0.));
            assert(bCompare(A(i, j), Transpose(A)(j, i), 0.));
            assert(bCompare(A(i, j), Transpose(A).GetCol(i)(j), 0.));
            assert(bCompare(A(i, j), Transpose(A).GetRow(j)(i), 0.));
        }
    }
}

Here is the call graph for this function:

void MatManip_test

(

int

NLoops

)

Definition at line 3645 of file matvectest.cc.

References bCompare(), Eye3, grad::fabs(), M_PI, grad::MatCrossVec(), CGR_Rot::MatG, grad::MatGVec(), CGR_Rot::MatR, grad::MatRVec(), NLoops, Vec3::Norm(), grad::sqrt(), RotManip::VecRot(), and grad::VecRotMat().

Referenced by main().

                               {
    const doublereal tol1 = 10. * std::numeric_limits<doublereal>::epsilon();
    const doublereal tol2 = 1e-4 * sqrt(tol1);
    const doublereal alpha = 1. - sqrt(tol2);
    
    srand(0);
    
    for (int k = 0; k < 3 * NLoops; ++k) {
        Vector<doublereal, 3> g0;
        Vec3 g0_;
        
        for (int i = 1; i <= 3; ++i) {
            g0_(i) = (2. * rand() / RAND_MAX - 1.);
        }

        g0_ *= (2. * rand() / RAND_MAX - 1.) * alpha * M_PI / g0_.Norm();

        if (k >= NLoops && k < 2 * NLoops) {
            g0_ *= sqrt(std::numeric_limits<doublereal>::epsilon());
        } else if (k >= 2 * NLoops) {
            g0_ *= std::numeric_limits<doublereal>::epsilon();
        }
        
        g0 = g0_;
        
        const Mat3x3 G1(CGR_Rot::MatG, g0_);
        const Mat3x3 R1(CGR_Rot::MatR, g0_);
        const Mat3x3 X1(1., g0_);
        const Matrix<doublereal, 3, 3> G2(MatGVec(g0));
        Matrix<doublereal, 3, 3> R2(MatRVec(g0));
        const Matrix<doublereal, 3, 3> X2(MatCrossVec(g0, 1.));
        
        for (index_type i = 1; i <= 3; ++i) {
            for (index_type j = 1; j <= 3; ++j) {
                assert(bCompare(G1(i, j), G2(i, j), tol1));
                assert(bCompare(R1(i, j), R2(i, j), tol1));
                assert(bCompare(X1(i, j), X2(i, j), tol1));
            }
        }

        if (k == 0) {
            R2 = ::Eye3;
        } else if (k == 1) {
            R2(1, 1) = -1;
            R2(1, 2) = 0;
            R2(1, 3) = 0;
            R2(2, 1) = 0;
            R2(2, 2) = -1;
            R2(2, 3) = 0;
            R2(3, 1) = 0;
            R2(3, 2) = 0;
            R2(3, 3) = 1;
        } else if (k == 2) {
            R2(1, 1) = -1;
            R2(1, 2) = 0;
            R2(1, 3) = 0;
            R2(2, 1) = 0;
            R2(2, 2) = 1;
            R2(2, 3) = 0;
            R2(3, 1) = 0;
            R2(3, 2) = 0;
            R2(3, 3) = -1;            
        } else if (k == 3) {
            R2(1, 1) = 1;
            R2(1, 2) = 0;
            R2(1, 3) = 0;
            R2(2, 1) = 0;
            R2(2, 2) = -1;
            R2(2, 3) = 0;
            R2(3, 1) = 0;
            R2(3, 2) = 0;
            R2(3, 3) = -1;
        } else if (k == 4) {            
            R2(1,1)=-2.2841125213377644e-01;
            R2(1,2)=9.5997603033895429e-01;
            R2(1,3)=1.6209355654480440e-01;
            R2(2,1)=9.5997603033895418e-01;
            R2(2,2)=1.9435901751267337e-01;
            R2(2,3)=2.0166951551033063e-01;
            R2(3,1)=1.6209355654480423e-01;
            R2(3,2)=2.0166951551033083e-01;
            R2(3,3)=-9.6594776537889704e-01;         
        }
        
        Mat3x3 R2_;

        for (index_type i = 1; i <= 3; ++i) {
            for (index_type j = 1; j <= 3; ++j) {
                R2_(i, j) = R2(i, j);
            }
        }
        
        const Vec3 g1(RotManip::VecRot(R2_));
        const Vector<doublereal, 3> g2(VecRotMat(R2));
        
        for (index_type i = 1; i <= 3; ++i) {
            assert(bCompare(g1(i), g2(i), std::numeric_limits<doublereal>::epsilon()) || bCompare(fabs(g1(i) - g2(i)), 2 * M_PI, std::numeric_limits<doublereal>::epsilon()));
        }
    }
}

Here is the call graph for this function:

void MatNxN_test	(	int	N,
		int	M
	)

Definition at line 3156 of file matvectest.cc.

References bCompare(), MatNxN::iGetNumCols(), MatNxN::iGetNumRows(), mbdyn_clock_time(), and grad::sqrt().

Referenced by main().

                               {
    MatNxN A(N, 0.);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<doublereal, DYNAMIC_SIZE> x1(A.iGetNumCols()), x2(A.iGetNumCols()), x(A.iGetNumCols());

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x1(i) = i;
        x2(i) = 10 * i;
    }
    
    Vector<doublereal, DYNAMIC_SIZE> b = A * x;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < M; ++i) {
        x = x1 * 3. + x2 * 5.;
        b = A * x;
    }

    std::cerr << "MatNxN: " << (mbdyn_clock_time() - start) / M << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        doublereal b_i = 0.;
        
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            b_i += A(i, j) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_SIZE>

void MatNxN_test_grad	(	int	iNumDeriv,
		int	N,
		int	M
	)

Definition at line 3404 of file matvectest.cc.

References bCompare(), dof, grad::DYNAMIC_SIZE, MatNxN::iGetNumCols(), MatNxN::iGetNumRows(), grad::MapVectorBase::LOCAL, mbdyn_clock_time(), and grad::sqrt().

{
    assert((N_SIZE == 0) || (N_SIZE >= iNumDeriv));
    
    LocalDofMap dof;
    
    MatNxN A(N, 0.);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> x1(A.iGetNumCols()), x2(A.iGetNumCols()), x(A.iGetNumCols());

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x1(i).SetValuePreserve(i);
        x1(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x1(i).SetDerivativeLocal(k, -1. * k);
        }
        
        x2(i).SetValuePreserve(10 * i);
        x2(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x2(i).SetDerivativeLocal(k, 10. * k);
        }
    }
    
    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> b = A * x;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < M; ++i) {
        x = x1 * 3. + x2 * 5.;
        b = A * x;
    }

    std::cerr << "MatNxN * Gradient: " << (mbdyn_clock_time() - start) / M << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        Gradient<N_SIZE> b_i;
        
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            b_i += A(i, j) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<index_type N_SIZE>

void MatNxNT_test_grad	(	int	iNumDeriv,
		int	N,
		int	M
	)

Definition at line 3349 of file matvectest.cc.

References bCompare(), dof, grad::DYNAMIC_SIZE, MatNxN::iGetNumCols(), MatNxN::iGetNumRows(), grad::MapVectorBase::LOCAL, mbdyn_clock_time(), grad::sqrt(), and grad::Transpose().

                                                    {
    assert((N_SIZE == 0) || (N_SIZE >= iNumDeriv));
    
    LocalDofMap dof;
    
    MatNxN A(N, 0.);

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
            A(i, j) = 100 * i + j;
        }
    }

    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> x1(A.iGetNumRows()), x2(A.iGetNumRows());

    for (index_type i = 1; i <= x1.iGetNumRows(); ++i) {
        x1(i).SetValuePreserve(i);
        x1(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x1(i).SetDerivativeLocal(k, -1. * k);
        }
        
        x2(i).SetValuePreserve(10 * i);
        x2(i).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::LOCAL, 0.);
        
        for (index_type k = 0; k < iNumDeriv; ++k) {
            x2(i).SetDerivativeLocal(k, 10. * k);
        }
    }
    
    Vector<Gradient<N_SIZE>, DYNAMIC_SIZE> b = Transpose(A) * x1;
    
    const double start = mbdyn_clock_time();
    
    for (int i = 0; i < M; ++i) {
        b = Transpose(A) * (x1 * 3. + x2 * 5.);
    }

    std::cerr << "Transpose(MatNxN) * Gradient: " << (mbdyn_clock_time() - start) / M << "s\n";

    const doublereal tol = sqrt(std::numeric_limits<doublereal>::epsilon());
    
    for (index_type i = 1; i <= A.iGetNumCols(); ++i) {
        Gradient<N_SIZE> b_i;
        
        for (index_type j = 1; j <= A.iGetNumRows(); ++j) {
            b_i += A(j, i) * (3. * x1(j) + 5. * x2(j));
        }

        assert(bCompare(b_i, b(i), tol));
    }
}

Here is the call graph for this function:

template<typename T , index_type N_rows>

T Norm_1 ( const Vector< T, N_rows > &  u )

inline

Definition at line 1274 of file matvectest.cc.

References grad::Dot(), and grad::sqrt().

                                   {
        return sqrt(Dot(u, u));
}

Here is the call graph for this function:

doublereal random1

(

)

Definition at line 88 of file matvectest.cc.

Referenced by callFunc2(), func(), and testMatVec3().

                     {
        return 2 * (doublereal(rand()) / RAND_MAX) - 1;
}

void testInv

(

)

Definition at line 2807 of file matvectest.cc.

References bCompare(), grad::Det(), grad::Inv(), grad::sqrt(), and grad::Tabular().

Referenced by main().

               {
        Matrix<doublereal, 2, 2> A;

        A(1, 1) =  0.658371336838182;
        A(1, 2) = 0.733036075010795;
        A(2, 1) = 0.483830962404444;
        A(2, 2) = 0.395950513263802;

        const doublereal detA = Det(A);
        const Matrix<doublereal, 2, 2> invA = Inv(A);
        const Matrix<doublereal, 2, 2> B1 = invA * A;
        const Matrix<doublereal, 2, 2> B2 = A * invA;
        std::cout << "A=" << Tabular(A) << std::endl;
        std::cout << "Inv(A)=" << Tabular(invA) << std::endl;
        std::cout << "Det(A)=" << detA << std::endl;
        std::cout << "Inv(A)*A=" << Tabular(B1) << std::endl;
        std::cout << "A * Inv(A)=" << Tabular(B2) << std::endl;

        const doublereal dTol = sqrt(std::numeric_limits<doublereal>::epsilon());

        assert(bCompare(B1(1, 1), 1., dTol));
        assert(bCompare(B2(1, 1), 1., dTol));
        assert(bCompare(B1(2, 2), 1., dTol));
        assert(bCompare(B2(2, 2), 1., dTol));
        assert(bCompare(B1(1, 2), 0., dTol));
        assert(bCompare(B2(1, 2), 0., dTol));
        assert(bCompare(B1(2, 1), 0., dTol));
        assert(bCompare(B2(2, 1), 0., dTol));

        assert(bCompare(invA(1, 1), -4.21299780160758, dTol));
        assert(bCompare(invA(2, 2), -7.00521126207687, dTol));
        assert(bCompare(invA(1, 2), 7.79965998039245, dTol));
        assert(bCompare(invA(2, 1), 5.14806449967026, dTol));
        assert(bCompare(detA, -0.0939830809103952, dTol));
}

Here is the call graph for this function:

template<index_type N>

void testMatVec

(

)

Definition at line 2883 of file matvectest.cc.

References bCompare().

                  {
        doublereal c_MatVecGradient[N], cd_MatVecGradient[N][N];
        doublereal c_MatVecDouble[N];

        std::cerr << "---------------------------\ntestMatVecGradient<" << N << ">()\n";
    testMatVecGradient<N>(c_MatVecGradient, cd_MatVecGradient);

    std::cerr << "---------------------------\ntestMatVecDouble<" << N << ">()\n";
    testMatVecDouble<N>(c_MatVecDouble);

#ifdef HAVE_BLITZ
        doublereal c_MatVecDoubleBlitz[N];
        doublereal c_MatVecGradientBlitz[N], cd_MatVecGradientBlitz[N][N];
    std::cerr << "---------------------------\ntestMatVecDoubleBlitz<" << N << ">()\n";
    testMatVecDoubleBlitz<N>(c_MatVecDoubleBlitz);

        std::cerr << "---------------------------\ntestMatVecGradientBlitz<" << N << ">()\n";
    testMatVecGradientBlitz<N>(c_MatVecGradientBlitz, cd_MatVecGradientBlitz);
#endif // HAVE_BLITZ

    const doublereal dTol = 10 * std::numeric_limits<scalar_deriv_type>::epsilon();

    for (int i = 0; i < N; ++i) {
        assert(bCompare(c_MatVecGradient[i], c_MatVecDouble[i], dTol));
#ifdef HAVE_BLITZ
        assert(bCompare(c_MatVecGradient[i], c_MatVecDoubleBlitz[i], dTol));
        assert(bCompare(c_MatVecGradient[i], c_MatVecGradientBlitz[i], dTol));

        for (int j = 0; j < N; ++j) {
                assert(bCompare(cd_MatVecGradient[i][j], cd_MatVecGradientBlitz[i][j]));
        }
#endif
    }
}

Here is the call graph for this function:

void testMatVec3

(

)

Definition at line 2632 of file matvectest.cc.

References MatCross, MatCrossCross, grad::MatCrossCrossVec(), grad::MatCrossVec(), CGR_Rot::MatG, grad::MatGVec(), NLoops, and random1().

Referenced by main().

                   {
        srand(0);

        for (int n = 0; n < NLoops; ++n) {
                Vector<doublereal, 3> g1, v1;
                Vec3 g2, v2;
                for (index_type i = 1; i <= 3; ++i) {
                        g2(i) = g1(i) = random1() * 1e-1;
                        v2(i) = v1(i) = random1() * 1e1;
                }

                Matrix<doublereal, 3, 3> R1(MatGVec(g1));
                Mat3x3 R2(CGR_Rot::MatG, g2);
                Matrix<doublereal, 3, 3> C1(MatCrossVec(g1));
                Mat3x3 C2(MatCross, g2);
                Matrix<doublereal, 3, 3> CC1(MatCrossCrossVec(g1));
                Mat3x3 CC2(MatCrossCross, g2, g2);
                Vector<doublereal, 3> CCv1 = CC1 * v1;
                Vec3 CCv2 = CC2 * v2;
                Vector<doublereal, 3> CCv3 = CC1 * v2;
                Vector<doublereal, 3> CCv4 = CC2 * v1;
                Matrix<doublereal, 3, 3> A1 = R1 * C1;
                Mat3x3 A2 = R2 * C2;
                Matrix<doublereal, 3, 3> A3 = R1 * C2;
                Matrix<doublereal, 3, 3> A4 = R2 * C1;
                Vector<doublereal, 3> v3(v1(1), v1(2), v1(3));
                Vector<doublereal, 2> v4(v1(1), v1(2));
        Matrix<doublereal, 3, 3> X1(MatCrossVec(g1, 1.));
        Mat3x3 X2(1., g2);

                if (n == 0) {
                        std::cout << "g1=" << std::endl << g1 << std::endl;
                        std::cout << "g2=" << std::endl << g2 << std::endl;
                        std::cout << "R1=" << std::endl << R1 << std::endl;
                        std::cout << "R2=" << std::endl << R2 << std::endl;
                        std::cout << "C1=" << std::endl << C1 << std::endl;
                        std::cout << "C2=" << std::endl << C2 << std::endl;
                        std::cout << "CC1=" << std::endl << CC1 << std::endl;
                        std::cout << "CC2=" << std::endl << CC2 << std::endl;
                        std::cout << "CCv1=" << std::endl << CCv1 << std::endl;
                        std::cout << "CCv2=" << std::endl << CCv2 << std::endl;
                        std::cout << "CCv3=" << std::endl << CCv3 << std::endl;
                        std::cout << "CCv4=" << std::endl << CCv4 << std::endl;
                        std::cout << "A1=" << std::endl << A1 << std::endl;
                        std::cout << "A2=" << std::endl << A2 << std::endl;
                        std::cout << "A3=" << std::endl << A3 << std::endl;
                        std::cout << "A4=" << std::endl << A4 << std::endl;
                }

        const doublereal dTol = 10*std::numeric_limits<scalar_deriv_type>::epsilon();

                for (index_type i = 1; i <= 3; ++i) {
                        for (index_type j = 1; j <= 3; ++j) {
                                assert(std::abs(R1(i, j) - R2(i, j)) < dTol);
                                assert(std::abs(C1(i, j) - C2(i, j)) < dTol);
                                assert(std::abs(CC1(i, j) - CC2(i, j)) < dTol);
                                assert(std::abs(A1(i, j) - A2(i, j)) < dTol);
                                assert(std::abs(A1(i, j) - A3(i, j)) < dTol);
                                assert(std::abs(A1(i, j) - A4(i, j)) < dTol);
                assert(std::abs(X1(i, j) - X2(i, j)) < dTol);
                        }
                        assert(std::abs(CCv1(i) - CCv2(i)) < dTol);
                        assert(std::abs(CCv3(i) - CCv1(i)) < dTol);
                        assert(std::abs(CCv4(i) - CCv1(i)) < dTol);
                        assert(v3(i) == v1(i));

                        if (i < 3) {
                                assert(v4(i) == v1(i));
                        }
                }
        }
}

Here is the call graph for this function:

template<index_type N_rows>

void testMatVecCopy

(

)

Definition at line 1750 of file matvectest.cc.

References grad::Gradient< N_SIZE >::dGetDerivativeGlobal(), grad::Gradient< N_SIZE >::dGetValue(), grad::MapVectorBase::GLOBAL, grad::Vector< T, N_rows >::iGetNumRows(), and grad::Vector< T, N_rows >::Reset().

                      {
        std::cerr << "---------------------------\ntestMatVecCopy<" << N_rows << ">()\n";

        LocalDofMap dofMap1;
        Vector<Gradient<0>, N_rows> v;

        for (int i = 1; i <= v.iGetNumRows(); ++i) {
                v(i).SetValuePreserve(i);
                v(i).DerivativeResizeReset(&dofMap1, i, MapVectorBase::GLOBAL, i);
        }

        std::cout << "v=" << std::endl << v << std::endl;

        for (int i = 2; i <= v.iGetNumRows() - 1; ++i) {
                LocalDofMap dofMap2;
                Gradient<4> g1(v(1), &dofMap2), gi(v(i), &dofMap2), gim1(v(i - 1), &dofMap2), gip1(v(i + 1), &dofMap2);
                Gradient<4> x = 1000 * g1 + 100 * gi + 10 * gip1 + gim1;
                std::cout << "x(" << i << ")=" << x << std::endl;

                assert(x.dGetValue() == 1000 + 100 * i + 10 * (i + 1) + i - 1);

                if (i == 2) {
                        assert(x.dGetDerivativeGlobal(1) == 1001.);
                } else {
                        assert(x.dGetDerivativeGlobal(1) == 1000.);
                        assert(x.dGetDerivativeGlobal(i - 1) == i - 1);
                }

                assert(x.dGetDerivativeGlobal(i) == 100. * i);
                assert(x.dGetDerivativeGlobal(i + 1) == 10. * (i + 1));
        }

        v.Reset();
        std::cout << "v=" << v << std::endl;
}

Here is the call graph for this function:

template<index_type N>

void testMatVecDouble

(

doublereal

c_C[N]

)

Definition at line 635 of file matvectest.cc.

References a, bCompare(), c, callFunc(), grad::Dot(), grad::sqrt(), and grad::Sum().

                                         {
    Vector<doublereal, N> a, b;
    
    for (index_type i = 0; i < N; ++i) {
        a(i + 1) = 100*(i + 1);
        b(i + 1) = 1000*(i + 10);
    }
    
    Vector<doublereal, N> c, c1;

    callFunc(0, a, b, c, c1);

    Vector<doublereal, N> d = -c;

    const doublereal dTol = sqrt(std::numeric_limits<scalar_deriv_type>::epsilon());

    for (int i = 1; i <= N; ++i) {
        assert(bCompare(c(i), c1(i), dTol));
        assert(d(i) == -c(i));
    }

    for (int i = 0; i < N; ++i) {
        c_C[i] = c(i + 1);
    }

    doublereal s1 = 0.;

    for (int i = 0; i < N; ++i) {
        s1 += c(i + 1);
    }

    doublereal s2 = Sum(c);

    assert(bCompare(s2, s1, dTol));

    doublereal s3 = Sum((a * 2 - b * 3) / 1.5 + c);

    doublereal s4 = 0.;

    for (int i = 1; i <= N; ++i) {
        s4 += (a(i) * 2 - b(i) * 3) / 1.5 + c(i);
    }

    assert(bCompare(s2, s1, dTol));
    assert(bCompare(s4, s3, dTol));

    doublereal d1 = Dot(a, b);
    doublereal d2 = 0.;

    for (int i = 1; i <= N; ++i) {
        d2 += a(i) * b(i);
    }

    assert(bCompare(d2, d1, std::numeric_limits<scalar_deriv_type>::epsilon()));

    std::cout << "s1=" << s1 << std::endl;
    std::cout << "s2=" << s2 << std::endl;
    std::cout << "s3=" << s3 << std::endl;
    std::cout << "s4=" << s4 << std::endl;
    std::cout << "d1=" << d1 << std::endl;
    std::cout << "d2=" << d2 << std::endl;

    const Vector<doublereal, N> d5 = (d + c + c1) * (3.5 / (c(1) + c(2)) * c(1) / 2. * (c(1) - c(3)) / c(2));

    d += c + c1;

    d *= 3.5;

    d /= c(1) + c(2);

    d *= c(1);

    d /= 2.;

    d *= c(1) - c(3);

    d /= c(2);

    for (int i = 1; i <= N; ++i) {
        assert(bCompare(d(i), d5(i), dTol));
    }

    std::cout << "d=" << d << std::endl;
    std::cout << "d5=" << d5 << std::endl;
}

Here is the call graph for this function:

void testMatVecDouble2

(

)

Definition at line 832 of file matvectest.cc.

References c, and callFunc2().

Referenced by main().

                         {
    Matrix<doublereal, 3, 3> A;
    Vector<doublereal, 3> b, c;

    for (index_type i = 0; i < 3; ++i) {
        b(i + 1)=(100*(i + 1));
        c(i + 1)=(1000*(i + 10));

        for (index_type j = 0; j < 3; ++j) {
            A(i + 1, j + 1) = (100*(i + 1) + j + 1);
        }
    }

    Vector<doublereal, 3> d, d_C, d_F;

    callFunc2(0, A, b, c, d, d_C, d_F);
}

Here is the call graph for this function:

template<index_type N>

void testMatVecGradient	(	doublereal	c_C[N],
		doublereal	cd_C[N][N]
	)

Definition at line 474 of file matvectest.cc.

References a, bCompare(), grad::Gradient< N_SIZE >::bIsEqual(), c, callFunc(), dof, grad::Dot(), grad::MapVectorBase::GLOBAL, grad::Gradient< N_SIZE >::SetValue(), grad::sqrt(), and grad::Sum().

                                                                  {
    LocalDofMap dof;
    Vector<Gradient<N>, N> a, b;
    
    for (index_type i = 0; i < N; ++i) {
        a(i + 1).SetValuePreserve(100*(i + 1));
        b(i + 1).SetValuePreserve(1000*(i + 10));

        a(i + 1).DerivativeResizeReset(&dof, i, MapVectorBase::GLOBAL, 0.);
        a(i + 1).SetDerivativeGlobal(i, -1. - 10.);
        b(i + 1).DerivativeResizeReset(&dof, i, MapVectorBase::GLOBAL, 0.);
        b(i + 1).SetDerivativeGlobal(i, -2. - 20.);
    }
    
    Vector<Gradient<N>, N> c, c1;

    callFunc(&dof, a, b, c, c1);

    Vector<Gradient<N>, N> d = -c;

    std::cout << "c=" << std::endl << c << std::endl;
    std::cout << "d=" << std::endl << d << std::endl;

    const doublereal dTol = sqrt(std::numeric_limits<scalar_deriv_type>::epsilon());

    for (int i = 1; i <= N; ++i) {
        assert(bCompare(c(i), c1(i), dTol));
        assert(bCompare(d(i), Gradient<N>(-c(i))));
    }

    for (int i = 0; i < N; ++i) {
        c_C[i] = c(i + 1).dGetValue();

        for (int j = 0; j < N; ++j) {
                cd_C[i][j] = c(i + 1).dGetDerivativeGlobal(j);
        }
    }

    Gradient<N> s1;

    for (int i = 1; i <= N; ++i) {
        s1 += c(i);
    }

    Gradient<N> s2 = Sum(c);

    Gradient<N> s4;

    for (int i = 1; i <= N; ++i) {
        s4 += (a(i) * 2 - b(i) * 3) / 1.5 + c(i);
    }

    Gradient<N> s3 = Sum((a * 2 - b * 3) / 1.5 + c);

    Gradient<N> d1 = Dot(a, b);
    Gradient<N> d2;

    for (int i = 1; i <= N; ++i) {
        d2 += a(i) * b(i);
    }

    Gradient<N> d3 = Dot(a + b, a - b);
    Gradient<N> d4;

    for (int i = 1; i <= N; ++i) {
        d4 += (a(i) + b(i)) * (a(i) - b(i));
    }

    assert(d3.bIsEqual(d4));

    d3 = Dot(a + b, b);

    d4.SetValue(0.);

    for (int i = 1; i <= N; ++i) {
        d4 += (a(i) + b(i)) * b(i);
    }

    assert(d3.bIsEqual(d4));

    d3 = Dot(b, a + b);

    std::cout << "s1=" << s1 << std::endl;
    std::cout << "s2=" << s2 << std::endl;
    std::cout << "s3=" << s3 << std::endl;
    std::cout << "s4=" << s4 << std::endl;
    std::cout << "d1=" << d1 << std::endl;
    std::cout << "d2=" << d2 << std::endl;
    std::cout << "d3=" << d3 << std::endl;
    std::cout << "d4=" << d4 << std::endl;

    assert(bCompare(d3, d4, dTol));
    assert(bCompare(s2, s1, dTol));
    assert(bCompare(s4, s3, dTol));
    assert(bCompare(d2, d1, dTol));

    const Vector<Gradient<N>, N> d5 = (d + c + c1) * (3.5 / (c(1) + c(2)) * c(1) / 2. * (c(1) - c(3)) / c(2));

    d += c + c1;

    d *= 3.5;

    d /= c(1) + c(2);

    d *= c(1);

    d /= 2.;

    d *= c(1) - c(3);

    d /= c(2);

    std::cout << "d=" << d << std::endl;
    std::cout << "d5=" << d5 << std::endl;

    for (int i = 1; i <= N; ++i) {
        assert(bCompare(d(i), d5(i), dTol));
    }
}

Here is the call graph for this function:

template<index_type N_SIZE>

void testMatVecGradient2

(

)

Definition at line 722 of file matvectest.cc.

References grad::Alias(), bCompare(), c, callFunc2(), dof, grad::Dot(), grad::MapVectorBase::GLOBAL, grad::sqrt(), and grad::Transpose().

                           {
    LocalDofMap dof;
    Matrix<Gradient<N_SIZE>, 3, 3> A;
    Vector<Gradient<N_SIZE>, 3> b, c;

    for (index_type i = 0; i < 3; ++i) {
        b(i + 1).SetValuePreserve(100*(i + 1));
        b(i + 1).DerivativeResizeReset(&dof, 0L, N_SIZE, MapVectorBase::GLOBAL, 0.);
        c(i + 1).SetValuePreserve(1000*(i + 10));
        c(i + 1).DerivativeResizeReset(&dof, 0L, N_SIZE, MapVectorBase::GLOBAL, 0.);

        for (index_type k = 0; k < N_SIZE; ++k) {
                c(i + 1).SetDerivativeGlobal(k, 3.5 * (k + 1) + 3.2 * (i + 1));
                b(i + 1).SetDerivativeGlobal(k, -17.4 * (k + 1) + 7.5 *(i + 1));
        }

        for (index_type j = 0; j < 3; ++j) {
            A(i + 1, j + 1).SetValuePreserve(100*(i + 1) + j + 1);
            A(i + 1, j + 1).DerivativeResizeReset(&dof, 0L, N_SIZE, MapVectorBase::GLOBAL, 0.);
            for (index_type k = 0; k < N_SIZE; ++k) {
                A(i + 1, j + 1).SetDerivativeGlobal(k, -9.2 * (k + 1.) + 10.5 * (i + 1) + 3.2 * (j + 1) + 7.5);
            }
        }
    }

    Vector<Gradient<N_SIZE>, 3> d, d_C, d_F;

    callFunc2(&dof, A, b, c, d, d_C, d_F);

    Matrix<Gradient<N_SIZE>, 3, 3> A_2 = A * 2.;
    Vector<Gradient<N_SIZE>, 3> b_2 = b * 2.;

    A = Alias(A) * 2.; // Test self assignment
    b = Alias(b) * 2.;

    for (index_type i = 1; i < 3; ++i) {
        assert(bCompare(b(i), b_2(i)));

        for (index_type j = 1; j < 3; ++j) {
                assert(bCompare(A(i, j), A_2(i, j)));
        }
    }

    const doublereal dTol = sqrt(std::numeric_limits<scalar_deriv_type>::epsilon());

    Matrix<Gradient<N_SIZE>, 3, 3> A_3 = (A + Transpose(A) * 0.75) * (2. / b(1));

    A += Transpose(Alias(A)) * 0.75;
    A *= 2.;
    A /= b(1);

    for (index_type i = 1; i < 3; ++i) {
        for (index_type j = 1; j < 3; ++j) {
                assert(bCompare(A(i, j), A_3(i, j), dTol));
        }
    }

    Matrix<Gradient<N_SIZE>, 3, 3> A_4 = A * A(2, 3);
    A *= Alias(A(2, 3));

    for (index_type i = 1; i < 3; ++i) {
        for (index_type j = 1; j < 3; ++j) {
                assert(bCompare(A(i, j), A_4(i, j), dTol));
        }
    }

    Vector<Gradient<N_SIZE>, 3> b_3 = b * b(2);

    b *= Alias(b(2));

    for (index_type i = 1; i <= 3; ++i) {
        assert(bCompare(b(i), b_3(i), dTol));
    }

    Vector<Gradient<N_SIZE>, 3> b_4 = b / sqrt(Dot(b, b));

    b /= sqrt(Dot(Alias(b), b));

    for (index_type i = 1; i <= 3; ++i) {
        assert(bCompare(b(i), b_4(i), dTol));
    }

    Vector<Gradient<N_SIZE>, 3> b_5 = ((b + b * 0.75 / b(1)) * (2. / (A(1, 3) + A(2, 3))) * A(2, 1) + Transpose(A).GetCol(3)) / A(1, 1) / 3.8;

    b += b * 0.75 / Alias(b(1));
    b *= 2.;
    b /= (A(1, 3) + A(2, 3));
    b *= A(2, 1);
    b += Transpose(A).GetCol(3);
    b /= A(1, 1);
    b /= 3.8;

    for (index_type i = 1; i <= 3; ++i) {
        assert(bCompare(b(i), b_5(i), dTol));
    }

    Vector<Gradient<N_SIZE>, 2> b23 = SubVector<2, 3>(b);
    Vector<Gradient<N_SIZE>, 2> b12 = SubVector<1, 2>(b);
    assert(bCompare(b23(1), b(2)));
    assert(bCompare(b23(2), b(3)));
    assert(bCompare(b12(1), b(1)));
    assert(bCompare(b12(2), b(2)));
    Vector<Gradient<N_SIZE>, 2> e = SubVector<2, 3>(b / 2.) + SubVector<1, 2>(b * 3.);
    assert(bCompare(e(1), Gradient<N_SIZE>(b(2) / 2. + b(1) * 3), dTol));
    assert(bCompare(e(2), Gradient<N_SIZE>(b(3) / 2. + b(2) * 3), dTol));

    A = Alias(A) * 0.5 + Transpose(A) * 3.5;
    b = Alias(b) * 2.3 + c * 5.0;
}

Here is the call graph for this function:

void testMatVecProduct

(

)

Definition at line 850 of file matvectest.cc.

References bCompare(), c, grad::Direct(), grad::Dot(), grad::Matrix< T, N_rows, N_cols >::GetCol(), grad::Matrix< T, N_rows, N_cols >::GetRow(), grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::VectorExpression< Expression, N_rows >::iGetNumRows(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), grad::Vector< T, N_rows >::iGetNumRows(), R, grad::Tabular(), and grad::Transpose().

Referenced by main().

                         {
    Matrix<doublereal, 3, 4> A;
    Matrix<doublereal, 3, 3> R;
    Matrix<doublereal, 4, 2> B;

    for (index_type i = 0; i < A.iGetNumRows(); ++i) {
        for (index_type j = 0; j < A.iGetNumCols(); ++j) {
                A(i + 1, j + 1) = 10 * (i + 1) + j + 1;
        }
    }

    for (index_type i = 0; i < R.iGetNumRows(); ++i) {
        for (index_type j = 0; j < R.iGetNumCols(); ++j) {
                R(i + 1, j + 1) = 10 * (i + 1) + j + 1;
        }
    }

    for (index_type i = 0; i < B.iGetNumRows(); ++i) {
        for (index_type j = 0; j < B.iGetNumCols(); ++j) {
                B(i + 1, j + 1) = (10 * (i + 1) + j + 1);
        }
    }

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

    Matrix<doublereal, 4, 3> A_T = Transpose(A);

    std::cout << "A^T=" << std::endl << A_T << std::endl;

    Matrix<doublereal, 4, 3> A_T2 = Transpose(Direct(A));

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
                assert(bCompare(A(i, j), A_T(j, i)));
                assert(bCompare(A(i, j), A_T2(j, i)));
        }
    }

    Matrix<doublereal, 3, 4> A2 = Transpose(Transpose(A));

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
                assert(bCompare(A2(i, j), A(i, j)));
        }
    }

    std::cout << "\nrows of A:" << std::endl;

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        Matrix<doublereal, 3, 4>::RowVectorType r = A.GetRow(i);
        std::cout << i - 1 << ": ";

        for (index_type j = 1; j <= r.iGetNumRows(); ++j) {
                assert(r(j) == A(i, j));
                std::cout << r(j) << " ";
        }

        std::cout << std::endl;
    }

    std::cout << "\ncolumns of A:" << std::endl;

    for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
        Matrix<doublereal, 3, 4>::ColumnVectorType c = A.GetCol(j);
        std::cout << j - 1 << ": ";

        for (index_type i = 1; i <= c.iGetNumRows(); ++i) {
                assert(c(i) == A(i, j));
                std::cout << c(i) << " ";
        }

        std::cout << std::endl;
    }

    Matrix<doublereal, 3, 4> D = -A;
    Matrix<doublereal, 3, 3> E = -(A * Transpose(A));


    for (int i = 1; i <= 3; ++i) {
        for (int j = 1; j <= 3; ++j) {
                assert(bCompare(D(i, j), -A(i, j)));
        }
    }

    assert(E(1, 1) == -630);
    assert(E(1, 2) == -1130);
    assert(E(1, 3) == -1630);
    assert(E(2, 1) == -1130);
    assert(E(2, 2) == -2030);
    assert(E(2, 3) == -2930);
    assert(E(3, 1) == -1630);
    assert(E(3, 2) == -2930);
    assert(E(3, 3) == -4230);

    Vector<doublereal, 4> x;
    Vector<doublereal, 3> b;

    for (index_type i = 0; i < x.iGetNumRows(); ++i) {
        x(i + 1) = 100 * (i + 1);
    }

    for (index_type i = 1; i <= b.iGetNumRows(); ++i) {
        b(i) = Dot(A.GetRow(i), x);
    }

    Vector<Gradient<0>, 3> grad_b = Direct(b);
    Vector<doublereal, 3> b2 = A * x;
    Vector<doublereal, 3> b3 = b + A * x + b;
    Vector<doublereal, 3> b4 = Direct(A) * Direct(x);
    Vector<doublereal, 3> b5 = Direct(A) * x;
    Vector<doublereal, 3> b6 = A * Direct(x);
    Vector<doublereal, 3> c = R * (A * x);
    Vector<doublereal, 4> d = Transpose(A) * b;
    Vector<doublereal, 4> e = Transpose(A) * (A * x);
    Vector<doublereal, 3> f = A * (Transpose(A) * (A * (x + x)));

    for (index_type i = 1; i <= d.iGetNumRows(); ++i) {
        assert(d(i) == e(i));
    }

    Matrix<doublereal, 3, 2> C = Direct(A) * Direct(B);
    Vector<doublereal, 2> g;
    g(1) = 523;
    g(2) = -786;
    Vector<doublereal, 3> h = A * (B * g);
    Vector<doublereal, 3> h2 = (Direct(A) * Direct(B)) * g;
    Vector<doublereal, 3> h3 = (Direct(A) * B) * g;
    Vector<doublereal, 3> h4 = (A * Direct(B)) * g;
    Vector<doublereal, 3> h5 = (A * B) * g;
    Vector<doublereal, 3> h6 = (A * B) * Direct(g);

    for (index_type i = 1; i <= h.iGetNumRows(); ++i) {
        assert(h(i) == h2(i));
        assert(h(i) == h3(i));
        assert(h(i) == h4(i));
        assert(h(i) == h5(i));
        assert(h(i) == h6(i));
    }

    std::cout << "B=" << std::endl << B << std::endl;
    std::cout << "A * B  = C" << std::endl << C << std::endl;

    std::cout << "R=" << std::endl << R << std::endl;
    std::cout << "x = " << std::endl << x << std::endl;
    std::cout << "A * x = b" << std::endl << b << std::endl;
    std::cout << "R * A * x = c" << std::endl << c << std::endl;
    std::cout << "A^T * b = d" << std::endl << d << std::endl;
    std::cout << "A^T * A * x = e" << std::endl << e << std::endl;
    std::cout << "A * A^T * A * 2 * x = f" << std::endl << f << std::endl;
    std::cout << "A * B * g = h" << std::endl << h << std::endl;

    for (index_type i = 1; i <= b2.iGetNumRows(); ++i) {
        assert(b2(i) == b(i));
        assert(b3(i) == 3 * b(i));
        assert(b4(i) == b(i));
        assert(b5(i) == b(i));
        assert(b6(i) == b(i));
    }

    assert(b(1) == 13000);
    assert(b(2) == 23000);
    assert(b(3) == 33000);

    assert(c(1) == 848000);
    assert(c(2) == 1538000);
    assert(c(3) == 2228000);

    assert(d(1) == 1649000.00);
    assert(d(2) == 1718000.00);
    assert(d(3) == 1787000.00);
    assert(d(4) == 1856000.00);

    assert(C(1, 1) == 1.35e3);
    assert(C(2, 1) == 2.39e3);
    assert(C(3, 1) == 3.43e3);
    assert(C(1, 2) == 1.4e3);
    assert(C(2, 2) == 2.48e3);
    assert(C(3, 2) == 3.56e3);

    assert(h(1) == -394350);
    assert(h(2) == -699310);
    assert(h(3) == -1004270);

    Vector<doublereal, 2> b23 = SubVector<2, 3>(b);
    Vector<doublereal, 2> b12 = SubVector<1, 2>(b);
    assert(bCompare(b23(1), b(2)));
    assert(bCompare(b23(2), b(3)));
    assert(bCompare(b12(1), b(1)));
    assert(bCompare(b12(2), b(2)));
    Vector<doublereal, 2> e123 = SubVector<2, 3>(b / 2.) + SubVector<1, 2>(b * 3.);
    assert(bCompare(e123(1), (b(2) / 2. + b(1) * 3)));
    assert(bCompare(e123(2), (b(3) / 2. + b(2) * 3)));

    Matrix<doublereal, 5, 7> F;

    for (index_type i = 1; i <= 5; ++i) {
        for (index_type j = 1; j <= 7; ++j) {
                F(i, j) = i * 10 + j;
        }
    }

    Matrix<doublereal, 2, 3> G = SubMatrix<3, 4, 5, 7>(F);

    for (index_type i = 1; i <= 2; ++i) {
        for (index_type j = 1; j <= 3; ++j) {
                assert(bCompare(G(i, j), F(i + 2, j + 4)));
        }
    }

    Vector<doublereal, 3> G1 = SubMatrix<4, 4, 2, 4>(F).GetRow(1);

    assert(bCompare(G1(1), F(4, 2)));
    assert(bCompare(G1(2), F(4, 3)));
    assert(bCompare(G1(3), F(4, 4)));

    Vector<doublereal, 2> G2 = SubMatrix<1, 2, 5, 7>(F).GetCol(1);

    assert(bCompare(G2(1), F(1, 5)));
    assert(bCompare(G2(2), F(2, 5)));

    Vector<doublereal, 2> G3 = SubMatrix<1, 2, 5, 7>(F).GetCol(2);

    assert(bCompare(G3(1), F(1, 6)));
    assert(bCompare(G3(2), F(2, 6)));


    Vector<doublereal, 2> G4 = SubMatrix<1, 2, 5, 7>(F).GetCol(3);

    assert(bCompare(G4(1), F(1, 7)));
    assert(bCompare(G4(2), F(2, 7)));

    Matrix<doublereal, 3, 4> M = SubMatrix<2, 4, 2, 5>(F) * SubMatrix<2, 5, 4, 7>(F);

    assert(bCompare(M(1, 1), 3716.));
    assert(bCompare(M(1, 2), 3810.));
    assert(bCompare(M(1, 3), 3904.));
    assert(bCompare(M(1, 4), 3998.));
    assert(bCompare(M(2, 1), 5276.));
    assert(bCompare(M(2, 2), 5410.));
    assert(bCompare(M(2, 3), 5544.));
    assert(bCompare(M(2, 4), 5678.));
    assert(bCompare(M(3, 1), 6836.));
    assert(bCompare(M(3, 2), 7010.));
    assert(bCompare(M(3, 3), 7184.));
    assert(bCompare(M(3, 4), 7358.));

    std::cout << "F=\n" << Tabular(F, 5) << std::endl;
    std::cout << "G=\n" << Tabular(G, 5) << std::endl;
    std::cout << "G1=\n" << G1 << std::endl;
    std::cout << "G2=\n" << G2 << std::endl;
    std::cout << "G3=\n" << G3 << std::endl;
    std::cout << "G4=\n" << G4 << std::endl;
    std::cout << "M=\n" << Tabular(M, 5) << std::endl;
}

Here is the call graph for this function:

void testMatVecProductGradient

(

)

Definition at line 1383 of file matvectest.cc.

References bCompare(), grad::Gradient< N_SIZE >::bIsEqual(), c, grad::Cross(), grad::dGetValue(), grad::Direct(), dof, grad::Dot(), grad::exp(), grad::Matrix< T, N_rows, N_cols >::GetCol(), grad::Matrix< T, N_rows, N_cols >::GetRow(), grad::MapVectorBase::GLOBAL, gradVecAssTest::I1, grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), grad::Vector< T, N_rows >::iGetNumRows(), grad::MatCrossVec(), grad::Norm(), testMatVecProductGradient_testData::oct_b, testMatVecProductGradient_testData::oct_c, testMatVecProductGradient_testData::oct_d, testMatVecProductGradient_testData::oct_e, testMatVecProductGradient_testData::oct_f, testMatVecProductGradient_testData::oct_g, testMatVecProductGradient_testData::oct_i, testMatVecProductGradient_testData::oct_j, testMatVecProductGradient_testData::oct_l, testMatVecProductGradient_testData::oct_m, testMatVecProductGradient_testData::oct_norm_g, testMatVecProductGradient_testData::oct_r, testMatVecProductGradient_testData::oct_s, testMatVecProductGradient_testData::oct_t, grad::pow(), grad::Matrix< T, N_rows, N_cols >::Reset(), grad::Vector< T, N_rows >::Reset(), testGradient(), and grad::Transpose().

Referenced by main().

                                 {
    Matrix<Gradient<0>, 3, 4> A;
    LocalDofMap dof;

    for (index_type i = 0; i < A.iGetNumRows(); ++i) {
        for (index_type j = 0; j < A.iGetNumCols(); ++j) {
                A(i + 1, j + 1).SetValue(10 * (i + 1) + j + 1);
                A(i + 1, j + 1).DerivativeResizeReset(&dof, j * A.iGetNumRows() + i + 1, MapVectorBase::GLOBAL, 0.);
                A(i + 1, j + 1).SetDerivativeGlobal(j * A.iGetNumRows() + i + 1, 1.);
        }
    }

    Matrix<Gradient<0>, 3, 4> B1;
    Matrix<doublereal, 3, 4> Ad;

    for (int i = 1; i <= 3; ++i) {
        for (int j = 1; j <= 4; ++j) {
                B1(i, j) = 3.5 * A(i, j);
                Ad(i, j) = A(i, j).dGetValue();
        }
    }

    Matrix<Gradient<0>, 3, 4> C1 = A + B1;
    Matrix<Gradient<0>, 3, 4> D1 = (A + B1) - (C1 + B1);
    Matrix<Gradient<0>, 3, 4> E1 = A + (C1 - B1);
    Matrix<Gradient<0>, 3, 4> F1 = (C1 - B1) - A;
    Matrix<Gradient<0>, 3, 4> G1 = A - B1;
    Matrix<Gradient<0>, 3, 4> H1 = A * 3.5 + B1 / 4.;
    Matrix<Gradient<0>, 3, 4> I1 = (A * 2. - B1 / 5.) * 3.5;
    Matrix<Gradient<0>, 3, 4> J1 = A * B1(1, 1) + I1 * H1(2, 2);
    Matrix<Gradient<0>, 3, 4> K1 = (I1 + J1) * H1(3, 4);

    const doublereal dTol = 10 * std::numeric_limits<scalar_deriv_type>::epsilon();

    for (int i = 1; i <= 3; ++i) {
        for (int j = 1; j <= 4; ++j) {
                assert(bCompare(C1(i, j), Gradient<0>(A(i, j) + B1(i, j)), dTol));
                assert(bCompare(D1(i, j), Gradient<0>(A(i, j) + B1(i, j) - C1(i, j) - B1(i, j)), dTol));
                assert(bCompare(E1(i, j), Gradient<0>(A(i, j) + C1(i, j) - B1(i, j)), dTol));
                assert(bCompare(F1(i, j), Gradient<0>(C1(i, j) - B1(i, j) - A(i, j)), dTol));
                assert(bCompare(G1(i, j), Gradient<0>(A(i, j) - B1(i, j)), dTol));
                assert(bCompare(H1(i, j), Gradient<0>(A(i, j) * 3.5 + B1(i, j) / 4.), dTol));
                assert(bCompare(I1(i, j), Gradient<0>((A(i, j) * 2. - B1(i, j) / 5.) * 3.5), dTol));
                assert(bCompare(J1(i, j), Gradient<0>(A(i, j) * B1(1, 1) + I1(i, j) * H1(2, 2)), dTol));
                assert(bCompare(K1(i, j), Gradient<0>((I1(i, j) + J1(i, j)) * H1(3, 4)), dTol));
        }
    }

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

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

    Matrix<Gradient<0>, 4, 3> A_T = Transpose(A);

    std::cout << "A^T=" << std::endl << A_T << std::endl;

    Matrix<Gradient<0>, 4, 3> A_T2 = Transpose(Direct(A));

    Matrix<Gradient<0>, 3, 4> B = -A;

    Matrix<Gradient<0>, 3, 3> D = -(A * Transpose(A));

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
                assert(bCompare(A(i, j), A_T(j, i)));
                assert(bCompare(A(i, j), A_T2(j, i)));
                assert(bCompare(B(i, j), Gradient<0>(-A(i, j))));
        }
    }

    Matrix<Gradient<0>, 3, 4> A2 = Transpose(Transpose(A));

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
                assert(bCompare(A2(i, j), A(i, j)));
        }
    }

    std::cout << "\nrows of A:" << std::endl;

    for (index_type i = 1; i <= A.iGetNumRows(); ++i) {
        Matrix<Gradient<0>, 3, 4>::RowVectorType r = A.GetRow(i);
        std::cout << i - 1 << ": ";

        for (index_type j = 1; j <= r.iGetNumRows(); ++j) {
                assert(bCompare(r(j), A(i, j), std::numeric_limits<scalar_deriv_type>::epsilon()));
                std::cout << r(j) << " ";
        }

        std::cout << std::endl;
    }

    std::cout << "\ncolumns of A:" << std::endl;

    for (index_type j = 1; j <= A.iGetNumCols(); ++j) {
        Matrix<Gradient<0>, 3, 4>::ColumnVectorType c = A.GetCol(j);
        std::cout << j - 1 << ": ";

        for (index_type i = 1; i <= c.iGetNumRows(); ++i) {
                assert(bCompare(c(i), A(i, j), std::numeric_limits<scalar_deriv_type>::epsilon()));
                std::cout << c(i) << " ";
        }

        std::cout << std::endl;
    }

    Vector<Gradient<0>, 4> x;
    Vector<Gradient<0>, 3> b;
    Vector<doublereal, 4> v;

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
        x(i) = A(1, i) * 10.;
        v(i) = x(i).dGetValue();
    }

    Vector<Gradient<0>, 3> b_1;
    Gradient<0> b1_2;
    b1_2 = A(1, 1) * x(1);
    b1_2 += A(1, 2) * x(2);
    b1_2 += A(1, 3) * x(3);
    b1_2 += A(1, 4) * x(4);

    b_1(1) = A(1, 1) * x(1) + A(1, 2) * x(2) + A(1, 3) * x(3) + A(1, 4) * x(4);

    assert(b1_2.bIsEqual(b_1(1)));

    b_1(2) = A(2, 1) * x(1) + A(2, 2) * x(2) + A(2, 3) * x(3) + A(2, 4) * x(4);
    b_1(3) = A(3, 1) * x(1) + A(3, 2) * x(2) + A(3, 3) * x(3) + A(3, 4) * x(4);


    using namespace testMatVecProductGradient_testData;
    testGradient(oct_b, b_1, dTol);

    for (index_type i = 1; i <= b.iGetNumRows(); ++i) {
        b(i) = Dot(A.GetRow(i), x);
    }

    testGradient(oct_b, b, dTol);

    Vector<Gradient<0>, 3> b2 = A * x;
    Vector<Gradient<0>, 3> b2d = Ad * x;
    Vector<Gradient<0>, 3> b2d1 = -Ad * x - Ad * x;
    Vector<Gradient<0>, 3> b3 = b + A * x + b2;
    Vector<Gradient<0>, 3> b4 = A * (x + x);
    Vector<Gradient<0>, 3> b5 = Direct(A) * Direct(x);
    Vector<Gradient<0>, 3> b6 = Direct(A) * x;
    Vector<Gradient<0>, 3> b7 = A * Direct(x);
    Vector<Gradient<0>, 3> b8 = A * v;
    Vector<Gradient<0>, 3> c = -(A * x);
    Vector<Gradient<0>, 3> d = Cross(c, b) * 5.;
    Vector<Gradient<0>, 3> e = Cross(c + d, b) / 2.;
    Vector<Gradient<0>, 3> f = Cross(e, d + c) - e;
    Vector<Gradient<0>, 3> g = Cross(d + e, f - c) + b;
    Matrix<Gradient<0>, 3, 3> skew_g(MatCrossVec(g));
    Matrix<Gradient<0>, 3, 3> skew_gmf(MatCrossVec(g - f));
    Vector<Gradient<0>, 3> gmf = g - f;
    Gradient<0> norm_Ax = Norm(A * x);
    Gradient<0> z = Norm(A * x) / Norm(Ad * x) + (Norm(A * x) / Norm(Ad * x) - Norm(A * x) / Norm(Ad * x)) * exp(-pow(Norm(x) / 5., 3./2.));
    doublereal zd = Norm(Ad * v) / Norm(Ad * v) + (Norm(Ad * v) / Norm(Ad * v) - Norm(Ad * v) / Norm(Ad * v)) * exp(-pow(Norm(v) / 5., 3./2.));

    for (int i = 1; i <= 3; ++i) {
        for (int j = 1; j <= 3; ++j) {
                if (i == j) {
                        assert(skew_g(i, j).bIsEqual(Gradient<0>()));
                        assert(skew_gmf(i, j).bIsEqual(Gradient<0>()));
                } else {
                        assert(skew_g(i, j).bIsEqual(-skew_g(j, i)));
                        assert(skew_gmf(i, j).bIsEqual(-skew_gmf(j, i)));
                }
        }
    }

    assert(skew_g(1, 2).bIsEqual(-g(3)));
    assert(skew_g(2, 1).bIsEqual(g(3)));
    assert(skew_g(1, 3).bIsEqual(g(2)));
    assert(skew_g(3, 1).bIsEqual(-g(2)));
    assert(skew_g(2, 3).bIsEqual(-g(1)));
    assert(skew_g(3, 2).bIsEqual(g(1)));

    assert(skew_gmf(1, 2).bIsEqual(-gmf(3)));
    assert(skew_gmf(2, 1).bIsEqual(gmf(3)));
    assert(skew_gmf(1, 3).bIsEqual(gmf(2)));
    assert(skew_gmf(3, 1).bIsEqual(-gmf(2)));
    assert(skew_gmf(2, 3).bIsEqual(-gmf(1)));
    assert(skew_gmf(3, 2).bIsEqual(gmf(1)));

    Vector<doublereal, 3> h;
    h(1) = 15.7;
    h(2) = -7.3;
    h(3) = 12.4;

    Vector<Gradient<0>, 3> i = Cross(d, h) * 5.7;
    Vector<Gradient<0>, 3> j = Cross(h, g) / 3.5;
    Vector<doublereal, 3> k = Cross(h, h) + h * 3.;
    Vector<Gradient<0>, 3> l = Cross(h + h, j) / 2.;
    Vector<Gradient<0>, 3> m = Cross(h, j - i) * 5.;
    Vector<doublereal, 3> n = Cross(h + h, h) + h;
    Vector<doublereal, 3> o = Cross(h, h + h) * 2.;
    Vector<doublereal, 3> p = Cross(h * 2.5, h/3.) - h;
    Gradient<0> r = Dot(h, g);
    Gradient<0> s = Dot(g, h);
    Gradient<0> t = Dot(g, g);
    doublereal u = Dot(h, h);
    Gradient<0> r1 = Dot(Direct(h), g);
    Gradient<0> r2 = Dot(h, Direct(g));
    Gradient<0> r3 = Dot(Direct(h), Direct(g));
    Gradient<0> norm_g1 = Norm(g);
    Gradient<0> norm_g2 = Norm(Cross(d + e, f - c) + b);

    assert(r.bIsEqual(r1));
    assert(r.bIsEqual(r2));
    assert(r.bIsEqual(r3));

    std::cout << "x = " << std::endl << x << std::endl;
    std::cout << "b = A * x" << std::endl << b << std::endl;
    std::cout << "c = -A * x" << std::endl << c << std::endl;
    std::cout << "d = cross(c, b) * 5" << std::endl << d << std::endl;
    std::cout << "e = cross(c + d, b) / 2" << std::endl << e << std::endl;
    std::cout << "f = cross(e, d + c) - e" << std::endl << f << std::endl;
    std::cout << "g = cross(d + e, f - c) + b" << std::endl << g << std::endl;
    std::cout << "h = " << std::endl << h << std::endl;
    std::cout << "i = cross(d, h) * 5.7" << std::endl << i << std::endl;
    std::cout << "j = cross(h, g) / 3.5" << std::endl << j << std::endl;
    std::cout << "k = cross(h, h) + h * 3" << std::endl << k << std::endl;
    std::cout << "l = cross(h + h, j) / 2" << std::endl << l << std::endl;
    std::cout << "m = cross(h, j - i) * 5" << std::endl << m << std::endl;
    std::cout << "n = cross(h + h, h) + h" << std::endl << n << std::endl;
    std::cout << "o = cross(h, h + h) * 2" << std::endl << o << std::endl;
    std::cout << "p = cross(h * 2.5, h / 3) - h" << std::endl << p << std::endl;
    std::cout << "r = Dot(h, g) = " << r << std::endl;
    std::cout << "s = Dot(g, h) = " << s << std::endl;
    std::cout << "t = Dot(g, g) = " << t << std::endl;
    std::cout << "u = Dot(h, h) = " << u << std::endl;
    std::cout << "norm_g1 = Norm(g) = " << norm_g1 << std::endl;
    std::cout << "norm_g2 = Norm(Cross(d + e, f - c) + b) = " << norm_g2 << std::endl;
    std::cout << "z=" << z << std::endl;
    std::cout << "zd=" << zd << std::endl;

    for (index_type i = 1; i <= b2.iGetNumRows(); ++i) {
        assert(bCompare(b2d(i).dGetValue(), b2(i).dGetValue(), dTol));
        assert(bCompare(b2d1(i).dGetValue(), -2 * b2(i).dGetValue(), dTol));
        assert(bCompare(b2(i), b(i), dTol));
        assert(bCompare(b3(i), Gradient<0>(3 * b(i)), dTol));
        assert(bCompare(b4(i), Gradient<0>(2 * b(i)), dTol));
        assert(bCompare(b5(i), b(i), dTol));
        assert(bCompare(b6(i), b(i), dTol));
        assert(bCompare(b7(i), b(i), dTol));
        assert(bCompare<scalar_deriv_type>(b8(i).dGetValue(), b(i).dGetValue(), dTol));
        assert(bCompare(c(i), Gradient<0>(-b(i)), dTol));
    }


    assert(bCompare(b(1).dGetValue(), 6300., dTol));
    assert(bCompare(b(2).dGetValue(), 11300., dTol));
    assert(bCompare(b(3).dGetValue(), 16300., dTol));

    assert(bCompare(D(1, 1).dGetValue(), -630., dTol));
    assert(bCompare(D(1, 2).dGetValue(), -1130., dTol));
    assert(bCompare(D(1, 3).dGetValue(), -1630., dTol));
    assert(bCompare(D(2, 1).dGetValue(), -1130., dTol));
    assert(bCompare(D(2, 2).dGetValue(), -2030., dTol));
    assert(bCompare(D(2, 3).dGetValue(), -2930., dTol));
    assert(bCompare(D(3, 1).dGetValue(), -1630., dTol));
    assert(bCompare(D(3, 2).dGetValue(), -2930., dTol));
    assert(bCompare(D(3, 3).dGetValue(), -4230., dTol));

    testGradient(oct_b, b, dTol);
    testGradient(oct_c, c, dTol);
    testGradient(oct_d, d, dTol);
    testGradient(oct_e, e, dTol);
    testGradient(oct_f, f, dTol);
    testGradient(oct_g, g, dTol);
    testGradient(oct_i, i, dTol);
    testGradient(oct_j, j, dTol);
    testGradient(oct_l, l, dTol);
    testGradient(oct_m, m, dTol);
    testGradient(oct_r, r, dTol);
    testGradient(oct_s, s, dTol);
    testGradient(oct_t, t, dTol);
    testGradient(oct_norm_g, norm_g1, dTol);
    testGradient(oct_norm_g, norm_g2, dTol);

    A.Reset();
    x.Reset();
    b.Reset();
    std::cout << "after reset ...\n";
    std::cout << "A=" << A << std::endl;
    std::cout << "x=" << x << std::endl;
    std::cout << "b=" << b << std::endl;
}

Here is the call graph for this function:

template<index_type N_SIZE>

void testMatVecProductGradient2	(	index_type	iNumDeriv,
		int	N
	)

Definition at line 1675 of file matvectest.cc.

References dof, grad::MapVectorBase::GLOBAL, grad::Matrix< T, N_rows, N_cols >::iGetNumCols(), grad::Matrix< T, N_rows, N_cols >::iGetNumRows(), and mbdyn_clock_time().

                                                             {
    Matrix<Gradient<N_SIZE>, 3, 5> A;
    Matrix<Gradient<N_SIZE>, 5, 7> B;
    LocalDofMap dof;

    srand(0);
    
    for (index_type i = 0; i < A.iGetNumRows(); ++i) {
        for (index_type j = 0; j < A.iGetNumCols(); ++j) {
                A(i + 1, j + 1).SetValue(10. * (i + 1) + j + 1 + rand());
                A(i + 1, j + 1).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::GLOBAL, 0.);
            for (int k = 0; k < iNumDeriv; ++k) {
                A(i + 1, j + 1).SetDerivativeGlobal(k, 1000. * (i + 1) + 100. * (j + 1) + k + 1 + rand());
            }
        }
    }
    
    for (index_type i = 0; i < B.iGetNumRows(); ++i) {
        for (index_type j = 0; j < B.iGetNumCols(); ++j) {
                B(i + 1, j + 1).SetValue(1000. * (i + 1) + 100. * (j + 1) + rand());
                B(i + 1, j + 1).DerivativeResizeReset(&dof, 0, iNumDeriv, MapVectorBase::GLOBAL, 0.);
            for (int k = 0; k < iNumDeriv; ++k) {
                B(i + 1, j + 1).SetDerivativeGlobal(k, 10000. * (i + 1) + 1000. * (j + 1) + 10. * (k + 1) + rand());
            }
        }
    }    

    Matrix<Gradient<N_SIZE>, 3, 7> C;
    
    double start = mbdyn_clock_time();
    
    for (int i = 0; i < N; ++i) {
        C = A * B;
    }

    double dt = (mbdyn_clock_time() - start) / N;
    
    std::cerr << "iMaxDerivatives=" << iNumDeriv << std::endl;
    std::cerr << "dt=" << dt << std::endl;
    std::cout << "A=\n" << A << std::endl;
    std::cout << "B=\n" << B << std::endl;
    std::cout << "C=\n" << C << std::endl;
}

Here is the call graph for this function:

void testScalarTypeTraits

(

)

Definition at line 92 of file matvectest.cc.

References grad::Gradient< N_SIZE >::iGetMaxDerivatives(), grad::RangeVector< T, N_SIZE >::iGetMaxSize(), and v0.

Referenced by main().

                            {
        typedef ScalarBinaryExpressionTraits<FuncPlus, doublereal, doublereal, doublereal>::ExpressionType Expr001;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, Gradient<0>, Gradient<0> >::ExpressionType Expr002;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, doublereal, Gradient<0> >::ExpressionType Expr003;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, Gradient<0>, doublereal>::ExpressionType Expr004;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, Expr003, Expr004>::ExpressionType Expr005;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, Expr005, Expr001>::ExpressionType Expr006;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, Expr005, Expr006>::ExpressionType Expr007;
        typedef ScalarBinaryExpressionTraits<FuncPlus, doublereal, Expr001, Expr001>::ExpressionType Expr008;
        typedef ScalarBinaryExpressionTraits<FuncPlus, doublereal, Expr001, Expr008>::ExpressionType Expr009;
        typedef ScalarBinaryExpressionTraits<FuncPlus, doublereal, Expr008, Expr009>::ExpressionType Expr010;
        typedef ScalarBinaryExpressionTraits<FuncPlus, Gradient<0>, Expr010, Expr007>::ExpressionType Expr011;
        typedef BasicScalarType<Expr001>::ScalarType T001;
        typedef BasicScalarType<Expr002>::ScalarType T002;
        typedef BasicScalarType<Expr003>::ScalarType T003;
        typedef BasicScalarType<Expr004>::ScalarType T004;
        typedef BasicScalarType<Expr005>::ScalarType T005;
        typedef BasicScalarType<Expr006>::ScalarType T006;
        typedef BasicScalarType<Expr007>::ScalarType T007;
        typedef BasicScalarType<Expr008>::ScalarType T008;
        typedef BasicScalarType<Expr009>::ScalarType T009;
        typedef BasicScalarType<Expr010>::ScalarType T010;
        typedef BasicScalarType<Expr011>::ScalarType T011;
        typedef CrossTraits<VectorDirectExpr<Vector<Gradient<0>, 3> >, VectorDirectExpr<Vector<Gradient<0>, 3> > >::ExpressionType C001;
        typedef CrossTraits<VectorDirectExpr<Vector<doublereal, 3> >, VectorDirectExpr<Vector<doublereal, 3> > >::ExpressionType C002;
        std::cout << "Expr001=" << typeid(Expr001).name() << std::endl;
        std::cout << "Expr002=" << typeid(Expr002).name() << std::endl;
        std::cout << "Expr003=" << typeid(Expr003).name() << std::endl;
        std::cout << "Expr004=" << typeid(Expr004).name() << std::endl;
        std::cout << "Expr005=" << typeid(Expr005).name() << std::endl;
        std::cout << "Expr006=" << typeid(Expr006).name() << std::endl;
        std::cout << "Expr007=" << typeid(Expr007).name() << std::endl;
        std::cout << "T001=" << typeid(T001).name() << std::endl;
        std::cout << "T002=" << typeid(T002).name() << std::endl;
        std::cout << "T003=" << typeid(T003).name() << std::endl;
        std::cout << "T004=" << typeid(T004).name() << std::endl;
        std::cout << "T005=" << typeid(T005).name() << std::endl;
        std::cout << "T006=" << typeid(T006).name() << std::endl;
        std::cout << "T007=" << typeid(T007).name() << std::endl;
        std::cout << "T008=" << typeid(T008).name() << std::endl;
        std::cout << "T009=" << typeid(T009).name() << std::endl;
        std::cout << "T010=" << typeid(T010).name() << std::endl;
        std::cout << "T011=" << typeid(T011).name() << std::endl;
        assert(typeid(T001) == typeid(doublereal));
        assert(typeid(T002) == typeid(Gradient<0>));
        assert(typeid(T003) == typeid(Gradient<0>));
        assert(typeid(T004) == typeid(Gradient<0>));
        assert(typeid(T005) == typeid(Gradient<0>));
        assert(typeid(T006) == typeid(Gradient<0>));
        assert(typeid(T007) == typeid(Gradient<0>));
        assert(typeid(T008) == typeid(doublereal));
        assert(typeid(T009) == typeid(doublereal));
        assert(typeid(T010) == typeid(doublereal));
        assert(typeid(T011) == typeid(Gradient<0>));

        const index_type N_rows = 3;
        typedef VectorExpression<VectorVectorVectorBinaryExpr<ScalarBinaryOperation<FuncPlus, ScalarTypeTraits<Gradient<3> >::DirectExpressionType, ScalarTypeTraits<Gradient<3> >::DirectExpressionType>, VectorDirectExpr<Vector<Gradient<3>, N_rows> >, VectorDirectExpr<Vector<Gradient<3>, N_rows> > >, N_rows> V001;
        typedef SumTraits<V001, N_rows, 2> Sum003;
        Sum003 s003;
        std::cout << typeid(s003).name() << std::endl;

        std::cout << "sizeof(Matrix<Gradient<12>, 3, 3>())="
                        << sizeof(Matrix<Gradient<12>, 3, 3>) << std::endl;

        std::cout << "sizeof(Vector<Gradient<12> >())="
                        << sizeof(Vector<Gradient<12>, 3>) << std::endl;

        std::cout << "sizeof(Matrix<Gradient<12>, 3, 3>() * Vector<Gradient<12>, 3>())="
                        << sizeof(Matrix<Gradient<12>, 3, 3>() * Vector<Gradient<12>, 3>()) << std::endl;

        RangeVector<doublereal, 3> v3;
        std::cout << "v3.iGetMaxSize()=" << v3.iGetMaxSize() << std::endl;
        RangeVector<doublereal, 0> v0;
        std::cout << "v0.iGetMaxSize()=" << v0.iGetMaxSize() << std::endl;

        std::cout << typeid(C001).name() << std::endl;
        std::cout << typeid(C002).name() << std::endl;

        Gradient<0> g0;
        std::cout << "g0 max derivatives: " << g0.iGetMaxDerivatives() << std::endl;
}

Here is the call graph for this function:

void testSolve

(

)

Definition at line 2843 of file matvectest.cc.

References Mat3x3::LDLSolve(), grad::Norm(), Mat3x3::Solve(), grad::sqrt(), and grad::Tabular().

Referenced by main().

                 {
        const Mat3x3 A(1.32972137393521,      0.61849905020148,     0.709385530146435,
                                   0.61849905020148,     0.290559435134808,     0.344471283014357,
                                   0.709385530146435,     0.344471283014357,     0.752776020323268);

        const Vec3 b(0.815664130323409,
                         0.255816061333836,
                         0.416955203644826);

        const doublereal dTol = sqrt(std::numeric_limits<doublereal>::epsilon());

        const Vec3 x1 = A.Solve(b);
        const Vec3 x2 = A.LDLSolve(b);
        const doublereal f1 = (A * x1 - b).Norm();
        const doublereal f2 = (A * x2 - b).Norm();
        std::cout << "A=" << Tabular(Matrix<doublereal, 3, 3>(A)) << std::endl;
        std::cout << "b=" << b << std::endl;
        std::cout << "x1=" << x1 << std::endl;
        std::cout << "x2=" << x2 << std::endl;
        std::cout << "f1=" << f1 << std::endl;
        std::cout << "f2=" << f2 << std::endl;
        assert(f1 < dTol);
        assert(f2 < dTol);
}

Here is the call graph for this function:

void testSubVecAss

(

)

Definition at line 2705 of file matvectest.cc.

References grad::GradientAssVec< doublereal >::AddItem(), grad::Gradient< N_SIZE >::DerivativeResizeReset(), SparseSubMatrixHandler::dGetCoef(), MySubVectorHandler::dGetCoef(), dof, grad::MapVectorBase::GLOBAL, SparseSubMatrixHandler::iGetColIndex(), SparseSubMatrixHandler::iGetNumRows(), SparseSubMatrixHandler::iGetRowIndex(), MySubVectorHandler::iGetSize(), grad::Gradient< N_SIZE >::SetDerivativeGlobal(), and grad::Gradient< N_SIZE >::SetValue().

Referenced by main().

                     {
        LocalDofMap dof;
        const integer N = 3;
        MySubVectorHandler vh(N);

        GradientAssVec<doublereal> WorkVec(vh);

        for (integer i = 1; i <= N; ++i) {
                WorkVec.AddItem(i, i * 10.);
        }

        SparseSubMatrixHandler mh(4*N);
        GradientAssVec<Gradient<4> > WorkMat(mh);

        for (integer i = 1; i <= N; ++i) {
                Gradient<4> g;
                g.SetValue(i * 10.);
                g.DerivativeResizeReset(&dof, 1, 5, MapVectorBase::GLOBAL, 0.);
                for (integer k = 1; k <= 4; ++k) {
                        g.SetDerivativeGlobal(k, i + k * 0.1);
                }

                WorkMat.AddItem(i, g);
        }

        std::cout << "WorkVec=" << std::endl;
        for (integer i = 1; i <= vh.iGetSize(); ++i) {
                std::cout << " " << vh.dGetCoef(i) << std::endl;
        }

        std::cout << std::endl;

        std::cout << "WorkMat=" << std::endl;

        for (integer i = 1; i <= mh.iGetNumRows(); ++i) {
                std::cout << " " << mh.iGetRowIndex(i)
                                  << " " << mh.iGetColIndex(i)
                                  << " " << mh.dGetCoef(i, 0) << std::endl;
        }

        std::cout << std::endl;
}

Here is the call graph for this function:

void testSubVecAssMatVec

(

)

Definition at line 2748 of file matvectest.cc.

References grad::GradientAssVec< doublereal >::AddItem(), grad::GradientAssVecBase::APPEND, SparseSubMatrixHandler::dGetCoef(), MySubVectorHandler::dGetCoef(), dof, grad::MapVectorBase::GLOBAL, SparseSubMatrixHandler::iGetColIndex(), SparseSubMatrixHandler::iGetNumRows(), SparseSubMatrixHandler::iGetRowIndex(), MySubVectorHandler::iGetRowIndex(), and MySubVectorHandler::iGetSize().

Referenced by main().

                           {
        LocalDofMap dof;
        const integer N = 3;
        MySubVectorHandler vh(2*N);

        Vector<doublereal, 3> v;

        for (integer i = 1; i <= 3; ++i) {
                v(i) = i * 10;
        }

        GradientAssVec<doublereal> WorkVec(vh);

        WorkVec.AddItem(1, v);

        SparseSubMatrixHandler mh(2*4*N);
        GradientAssVec<Gradient<4> > WorkMat(mh);

        Vector<Gradient<4>, 3> g;

        for (integer i = 1; i <= 3; ++i) {
                g(i).SetValue(i * 10.);
                g(i).DerivativeResizeReset(&dof, 1, 5, MapVectorBase::GLOBAL, 0.);
                for (integer k = 1; k <= 4; ++k) {
                        g(i).SetDerivativeGlobal(k, i + k * 0.1);
                }
        }

        WorkMat.AddItem(1, g);

        GradientAssVec<doublereal> WorkVec2(vh, GradientAssVecBase::APPEND);
        Vector<doublereal, 3> v2 = v * 2.;
        WorkVec2.AddItem(4, v2);

        GradientAssVec<Gradient<4> > WorkMat2(mh, GradientAssVecBase::APPEND);
        Vector<Gradient<4>, 3> g2 = g * 2.;
        WorkMat2.AddItem(4, g2);

        std::cout << "v=" << std::endl << v << std::endl;
        std::cout << "WorkVec=" << std::endl;
        for (integer i = 1; i <= vh.iGetSize(); ++i) {
                std::cout << " " << vh.iGetRowIndex(i) << " " << vh.dGetCoef(i) << std::endl;
        }

        std::cout << std::endl;

        std::cout << "g=" << std::endl << g << std::endl;

        std::cout << "WorkMat=" << std::endl;

        for (integer i = 1; i <= mh.iGetNumRows(); ++i) {
                std::cout << " " << mh.iGetRowIndex(i)
                                  << " " << mh.iGetColIndex(i)
                                  << " " << mh.dGetCoef(i, 0) << std::endl;
        }

        std::cout << std::endl;
}

Here is the call graph for this function:

template<index_type iRowCount, index_type iMaxDeriv, typename Function , typename FunctionF77 >

void testVecOp	(	const int	M,
		const int	N,
		Function	f,
		FunctionF77	f77,
		const char *	function
	)

Definition at line 1312 of file matvectest.cc.

References grad::dGetValue(), dof, grad::MapVectorBase::GLOBAL, grad::Vector< T, N_rows >::iGetNumRows(), and mbdyn_clock_time().

                                                                                            {
    Vector<Gradient<iMaxDeriv>, iRowCount> x, y, z;
    LocalDofMap dof;

    for (index_type i = 1; i <= x.iGetNumRows(); ++i) {
                x(i).SetValuePreserve(i * 100);
                x(i).DerivativeResizeReset(&dof, 0, N, MapVectorBase::GLOBAL, 0.);
                y(i).SetValuePreserve(i);
                y(i).DerivativeResizeReset(&dof, 0, N, MapVectorBase::GLOBAL, 0.);

        for (int j = 0; j < N; ++j) {
                        x(i).SetDerivativeGlobal(j, (j + 1));
                        y(i).SetDerivativeGlobal(j, (j + 1));
                }
    }

    double start = mbdyn_clock_time();

    for (int i = 0; i < M; ++i) {
        f(x, y, z);
    }

    const double dt = mbdyn_clock_time() - start;

    std::cerr << "C++: testVecAdd<" << iRowCount << "," << iMaxDeriv << ">(" << M << ", " << N << ", \"" << function << "\"):" << dt << std::endl;

    doublereal xF[iRowCount], yF[iRowCount], zF[iRowCount];
    doublereal *xdF = new doublereal[iRowCount*N];
    doublereal *ydF = new doublereal[iRowCount*N];
    doublereal *zdF = new doublereal[iRowCount*N];

    for (int i = 0; i < iRowCount; ++i) {
        xF[i] = x(i + 1).dGetValue();
        yF[i] = y(i + 1).dGetValue();
        zF[i] = 0.;

        for (int j = 0; j < N; ++j) {
                xdF[i * N + j] = x(i + 1).dGetDerivativeLocal(j);
                ydF[i * N + j] = y(i + 1).dGetDerivativeLocal(j);
                zdF[i * N + j] = 0.;
        }
    }

    start = mbdyn_clock_time();

    for (int i = 0; i < M; ++i) {
        f77(xF, xdF, yF, ydF, zF, zdF, iRowCount, N);
    }

    const double dtF77 = mbdyn_clock_time() - start;

    std::cerr << "F77: testVecAdd<" << iRowCount << "," << iMaxDeriv << ">(" << M << ", " << N << ", \"" << function << "\"):" << dtF77 << std::endl;
    std::cerr << "overhead=" << dt/std::max(std::numeric_limits<doublereal>::epsilon(), dtF77) << std::endl;

    for (int i = 0; i < iRowCount; ++i) {
        assert(xF[i] == x(i + 1).dGetValue());
        assert(yF[i] == y(i + 1).dGetValue());
        assert(zF[i] == z(i + 1).dGetValue());

        for (int j = 0; j < N; ++j) {
                assert(xdF[i * N + j] == x(i + 1).dGetDerivativeLocal(j));
                assert(ydF[i * N + j] == y(i + 1).dGetDerivativeLocal(j));
                assert(zdF[i * N + j] == z(i + 1).dGetDerivativeLocal(j));
        }
    }

    delete [] xdF;
    delete [] ydF;
    delete [] zdF;
}

Here is the call graph for this function:

void tic

(

)

Definition at line 2874 of file matvectest.cc.

References mbdyn_clock_time().

Referenced by callFunc(), func2(), and func2ad().

           {
    dStartTime = mbdyn_clock_time();
}

Here is the call graph for this function:

void tic

(

doublereal &

dTime

)

Definition at line 2870 of file matvectest.cc.

References mbdyn_clock_time().

                            {
        dTime = mbdyn_clock_time();
}

Here is the call graph for this function:

doublereal toc

(

)

Definition at line 2878 of file matvectest.cc.

References dStartTime, and mbdyn_clock_time().

Referenced by callFunc().

                 {
        return mbdyn_clock_time() - dStartTime;
}

Here is the call graph for this function:

Variable Documentation

doublereal dStartTime

Definition at line 2868 of file matvectest.cc.

Referenced by toc().

const integer nbdirsmax = 12

Definition at line 277 of file matvectest.cc.

Referenced by func2ad().

int NLoops = 1

Definition at line 82 of file matvectest.cc.

Referenced by callFunc(), callFunc2(), main(), MatManip_test(), and testMatVec3().

int NLoopsAss = 1

Definition at line 83 of file matvectest.cc.