#include <cassert> #include <random> #include <hpc/matvec/copy.h> #include <hpc/matvec/gematrix.h> #include <hpc/matvec/mm.h> #include <hpc/matvec/print.h> #include <hpc/ulmblas/trlsm.h> #include <hpc/ulmblas/trusm.h> #include <hpc/ulmlapack/getf2.h> #include <hpc/ulmlapack/swap.h> // // Random initializer for general matrices: real and complex valued // template <typename Index, typename T> void randomInit(Index m, Index n, T *A, Index incRowA, Index incColA) { std::random_device random; std::default_random_engine mt(random()); std::uniform_real_distribution<T> uniform(-100,100); for (Index i=0; i<m; ++i) { for (Index j=0; j<n; ++j) { A[i*incRowA+j*incColA] = uniform(mt); } } } int main() { using namespace hpc::matvec; typedef double T; typedef std::size_t Index; const Index m = 600; const Index n = 800; GeMatrix<T> A(m, m); GeMatrix<T> X(m, n); GeMatrix<T> B(m, n); GeMatrix<Index> p(m, 1); randomInit(m, m, A.data, A.incRow, A.incCol); randomInit(m, n, X.data, X.incRow, X.incCol); // // A*X = B // mm(T(1), A, X, T(0), B); // // Factorize A=P*L*U // std::size_t info = hpc::ulmlapack::getf2(A.numRows, A.numCols, A.data, A.incRow, A.incCol, p.data, A.incRow); std::printf("getf2 returned: info = %ld\n", info); // // Solve P*L*U*X = B // hpc::ulmlapack::swap(B.numRows, B.numCols, B.data, B.incRow, B.incCol, Index(0), B.numRows, p.data, p.incRow); hpc::ulmblas::trlsm(m, n, T(1), true, A.data, A.incRow, A.incCol, B.data, B.incRow, B.incCol); hpc::ulmblas::trusm(m, n, T(1), false, A.data, A.incRow, A.incCol, B.data, B.incRow, B.incCol); // // Compute residual // double res = 0; for (Index j=0; j<n; ++j) { for (Index i=0; i<m; ++i) { res += std::abs(B(i,j)-X(i,j)); } } std::printf("res=%e\n", res); } |