1 // jacobi.cpp: Jacobi iterative method
2 //
3 // Copyright (C) 2017-2021 Stillwater Supercomputing, Inc.
4 // Authors: Theodore Omtzigt, Allan Leal
5 //
6 // This file is part of the universal numbers project, which is released under an MIT Open Source license.
7 #ifdef _MSC_VER
8 #pragma warning(disable : 4514) // unreferenced inline function has been removed
9 #pragma warning(disable : 4710) // 'int sprintf_s(char *const ,const size_t,const char *const ,...)': function not inlined
10 #pragma warning(disable : 4820) // 'sw::universal::value<23>': '3' bytes padding added after data member 'sw::universal::value<23>::_sign'
11 #pragma warning(disable : 5045) // Compiler will insert Spectre mitigation for memory load if /Qspectre switch specified
12 #endif
13
14 // standard library
15 #include <limits>
16 // Configure the posit library with arithmetic exceptions
17 // enable posit arithmetic exceptions
18 #define POSIT_THROW_ARITHMETIC_EXCEPTION 1
19 #include <universal/number/posit/posit.hpp>
20 #include <universal/blas/blas.hpp>
21 #include <universal/blas/generators.hpp>
22 #include <universal/blas/solvers/jacobi.hpp>
23
24
25 // specialized for native floating-point
26 template<typename Scalar,
27 typename = typename std::enable_if<std::is_floating_point<Scalar>::value>::type>
28 Scalar normL1(const sw::universal::blas::vector<float>& v) {
29 float L1Norm{ 0 };
30 for (auto e : v) {
31 L1Norm += std::abs(e);
32 }
33 return L1Norm;
34 }
35
main(int argc,char ** argv)36 int main(int argc, char** argv)
37 try {
38 using namespace sw::universal;
39 using namespace sw::universal::blas;
40
41 // constexpr size_t nbits = 16;
42 // constexpr size_t es = 1;
43 // using Scalar = posit<nbits, es>;
44 using Scalar = float;
45 using Matrix = sw::universal::blas::matrix<Scalar>;
46 using Vector = sw::universal::blas::vector<Scalar>;
47
48 if (argc == 1) std::cout << argv[0] << '\n';
49 int nrOfFailedTestCases = 0;
50
51 // Initialize 'A' 'b' & intial guess 'x' * _
52 Matrix A = {
53 { 5, -2, 3, 0},
54 {-3, 9, 1, -2},
55 { 2, -1, -7, 1},
56 { 4, 3, -5, 7} };
57 Vector b = { -1, 2, 3, 0.5 };
58 Vector x = { 0, 0, 0, 0 };
59
60 std::cout << A << '\n';
61 std::cout << b << '\n';
62 size_t iterations = Jacobi(A, b, x);
63 std::cout << "solution in " << iterations << " iterations\n";
64 std::cout << "solution is " << x << '\n';
65 std::cout << A * x << " = " << b << '\n';
66
67 return (nrOfFailedTestCases > 0 ? EXIT_FAILURE : EXIT_SUCCESS);
68 }
69 catch (char const* msg) {
70 std::cerr << "Caught exception: " << msg << std::endl;
71 return EXIT_FAILURE;
72 }
73 catch (const sw::universal::posit_arithmetic_exception& err) {
74 std::cerr << "Uncaught posit arithmetic exception: " << err.what() << std::endl;
75 return EXIT_FAILURE;
76 }
77 catch (const sw::universal::quire_exception& err) {
78 std::cerr << "Uncaught quire exception: " << err.what() << std::endl;
79 return EXIT_FAILURE;
80 }
81 catch (const sw::universal::posit_internal_exception& err) {
82 std::cerr << "Uncaught posit internal exception: " << err.what() << std::endl;
83 return EXIT_FAILURE;
84 }
85 catch (const std::runtime_error& err) {
86 std::cerr << "Uncaught runtime exception: " << err.what() << std::endl;
87 return EXIT_FAILURE;
88 }
89 catch (...) {
90 std::cerr << "Caught unknown exception" << std::endl;
91 return EXIT_FAILURE;
92 }
93