// license:GPL-2.0+ // copyright-holders:Couriersud /* * mat_cr.h * * Compressed row format matrices * */ #ifndef MAT_CR_H_ #define MAT_CR_H_ #include #include "../plib/pconfig.h" #include "../plib/palloc.h" template struct mat_cr_t { typedef C index_type; typedef T value_type; C diag[N]; // diagonal index pointer n C ia[N+1]; // row index pointer n + 1 C ja[N*N]; // column index array nz_num, initially (n * n) T A[N*N]; // Matrix elements nz_num, initially (n * n) std::size_t size; std::size_t nz_num; explicit mat_cr_t(const std::size_t n) : size(n) , nz_num(0) { #if 0 #if 0 ia = plib::palloc_array(n + 1); ja = plib::palloc_array(n * n); diag = plib::palloc_array(n); #else diag = plib::palloc_array(n + (n + 1) + n * n); ia = diag + n; ja = ia + (n+1); A = plib::palloc_array(n * n); #endif #endif } ~mat_cr_t() { #if 0 plib::pfree_array(diag); #if 0 plib::pfree_array(ia); plib::pfree_array(ja); #endif plib::pfree_array(A); #endif } void set_scalar(const T scalar) { for (std::size_t i=0, e=nz_num; i LUx = r * * ==> Ux = L⁻¹ r = w * * ==> r = Lw * * This can be solved for w using backwards elimination in L. * * Now Ux = w * * This can be solved for x using backwards elimination in U. * */ for (std::size_t i = 1; ia[i] < nz_num; ++i ) { T tmp = 0.0; const std::size_t j1 = ia[i]; const std::size_t j2 = diag[i]; for (std::size_t j = j1; j < j2; ++j ) tmp += LU[j] * r[ja[j]]; r[i] -= tmp; } // i now is equal to n; for (std::size_t i = size; i-- > 0; ) { T tmp = 0.0; const std::size_t di = diag[i]; const std::size_t j2 = ia[i+1]; for (std::size_t j = di + 1; j < j2; j++ ) tmp += LU[j] * r[ja[j]]; r[i] = (r[i] - tmp) / LU[di]; } } }; #endif /* MAT_CR_H_ */