LapackLinearAlgebra.hpp File Reference

LAPACK-backed fast paths for the dynamic, double-precision dense routines, dispatched by overload. More...

#include <algorithm>
#include <vector>
#include "config.h"
#include "math/CholeskyDecomposition.hpp"
#include "math/DynamicRectMatrix.hpp"
#include "math/LU.hpp"

Go to the source code of this file.

Detailed Description

LAPACK-backed fast paths for the dynamic, double-precision dense routines, dispatched by overload.

SKETCH / PROPOSAL – conventions worked out, not yet build-tested.

The generic templates in LU.hpp and CholeskyDecomposition.hpp are kept as the fallback for small fixed-size matrices (Mat3x3d, Mat6x6d in HydroProp, ...) and for non-RealType element types, where the packing + call overhead of going out to LAPACK would lose to the inlined VTK/Crout code. Here we add non-template overloads* for exactly DynamicRectMatrix<RealType>: a non-template overload is preferred over a function template for an exact type match, so any call invertMatrix(A, AI) / CholeskyDecomposition(A, L) with A,L of that concrete type routes here, while every fixed-size caller keeps the generic template unchanged. No edits to the generic headers.

Include this header instead of (or after) LU.hpp/CholeskyDecomposition.hpp at the call site (e.g. in RPYMobility.cpp) so the overloads are visible. With HAVE_LAPACK undefined it simply forwards to the generics.

ASSUMPTIONS to confirm against the tree:

• DynamicRectMatrix<RealType> stores its data contiguously and exposes a raw pointer (assumed getArray() below); if the accessor differs, only the two a = ...getArray() lines change.
• Storage is row-major (OpenMD convention). Since the matrices we factor here (the RPY mobility M and resistance R) are symmetric, the row/column major mismatch with LAPACK is harmless, and for the general inverse the "invert the transpose, read back transposed" identity makes a plain dgetrf/dgetri on the row-major buffer return the correct inverse.
• RealType is double. For a single-precision build (RealType == float), dispatch the s* entry points instead (a thin if-constexpr or a pair of overloaded extern "C" wrappers); only the BLAS/LAPACK names change.

Definition in file LapackLinearAlgebra.hpp.