Skip to content
Permalink
Browse files

Merge pull request #13376 from lindsayad/densematrix_dualreal_13146

Enable DenseMatrix<DualReal>
  • Loading branch information...
moosebuild committed May 7, 2019
2 parents 34bcba7 + 28552f0 commit e7845808d73bab20786d5721a9ffa66574f1662b
@@ -0,0 +1,17 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "libmesh/dense_matrix.h"

#include "metaphysicl/dualnumberarray.h"

template <>
libMesh::DenseMatrix<DualReal>::DenseMatrix(const unsigned int new_m, const unsigned int new_n);
@@ -0,0 +1,116 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "DualReal.h"
#include "MooseError.h"

#include "libmesh/dense_matrix_base_impl.h"
#include "libmesh/dense_matrix_impl.h"

#include "metaphysicl/dualnumberarray.h"

namespace libMesh
{
template <>
void
DenseMatrix<DualReal>::_multiply_blas(const DenseMatrixBase<DualReal> &, _BLAS_Multiply_Flag)
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_lu_decompose_lapack()
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_svd_lapack(DenseVector<Real> &)
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_svd_lapack(DenseVector<Real> &,
DenseMatrix<Number> &,
DenseMatrix<Number> &)
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_svd_helper(
char, char, std::vector<Real> &, std::vector<Number> &, std::vector<Number> &)
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_svd_solve_lapack(const DenseVector<DualReal> & /*rhs*/,
DenseVector<DualReal> & /*x*/,
Real /*rcond*/) const
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_evd_lapack(DenseVector<DualReal> &,
DenseVector<DualReal> &,
DenseMatrix<DualReal> *,
DenseMatrix<DualReal> *)
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_lu_back_substitute_lapack(const DenseVector<DualReal> &,
DenseVector<DualReal> &)
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
void
DenseMatrix<DualReal>::_matvec_blas(
DualReal, DualReal, DenseVector<DualReal> &, const DenseVector<DualReal> &, bool) const
{
mooseError("No blas/lapack support for type ", demangle(typeid(DualReal).name()));
}

template <>
DenseMatrix<DualReal>::DenseMatrix(const unsigned int new_m, const unsigned int new_n)
: DenseMatrixBase<DualReal>(new_m, new_n),
use_blas_lapack(false),
_val(),
_decomposition_type(NONE)
{
this->resize(new_m, new_n);
}

template class DenseMatrixBase<DualReal>;
template class DenseMatrix<DualReal>;

// Instantiations of these template methods needed for unit test
template void DenseMatrix<DualReal>::vector_mult_add(DenseVector<DualReal> &,
const int,
const DenseVector<DualReal> &) const;
template void DenseMatrix<DualReal>::vector_mult_add(DenseVector<DualReal> &,
const double,
const DenseVector<DualReal> &) const;

template void DenseMatrix<DualReal>::cholesky_solve(const DenseVector<DualReal> & b,
DenseVector<DualReal> & x);
}
@@ -0,0 +1,64 @@
//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "DualReal.h"
#include "DenseMatrix.h"

#include "libmesh/dense_vector.h"

#include "metaphysicl/dualnumberarray.h"

#include "gtest/gtest.h"

using namespace libMesh;

TEST(DenseMatrixDualReal, TryOperations)
{
DenseMatrix<DualReal> a(3, 3);
DenseMatrix<DualReal> b(3, 3);
a.left_multiply(b);
a.left_multiply_transpose(b);
a.right_multiply(b);
a.right_multiply_transpose(b);

DenseVector<DualReal> vec(3);
DenseVector<DualReal> dest(3);

a.vector_mult(dest, vec);
a.vector_mult_transpose(dest, vec);
a.vector_mult_add(dest, 1, vec);
a.outer_product(vec, dest);

DenseMatrix<DualReal> sub(2, 2);
a.get_principal_submatrix(2, 2, sub);

vec(0) = 1;
vec(1) = 2;
vec(2) = 3;

a(0, 0) = 1;
a(1, 1) = 1;
a(2, 2) = 1;

a.lu_solve(vec, dest);
EXPECT_NEAR(dest(0).value(), vec(0).value(), TOLERANCE);
EXPECT_NEAR(dest(1).value(), vec(1).value(), TOLERANCE);
EXPECT_NEAR(dest(2).value(), vec(2).value(), TOLERANCE);

EXPECT_NEAR(1, a.det().value(), TOLERANCE);

b(0, 0) = 1;
b(1, 1) = 1;
b(2, 2) = 1;

b.cholesky_solve(vec, dest);
EXPECT_NEAR(dest(0).value(), vec(0).value(), TOLERANCE);
EXPECT_NEAR(dest(1).value(), vec(1).value(), TOLERANCE);
EXPECT_NEAR(dest(2).value(), vec(2).value(), TOLERANCE);
}

0 comments on commit e784580

Please sign in to comment.
You can’t perform that action at this time.