520 dense cholesky integration (#72)

run like this:
make test_alp_cholesky_ndebug_alp_reference &&  tests/unit/alp_cholesky_ndebug_alp_reference  -fname datasets/spd5.mtx
or 
make test_alp_cholesky_ndebug_alp_reference &&  tests/unit/alp_cholesky_ndebug_alp_reference  -n 100

The first version reads  dense symmetric real matrix in Matrix Market  Array format,  the second version  generated a random SPD matrix internally. Once input matrix is present the code calculates Cholesky decomposition and checks the result. 


* views fix

* test cholesky

* First numerically correct version

* Add missing structure type traits

* Use utils print functions

* Add preliminary checks for cholesky result

* Fix a bug in print_matrix

* Fix band upper limits as the end is not included

TODO: Check other structures

* file reader enabled

currently works on general input matrices:
make test_alp_cholesky_ndebug_alp_reference &&  tests/unit/alp_cholesky_ndebug_alp_reference  -fname ../datasets/spd5_convertotodense.mtx

But it  should also work for symmetric input matrices
make test_alp_cholesky_ndebug_alp_reference &&  tests/unit/alp_cholesky_ndebug_alp_reference  -fname ../datasets/spd5_convertotodense_symmetric.mtx

* clean prints and DEBUG

* Unit test with check or the solution by
calculating the Frobenius norm of H - L L^t
where H is the input matrix and L is Cholesky solution

H is expcted to be (currently real) symmetric matrix in
Matrix Market Array format

The code is currently generating Upper instead of lower triangular matrix,
which is solved by using the transpose view.

* remove temporary parser in implement it in a separate branch

* review and internal matrix generation

* remove temp matrix files

* spaces

Co-authored-by: Denis Jelovina <denis.jelovina@huawei.com>
Co-authored-by: Vladimir Dimic <vladimir.dimic@huawei.com>

Author: 3 people

include/alp/algorithms/cholesky.hpp

@@ -16,9 +16,14 @@
 
 #include <cmath>
 #include <iostream>
+#include <iomanip>
 
 #include <alp.hpp>
 
+#include "../../../tests/utils/print_alp_containers.hpp"
+
+#define TEMP_DISABLE
+
 namespace alp {
 
 	namespace algorithms {
@@ -42,9 +47,9 @@ namespace alp {
 			typename Ring = Semiring< operators::add< D >, operators::mul< D >, identities::zero, identities::one >,
 			typename Minus = operators::subtract< D >,
 			typename Divide = operators::divide< D > >
-		RC cholesky_lowtr(
-			Matrix< D, structures::LowerTriangular, Dense > &L,
-			const Matrix< D, structures::SymmetricPositiveDefinite, Dense > &H,
+		RC cholesky_uptr(
+			Matrix< D, structures::UpperTriangular, Dense > &L,
+			const Matrix< D, structures::Symmetric, Dense > &H,
 			const Ring & ring = Ring(),
 			const Minus & minus = Minus(),
 			const Divide & divide = Divide() ) {
@@ -55,54 +60,116 @@ namespace alp {
 			const size_t n = nrows( H );
 
 			// Out of place specification of the operation
-			Matrix< D, structures::SymmetricPositiveDefinite, Dense > LL( n, n );
+			Matrix< D, structures::Symmetric, Dense > LL( n, n );
 			rc = set( LL, H );
 
-			for( size_t k = 0; k < n ; ++k ) {
+			if( rc != SUCCESS ) {
+				std::cerr << " set( LL, H ) failed\n";
+				return rc;
+			}
+#ifdef DEBUG
+			print_matrix( " -- LL --  " , LL );
+#endif
 
-				const size_t m = n - k - 1;
+			for( size_t k = 0; k < n ; ++k ) {
+#ifdef DEBUG
+				std::cout << "============ Iteration " << k << " ============" << std::endl;
+#endif
 
-				auto v = get_view( LL , utils::range( k, n) , k );
+				auto a = get_view( LL, k, utils::range( k, n ) );
+#ifdef DEBUG
+				print_vector( " -- a --  " , a );
+#endif
 
 				// L[ k, k ] = alpha = sqrt( LL[ k, k ] )
 				Scalar< D > alpha;
 				rc = eWiseLambda(
-							[ &v, &alpha, &ring ]( const size_t i ) {
-								if ( i == 0 ) {
-									v[ i ] = alpha = std::sqrt( v[ i ] );
-								}
-							},
-							v );
-
+					[ &alpha, &ring ]( const size_t i, D &val ) {
+						if ( i == 0 ) {
+#ifdef TEMP_DISABLE
+							internal::setInitialized( alpha, true );
+							*alpha = std::sqrt( val );
+#else
+							(void)set( alpha, std::sqrt( val ) );
+#endif
+							val = *alpha;
+						}
+					},
+					a
+				);
+
+#ifdef DEBUG
+				std::cout << "alpha " << *alpha << std::endl;
+				if( rc != SUCCESS ) {
+					std::cerr << " eWiseLambda( lambda, view ) (0) failed\n";
+					return rc;
+				}
+#endif
+
+				auto v = get_view( LL, k, utils::range( k + 1, n ) );
+#ifdef DEBUG
+				print_vector( " -- v --  " , v );
+#endif
 				// LL[ k + 1: , k ] = LL[ k + 1: , k ] / alpha
 				rc = eWiseLambda(
-							[ &v, &alpha, &divide ]( const size_t i ) {
-								if ( i > 0 ) {
-									foldl( v[ i ], alpha, divide );
-								}
-							},
-							v );
-
+					[ &alpha, &divide ]( const size_t i, D &val ) {
+						(void)i;
+						internal::foldl( val, *alpha, divide );
+					},
+					v
+				);
+
+#ifdef DEBUG
+				if( rc != SUCCESS ) {
+					std::cerr << " eWiseLambda( lambda, view ) (1) failed\n";
+					return rc;
+				}
+#endif
 
 				// LL[ k+1: , k+1: ] -= v*v^T
 				auto LLprim = get_view( LL, utils::range( k + 1, n ), utils::range( k + 1, n ) );
 
-				vvt = outer( v, ring.getMultiplicativeOperator() );
-				rc = foldl( LLprim, vvt, minus );
-
+				auto vvt = outer( v, ring.getMultiplicativeOperator() );
+#ifdef DEBUG
+				print_vector( " -- v --  " , v );
+				print_matrix( " vvt ", vvt );
+#endif
+
+				rc = alp::eWiseLambda(
+					[ &vvt, &minus, &divide ]( const size_t i, const size_t j, D &val ) {
+						internal::foldl(
+							val,
+							internal::access( vvt, internal::getStorageIndex( vvt, i, j ) ),
+							minus
+						);
+					},
+					LLprim
+				);
+#ifdef DEBUG
+				if( rc != SUCCESS ) {
+					std::cerr << " eWiseLambda( lambda, view ) (2) failed\n";
+					return rc;
+				}
+#endif
 			}
 
 			// Finally collect output into L matrix and return
 			for( size_t k = 0; k < n ; ++k ) {
 
 				// L[ k: , k ] = LL[ k: , k ]
-				auto vL  = get_view( L  , utils::range( k, n) , k );
-				auto vLL = get_view( LL , utils::range( k, n) , k );
+				auto vL  = get_view( L, k, utils::range( k, n )  );
+				auto vLL = get_view( LL, k, utils::range( k, n )  );
 
 				rc = set( vL, vLL );
-
+#ifdef DEBUG
+				if( rc != SUCCESS ) {
+					std::cerr << " set( view, view ) failed\n";
+					return rc;
+				}
+#endif
 			}
 
+			assert( rc == SUCCESS );
 			return rc;
 		}
 

include/alp/structures.hpp

@@ -289,6 +289,13 @@ namespace alp {
 			using inferred_structures = std::tuple< General >;
 		};
 
+		struct Square: BaseStructure {
+
+			typedef std::tuple< OpenInterval > band_intervals;
+
+			using inferred_structures = tuple_cat< std::tuple< Square >, General::inferred_structures >::type;
+		};
+
 		/**
 		 * @brief Static and runtime check to determine if a matrix view of structure TargetStructure
 		 * 		  and index mapping functions (IMFs) \a imf_r and \a imf_c can be defined over \a SourceStructure.
@@ -297,7 +304,7 @@ namespace alp {
 		 * @tparam TargetStructure The underlying structure of the target view.
 		 * @param imf_r            The IMF applied to the rows of the source matrix.
 		 * @param imf_c            The IMF applied to the columns of the source matrix.
-	
+
 		 * @return \a false if the function can determined that the new view may alter underlying assumptions
 		 * 			associated with the source structure \a SourceStructure; \a true otherwise.
 		 */
@@ -337,13 +344,6 @@ namespace alp {
 			};
 		};
 
-		struct Square: BaseStructure {
-
-			typedef std::tuple< OpenInterval > band_intervals;
-
-			using inferred_structures = tuple_cat< std::tuple< Square >, General::inferred_structures >::type;
-		};
-
 
 		template<>
 		struct isInstantiable< General, Square > {
@@ -425,6 +425,22 @@ namespace alp {
 			using inferred_structures = tuple_cat< std::tuple< Symmetric >, Square::inferred_structures >::type;
 		};
 
+		template<>
+		struct isInstantiable< General, Symmetric > {
+			template< typename ImfR, typename ImfC >
+			static bool check( const ImfR &imf_r, const ImfC &imf_c ) {
+				return ( ( imf_r.n == imf_c.n ) && ( imf_r.s == imf_c.s ) && ( imf_r.b == imf_c.b ) );
+			};
+		};
+
+		template<>
+		struct isInstantiable< Symmetric, Symmetric > {
+			template< typename ImfR, typename ImfC >
+			static bool check( const ImfR &imf_r, const ImfC &imf_c ) {
+				return imf_r.isSame(imf_c);
+			};
+		};
+
 		struct SymmetricPositiveDefinite: BaseStructure {
 
 			typedef std::tuple< OpenInterval > band_intervals;
@@ -444,14 +460,14 @@ namespace alp {
 
 		struct LowerTrapezoidal: BaseStructure {
 
-			typedef std::tuple< LeftOpenInterval< 0 > > band_intervals;
+			typedef std::tuple< LeftOpenInterval< 1 > > band_intervals;
 
 			using inferred_structures = tuple_cat< std::tuple< LowerTrapezoidal >, Trapezoidal::inferred_structures >::type;
 		};
 
 		struct LowerTriangular: BaseStructure {
 
-			typedef std::tuple< LeftOpenInterval< 0 > > band_intervals;
+			typedef std::tuple< LeftOpenInterval< 1 > > band_intervals;
 
 			using inferred_structures = tuple_cat< std::tuple< LowerTriangular >, Triangular::inferred_structures, LowerTrapezoidal::inferred_structures >::type;
 		};

include/alp/utils/parser/MatrixFileReaderBase.hpp

@@ -263,6 +263,16 @@ namespace alp {
 					return properties._n;
 				}
 
+				/** Check if matrix is symmetric. */
+				bool isSymmetric() const noexcept {
+					return properties._symmetry == MatrixFileProperties::MMsymmetries::SYMMETRIC;
+				}
+
+				/** Check if matrix is symmetric. */
+				bool isGeneral() const noexcept {
+					return properties._symmetry == MatrixFileProperties::MMsymmetries::GENERAL;
+				}
+
 				/**
 				 * If known, returns the number of nonzeroes contained in the matrix file.
 				 *

tests/unit/CMakeLists.txt

@@ -38,6 +38,10 @@ add_grb_executables( add15m add15m.cpp
 #	BACKENDS alp_reference
 #)
 
+add_grb_executables( alp_cholesky alp_cholesky.cpp
+	BACKENDS alp_reference
+)
+
 add_grb_executables( alp_views alp_views.cpp
 	BACKENDS alp_reference
 )