OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed: Moved the global matrix-vector operations to MatVec.C

Browse Source 
They belong better here, mostly because we want to retain matrix.h
a pure template file (no instantiations there).
Also added operator*(Vector&,Matrix&) assuming the transpose matrix.
This commit is contained in:
Knut Morten Okstad 2017-05-02 20:50:25 +02:00
parent d108607d35
commit 2e368af97a
4 changed files with 118 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
src/LinAlg/MatVec.C
View File

@ -7,7 +7,7 @@
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Index 1-based matrices and vectors for algebraic operations.
//! \brief Global algebraic operations on index 1-based matrices and vectors.
//!
//! The two transformation functions defined in this file are rewritten versions
//! of the Fortran subroutines MATTRA and VECTRA from the SAM library.
@ -24,6 +24,67 @@ int utl::nval_per_line = 6;
double utl::zero_print_tol = 1.0e-6;
// Global algebraic matrix-vector operators.
Vector utl::operator* (const Vector& X, Real c)
{
Vector result(X);
return result *= c;
}
Vector utl::operator+ (const Vector& X, const Vector& Y)
{
Vector result(X);
return result.add(Y);
}
Vector utl::operator- (const Vector& X, const Vector& Y)
{
Vector result(X);
return result.add(Y,Real(-1));
}
Matrix utl::operator* (const Matrix& A, Real c)
{
Matrix B(A);
return B.multiply(c);
}
Vector utl::operator* (const Matrix& A, const Vector& X)
{
Vector Y;
A.multiply(X,Y);
return Y;
}
Vector utl::operator* (const Vector& X, const Matrix& A)
{
Vector Y;
A.multiply(X,Y,true);
return Y;
}
Matrix utl::operator* (const Matrix& A, const Matrix& B)
{
Matrix C;
C.multiply(A,B);
return C;
}
/*!
The following matrix multiplication is performed by this function:
\f[ {\bf A} = {\bf T}^T{\bf A}{\bf T} \f]
where \b A is a full, symmetric matrix, and \b T is an identity matrix
with the nodal sub-matrix \b Tn inserted on the diagonal.
*/
bool utl::transform (Matrix& A, const Matrix& T, size_t K)
{
size_t M = A.rows();
@ -74,6 +135,11 @@ bool utl::transform (Matrix& A, const Matrix& T, size_t K)
}
/*!
The vector \b V is pre-multiplied with the transformation matrix \b T which is
the identity matrix with the nodal sub-matrix \b Tn inserted on the diagonal.
*/
bool utl::transform (Vector& V, const Matrix& T, size_t K, bool transpose)
{
size_t M = V.size();
@ -112,6 +178,13 @@ bool utl::invert (Matrix& A)
return A.inverse() > 0.0;
#ifdef HAS_BLAS
if (sizeof(Real) != sizeof(double))
{
std::cerr <<" *** utl::invert: Available for double precision matrices"
<<" only."<< std::endl;
return false;
}
// Use LAPack/BLAS for larger matrices
int INFO, N = A.rows() > A.cols() ? A.rows() : A.cols();
int* IPIV = new int[N];

52
src/LinAlg/MatVec.h
View File

@ -7,7 +7,7 @@
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Index 1-based matrices and vectors for algebraic operations.
//! \brief Global algebraic operations on index 1-based matrices and vectors.
//!
//==============================================================================
@ -16,6 +16,8 @@
#include "matrix.h"
// Some convenience type definitions:
//! A real-valued vector with algebraic operations
typedef utl::vector<Real> Vector;
//! A real-valued matrix with algebraic operations
@ -31,24 +33,58 @@ typedef std::vector<RealArray> Real2DMat;
typedef std::vector<Real2DMat> Real3DMat;
//! An array of real-valued vectors with algebraic operations
typedef std::vector<Vector> Vectors;
//! An array of real-valued matrices with algebraic operations
typedef std::vector<Matrix> Matrices;
namespace utl
{
//! \brief Multiplication of a vector and a scalar.
//! \return \f$ {\bf Y} = c {\bf X} \f$
Vector operator*(const Vector& X, Real c);
//! \brief Multiplication of a scalar and a vector.
//! \return \f$ {\bf Y} = c {\bf X} \f$
inline Vector operator*(Real c, const Vector& X) { return X*c; }
//! \brief Division of a vector by a scalar.
//! \return \f$ {\bf Y} = \frac{1}{d} {\bf X} \f$
inline Vector operator/(const Vector& X, Real d) { return X*(Real(1)/d); }
//! \brief Addition of two vectors.
//! \return \f$ {\bf Z} = {\bf X} + {\bf Y} \f$
Vector operator+(const Vector& X, const Vector& Y);
//! \brief Subtraction of two vectors.
//! \return \f$ {\bf Z} = {\bf X} - {\bf Y} \f$
Vector operator-(const Vector& X, const Vector& Y);
//! \brief Multiplication of a matrix and a scalar.
//! \return \f$ {\bf B} = c {\bf A} \f$
Matrix operator*(const Matrix& A, Real c);
//! \brief Multiplication of a scalar and a matrix.
//! \return \f$ {\bf B} = c {\bf A} \f$
inline Matrix operator*(Real c, const Matrix& A) { return A*c; }
//! \brief Dot product of two vectors.
//! \return \f$ a = {\bf X}^T {\bf Y} \f$
inline Real operator*(const Vector& X, const Vector& Y) { return X.dot(Y); }
//! \brief Multiplication of a matrix and a vector.
//! \return \f$ {\bf Y} = {\bf A} {\bf X} \f$
Vector operator*(const Matrix& A, const Vector& X);
//! \brief Multiplication of a vector and a matrix.
//! \return \f$ {\bf Y} = {\bf A}^T {\bf X} \f$
Vector operator*(const Vector& X, const Matrix& A);
//! \brief Multiplication of two matrices.
//! \return \f$ {\bf C} = {\bf A} {\bf B} \f$
Matrix operator*(const Matrix& A, const Matrix& B);
//! \brief Congruence transformation of a symmetric matrix.
//! \details The following matrix multiplication is performed:
//! \f[ {\bf A} = {\bf T}^T{\bf A}{\bf T} \f]
//! where \b A is a full, symmetric matrix and \b T is an identity matrix
//! with the nodal sub-matrix \b Tn inserted on the diagonal.
//! \param A The matrix to be transformed
//! \param[in] Tn Nodal transformation matrix
//! \param[in] k Index telling where to insert \b Tn on the diagonal of \b T
bool transform(Matrix& A, const Matrix& Tn, size_t k);
//! \brief Congruence transformation of a vector.
//! \details The vector \b V is pre-multiplied with the transformation matrix
//! \b T which is the identity matrix with the nodal sub-matrix \b Tn
//! inserted on the diagonal.
//! \param V The vector to be transformed
//! \param[in] Tn Nodal transformation matrix
//! \param[in] k Index telling where to insert \b Tn on the diagonal of \b T

62
src/LinAlg/matrix.C
View File

@ -1,62 +0,0 @@
#include "matrix.h"
utl::vector<Real> utl::operator*(const utl::vector<Real>& X, Real c)
{
utl::vector<Real> result(X);
return result *= c;
}
utl::vector<Real> utl::operator*(Real c, const utl::vector<Real>& X)
{
utl::vector<Real> result(X);
return result *= c;
}
Real utl::operator*(const utl::vector<Real>& X, const utl::vector<Real>& Y)
{
return X.dot(Y);
}
utl::vector<Real> utl::operator/(const utl::vector<Real>& X, Real d)
{
utl::vector<Real> result(X);
return result /= d;
}
utl::vector<Real> utl::operator+(const utl::vector<Real>& X, const utl::vector<Real>& Y)
{
utl::vector<Real> result(X);
return result += Y;
}
utl::vector<Real> utl::operator-(const utl::vector<Real>& X, const utl::vector<Real>& Y)
{
utl::vector<Real> result(X);
return result -= Y;
}
utl::matrix<Real> utl::operator*(const utl::matrix<Real>& A, Real c)
{
utl::matrix<Real> B(A);
return B.multiply(c);
}
utl::matrix<Real> utl::operator*(Real c, const utl::matrix<Real>& A)
{
utl::matrix<Real> B(A);
return B.multiply(c);
}
utl::vector<Real> utl::operator*(const utl::matrix<Real>& A, const utl::vector<Real>& X)
{
utl::vector<Real> Y;
A.multiply(X,Y);
return Y;
}
utl::matrix<Real> utl::operator*(const utl::matrix<Real>& A, const utl::matrix<Real>& B)
{
utl::matrix<Real> C;
C.multiply(A, B);
return C;
}

38
src/LinAlg/matrix.h
View File

@ -1626,44 +1626,6 @@ namespace utl //! General utility classes and functions.
//=== Global operators ===================================================
//============================================================================
//! \brief Multiplication of a vector and a scalar.
//! \return \f$ {\bf Y} = c {\bf X} \f$
vector<Real> operator*(const vector<Real>& X, Real c);
//! \brief Multiplication of a scalar and a vector.
//! \return \f$ {\bf Y} = c {\bf X} \f$
vector<Real> operator*(Real c, const vector<Real>& X) ;
//! \brief Division of a vector by a scalar.
//! \return \f$ {\bf Y} = \frac{1}{d} {\bf X} \f$
vector<Real> operator/(const vector<Real>& X, Real d);
//! \brief Addition of two vectors.
//! \return \f$ {\bf Z} = {\bf X} + {\bf Y} \f$
vector<Real> operator+(const vector<Real>& X, const vector<Real>& Y);
//! \brief Subtraction of two vectors.
//! \return \f$ {\bf Z} = {\bf X} - {\bf Y} \f$
vector<Real> operator-(const vector<Real>& X, const vector<Real>& Y);
//! \brief Multiplication of a matrix and a scalar.
//! \return \f$ {\bf B} = c {\bf A} \f$
matrix<Real> operator*(const matrix<Real>& A, Real c);
//! \brief Multiplication of a scalar and a matrix.
//! \return \f$ {\bf B} = c {\bf A} \f$
matrix<Real> operator*(Real c, const matrix<Real>& A);
//! \brief Dot product of two vectors.
//! \return \f$ a = {\bf X^T} {\bf Y} \f$
Real operator*(const vector<Real>& X, const vector<Real>& Y);
//! \brief Multiplication of a matrix and a vector.
//! \return \f$ {\bf Y} = {\bf A} {\bf X} \f$
vector<Real> operator*(const matrix<Real>& A, const vector<Real>& X);
//! \brief Multiplication of two matrices.
//! \return \f$ {\bf C} = {\bf A} {\bf B} \f$
matrix<Real> operator*(const matrix<Real>& A, const matrix<Real>& B);
extern double zero_print_tol; //!< Zero tolerance for printing numbers
extern int nval_per_line; //!< Number of values to print per line