Content |
High-Level BLAS
Für folgende Funktionen im Namensraum hpc::matvec sollen High-Level BLAS-Wrapper implementiert werden:
-
i = iamax(x) zum Bestimmen des betragsgrößten Elements,
-
swap(x, y) zum Vertauschen von Vektoren,
-
scal(alpha, x) zum Skalieren von Vektoren und
-
r(alpha, x, y, A) für die Rang-1 Operation \(A \leftarrow \alpha x\,y^T\).
Am Beispiel von scal soll gezeigt werden, wie dabei die Funktionalität lediglich an die low-level Implementierung in hpc::ulmblas weitergereicht wird:
#ifndef HPC_MATVEC_SCAL_H #define HPC_MATVEC_SCAL_H 1 #include <cassert> #include <type_traits> #include <hpc/ulmblas/gescal.h> #include <hpc/ulmblas/scal.h> #include <hpc/matvec/isdensevector.h> #include <hpc/matvec/isgematrix.h> namespace hpc { namespace matvec { template <typename Alpha, typename MA> typename std::enable_if<IsGeMatrix<MA>::value, void>::type scal(const Alpha &alpha, MA &A) { typedef typename MA::Index Index; const Index m = A.numRows; const Index n = A.numCols; ulmblas::gescal(m, n, alpha, A.data, A.incRow, A.incCol); } template <typename Alpha, typename VX> typename std::enable_if<IsDenseVector<VX>::value, void>::type scal(const Alpha &alpha, VX &x) { ulmblas::scal(x.length, alpha, x.data, x.inc); } } } // namespace matvec, hpc #endif // HPC_MATVEC_SCAL_H
Aufgabe
Implementiert die restlichen BLAS-Funktionen. Zum Testen könnt Ihr folgendes Programm benutzen:
#include <cassert> #include <random> #include <type_traits> #include <hpc/matvec/densevector.h> #include <hpc/matvec/iamax.h> #include <hpc/matvec/r.h> #include <hpc/matvec/scal.h> #include <hpc/matvec/swap.h> #include <hpc/matvec/print.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); } } } template <typename VX> typename std::enable_if<hpc::matvec::IsDenseVector<VX>::value, void>::type randomInit(VX &x) { typedef typename VX::Index Index; randomInit(x.length, Index(1), x.data, x.inc, Index(1)); } template <typename MA> typename std::enable_if<hpc::matvec::IsGeMatrix<MA>::value, void>::type randomInit(MA &A) { randomInit(A.numRows, A.numCols, A.data, A.incRow, A.incCol); } //------------------------------------------------------------------------------ int main() { using namespace hpc::matvec; typedef double T; typedef std::size_t Index; GeMatrix<T, Index> A(8,10); auto x = A.row(2); auto y = x(3, 7); auto z = y(0, 4, 2); randomInit(A); print(A, "A"); print(x, "x"); print(y, "y"); print(z, "z"); printf("iamax(x) = %ld\n", iamax(x)); printf("iamax(y) = %ld\n", iamax(y)); printf("iamax(z) = %ld\n", iamax(z)); scal(2.5, y); auto x1 = x(0,5); auto x2 = x(5,5); swap(x1, x2); print(A, "A"); print(x, "x"); print(y, "y"); print(z, "z"); auto a10 = A.col(0)(1,7); auto a01 = A.row(0)(1,9); auto A11 = A(1,1,7,9); print(a01, "a01"); print(a10, "a10"); print(A11, "A11"); r(1.5, a10, a01, A11); print(A, "A"); }