diff --git a/src/common/KokkosKernels_IOUtils.hpp b/src/common/KokkosKernels_IOUtils.hpp index a0756651ee..39d17f7ae5 100644 --- a/src/common/KokkosKernels_IOUtils.hpp +++ b/src/common/KokkosKernels_IOUtils.hpp @@ -49,7 +49,6 @@ #include #include #include -#include #include #include #ifndef _KOKKOSKERNELSIOUTILS_HPP @@ -57,219 +56,221 @@ #include "Kokkos_ArithTraits.hpp" #include -#include "Kokkos_Random.hpp" #include "KokkosKernels_SimpleUtils.hpp" #include -namespace KokkosKernels { +namespace KokkosKernels{ -namespace Impl { -// Get the interval for Kokkos::fill_random -// For real, interval is (-mag, mag) -// For complex, both real and imaginary parts will have interval (-mag, mag) -template -inline void getRandomBounds(double mag, Scalar &start, Scalar &end) { - start = -mag * Kokkos::ArithTraits::one(); - end = mag * Kokkos::ArithTraits::one(); -} - -template <> -inline void getRandomBounds(double mag, Kokkos::complex &start, - Kokkos::complex &end) { - start = Kokkos::complex(-mag, -mag); - end = Kokkos::complex(mag, mag); -} +namespace Impl{ -template <> -inline void getRandomBounds(double mag, Kokkos::complex &start, - Kokkos::complex &end) { - start = Kokkos::complex(-mag, -mag); - end = Kokkos::complex(mag, mag); -} -// MD: Bases on Christian's sparseMatrix_generate function in test_crsmatrix.cpp -// file. -template -void kk_sparseMatrix_generate(OrdinalType nrows, OrdinalType ncols, - SizeType &nnz, OrdinalType row_size_variance, - OrdinalType bandwidth, ScalarType *&values, - SizeType *&rowPtr, OrdinalType *&colInd) { - rowPtr = new SizeType[nrows + 1]; +//MD: Bases on Christian's sparseMatrix_generate function in test_crsmatrix.cpp file. +template< typename ScalarType , typename OrdinalType, typename SizeType> +void kk_sparseMatrix_generate( + OrdinalType nrows, + OrdinalType ncols, + SizeType &nnz, + OrdinalType row_size_variance, + OrdinalType bandwidth, + ScalarType* &values, + SizeType* &rowPtr, + OrdinalType* &colInd) +{ + rowPtr = new SizeType[nrows+1]; - OrdinalType elements_per_row = nrows ? nnz / nrows : 0; + OrdinalType elements_per_row = nrows ? nnz/nrows : 0; srand(13721); rowPtr[0] = 0; - for (int row = 0; row < nrows; row++) { - int varianz = (1.0 * rand() / RAND_MAX - 0.5) * row_size_variance; + for(int row=0;row 0.66 * ncols) numRowEntries = 0.66 * ncols; - rowPtr[row + 1] = rowPtr[row] + numRowEntries; - } - nnz = rowPtr[nrows]; + if(numRowEntries < 0) + numRowEntries = 0; + rowPtr[row+1] = rowPtr[row] + numRowEntries; + } + nnz = rowPtr[nrows]; values = new ScalarType[nnz]; colInd = new OrdinalType[nnz]; - for (OrdinalType row = 0; row < nrows; row++) { - for (SizeType k = rowPtr[row]; k < rowPtr[row + 1]; ++k) { - while (true) { - OrdinalType pos = (1.0 * rand() / RAND_MAX - 0.5) * bandwidth + row; - while (pos < 0) pos += ncols; - while (pos >= ncols) pos -= ncols; + for(OrdinalType row=0;row=ncols) pos-=ncols; bool is_already_in_the_row = false; - for (SizeType j = rowPtr[row]; j < k; j++) { - if (colInd[j] == pos) { + for(SizeType j = rowPtr[row] ; j valuesView(values, nnz); - ScalarType randStart, randEnd; - getRandomBounds(50.0, randStart, randEnd); - Kokkos::Random_XorShift64_Pool pool(13718); - Kokkos::fill_random(valuesView, pool, randStart, randEnd); } -template +template< typename ScalarType , typename OrdinalType, typename SizeType> void kk_sparseMatrix_generate_lower_upper_triangle( - char uplo, OrdinalType nrows, OrdinalType ncols, SizeType &nnz, - OrdinalType /*row_size_variance*/, OrdinalType /*bandwidth*/, - ScalarType *&values, SizeType *&rowPtr, OrdinalType *&colInd) { - rowPtr = new SizeType[nrows + 1]; - - // OrdinalType elements_per_row = nnz/nrows; + char uplo, + OrdinalType nrows, + OrdinalType ncols, + SizeType &nnz, + OrdinalType /*row_size_variance*/, + OrdinalType /*bandwidth*/, + ScalarType* &values, + SizeType* &rowPtr, + OrdinalType* &colInd) +{ + rowPtr = new SizeType[nrows+1]; + + //OrdinalType elements_per_row = nnz/nrows; srand(13721); rowPtr[0] = 0; - for (int row = 0; row < nrows; row++) { - if (uplo == 'L') - rowPtr[row + 1] = rowPtr[row] + row + 1; + for(int row=0;row +template< typename ScalarType , typename OrdinalType, typename SizeType> void kk_diagonally_dominant_sparseMatrix_generate( - OrdinalType nrows, OrdinalType ncols, SizeType &nnz, - OrdinalType row_size_variance, OrdinalType bandwidth, ScalarType *&values, - SizeType *&rowPtr, OrdinalType *&colInd, - ScalarType diagDominance = 10 * Kokkos::ArithTraits::one()) { - rowPtr = new SizeType[nrows + 1]; - - OrdinalType elements_per_row = nnz / nrows; + OrdinalType nrows, + OrdinalType ncols, + SizeType &nnz, + OrdinalType row_size_variance, + OrdinalType bandwidth, + ScalarType* &values, + SizeType* &rowPtr, + OrdinalType* &colInd, + ScalarType diagDominance = 10 * Kokkos::ArithTraits::one()) +{ + rowPtr = new SizeType[nrows+1]; + + OrdinalType elements_per_row = nnz/nrows; srand(13721); rowPtr[0] = 0; - for (int row = 0; row < nrows; row++) { - int varianz = (1.0 * rand() / RAND_MAX - 0.5) * row_size_variance; - if (varianz < 1) varianz = 1; - if (varianz > 0.75 * ncols) varianz = 0.75 * ncols; - rowPtr[row + 1] = rowPtr[row] + elements_per_row + varianz; - if (rowPtr[row + 1] <= rowPtr[row]) // This makes sure that there is - rowPtr[row + 1] = rowPtr[row] + 1; // at least one nonzero in the row - } - nnz = rowPtr[nrows]; + for(int row=0;row entriesInRow; - // We always add the diagonal entry (after this loop) - entriesInRow.insert(row); - for (SizeType k = rowPtr[row]; k < rowPtr[row + 1] - 1; k++) { - while (true) { - OrdinalType pos = (1.0 * rand() / RAND_MAX - 0.5) * bandwidth + row; - while (pos < 0) pos += ncols; - while (pos >= ncols) pos -= ncols; - - if (entriesInRow.find(pos) == entriesInRow.end()) { - entriesInRow.insert(pos); - colInd[k] = pos; - values[k] = 100.0 * rand() / RAND_MAX - 50.0; - total_values += - Kokkos::Details::ArithTraits::abs(values[k]); + for(SizeType k=rowPtr[row]; k= ncols) + pos -= ncols; + + bool is_already_in_the_row = false; + if(pos == row) + is_already_in_the_row = true; + else + { + for(SizeType j = rowPtr[row] ; j::abs(values[k]); break; } } } - colInd[rowPtr[row + 1] - 1] = row; - values[rowPtr[row + 1] - 1] = total_values * diagDominance; + colInd[rowPtr[row+1] - 1]= row; + values[rowPtr[row+1] - 1] = total_values * diagDominance; } } + // This function creates a diagonal sparse matrix for testing matrix operations. // The elements on the diagonal are 1, 2, ..., n-1, n. // If "invert" is true, it will return the inverse of the above diagonal matrix. template crsMat_t kk_generate_diag_matrix(typename crsMat_t::const_ordinal_type n, - const bool invert = false) { + const bool invert = false){ typedef typename crsMat_t::ordinal_type ot; typedef typename crsMat_t::StaticCrsGraphType graph_t; typedef typename graph_t::row_map_type::non_const_type row_map_view_t; - typedef typename graph_t::entries_type::non_const_type cols_view_t; + typedef typename graph_t::entries_type::non_const_type cols_view_t; typedef typename crsMat_t::values_type::non_const_type values_view_t; typedef typename row_map_view_t::non_const_value_type size_type; typedef typename cols_view_t::non_const_value_type lno_t; typedef typename values_view_t::non_const_value_type scalar_t; - row_map_view_t rowmap_view("rowmap_view", n + 1); + row_map_view_t rowmap_view("rowmap_view", n+1); cols_view_t columns_view("colsmap_view", n); values_view_t values_view("values_view", n); { - typename row_map_view_t::HostMirror hr = - Kokkos::create_mirror_view(rowmap_view); - typename cols_view_t::HostMirror hc = - Kokkos::create_mirror_view(columns_view); - typename values_view_t::HostMirror hv = - Kokkos::create_mirror_view(values_view); - - for (lno_t i = 0; i <= n; ++i) { + typename row_map_view_t::HostMirror hr = Kokkos::create_mirror_view (rowmap_view); + typename cols_view_t::HostMirror hc = Kokkos::create_mirror_view (columns_view); + typename values_view_t::HostMirror hv = Kokkos::create_mirror_view (values_view); + + for (lno_t i = 0; i <= n; ++i){ hr(i) = size_type(i); } - for (ot i = 0; i < n; ++i) { + for (ot i = 0; i < n; ++i){ hc(i) = lno_t(i); - if (invert) { - hv(i) = scalar_t(1.0) / (scalar_t(i + 1)); - } else { + if(invert){ + hv(i) = scalar_t(1.0)/(scalar_t(i + 1)); + } + else{ hv(i) = scalar_t(i + 1); } } - Kokkos::deep_copy(rowmap_view, hr); - Kokkos::deep_copy(columns_view, hc); - Kokkos::deep_copy(values_view, hv); + Kokkos::deep_copy (rowmap_view , hr); + Kokkos::deep_copy (columns_view , hc); + Kokkos::deep_copy (values_view , hv); } - graph_t static_graph(columns_view, rowmap_view); + graph_t static_graph (columns_view, rowmap_view); crsMat_t crsmat("CrsMatrix", n, values_view, static_graph); return crsmat; } @@ -282,109 +283,106 @@ crsMat_t kk_generate_diagonally_dominant_sparse_matrix( typename crsMat_t::const_ordinal_type row_size_variance, typename crsMat_t::const_ordinal_type bandwidth, typename crsMat_t::const_value_type diagDominance = - 10 * Kokkos::ArithTraits::one()) { + 10 * Kokkos::ArithTraits::one()) +{ typedef typename crsMat_t::StaticCrsGraphType graph_t; typedef typename graph_t::row_map_type::non_const_type row_map_view_t; - typedef typename graph_t::entries_type::non_const_type cols_view_t; + typedef typename graph_t::entries_type::non_const_type cols_view_t; typedef typename crsMat_t::values_type::non_const_type values_view_t; + typedef typename row_map_view_t::non_const_value_type size_type; typedef typename cols_view_t::non_const_value_type lno_t; typedef typename values_view_t::non_const_value_type scalar_t; lno_t *adj; - size_type *xadj; //, nnzA; + size_type *xadj;//, nnzA; scalar_t *values; kk_diagonally_dominant_sparseMatrix_generate( - nrows, ncols, nnz, row_size_variance, bandwidth, values, xadj, adj, - diagDominance); + nrows, ncols, nnz, row_size_variance, bandwidth, + values, xadj, adj, diagDominance); - row_map_view_t rowmap_view("rowmap_view", nrows + 1); + row_map_view_t rowmap_view("rowmap_view", nrows+1); cols_view_t columns_view("colsmap_view", nnz); values_view_t values_view("values_view", nnz); { - typename row_map_view_t::HostMirror hr = - Kokkos::create_mirror_view(rowmap_view); - typename cols_view_t::HostMirror hc = - Kokkos::create_mirror_view(columns_view); - typename values_view_t::HostMirror hv = - Kokkos::create_mirror_view(values_view); - - for (lno_t i = 0; i <= nrows; ++i) { + typename row_map_view_t::HostMirror hr = Kokkos::create_mirror_view (rowmap_view); + typename cols_view_t::HostMirror hc = Kokkos::create_mirror_view (columns_view); + typename values_view_t::HostMirror hv = Kokkos::create_mirror_view (values_view); + + for (lno_t i = 0; i <= nrows; ++i){ hr(i) = xadj[i]; } - for (size_type i = 0; i < nnz; ++i) { + for (size_type i = 0; i < nnz; ++i){ hc(i) = adj[i]; hv(i) = values[i]; } - Kokkos::deep_copy(rowmap_view, hr); - Kokkos::deep_copy(columns_view, hc); - Kokkos::deep_copy(values_view, hv); + Kokkos::deep_copy (rowmap_view , hr); + Kokkos::deep_copy (columns_view , hc); + Kokkos::deep_copy (values_view , hv); } - graph_t static_graph(columns_view, rowmap_view); + graph_t static_graph (columns_view, rowmap_view); crsMat_t crsmat("CrsMatrix", ncols, values_view, static_graph); - delete[] xadj; - delete[] adj; - delete[] values; + delete [] xadj; delete [] adj; delete [] values; return crsmat; } template crsMat_t kk_generate_triangular_sparse_matrix( - char uplo, typename crsMat_t::const_ordinal_type nrows, + char uplo, + typename crsMat_t::const_ordinal_type nrows, typename crsMat_t::const_ordinal_type ncols, typename crsMat_t::non_const_size_type &nnz, typename crsMat_t::const_ordinal_type row_size_variance, - typename crsMat_t::const_ordinal_type bandwidth) { + typename crsMat_t::const_ordinal_type bandwidth){ + typedef typename crsMat_t::StaticCrsGraphType graph_t; typedef typename graph_t::row_map_type::non_const_type row_map_view_t; - typedef typename graph_t::entries_type::non_const_type cols_view_t; + typedef typename graph_t::entries_type::non_const_type cols_view_t; typedef typename crsMat_t::values_type::non_const_type values_view_t; + typedef typename row_map_view_t::non_const_value_type size_type; typedef typename cols_view_t::non_const_value_type lno_t; typedef typename values_view_t::non_const_value_type scalar_t; lno_t *adj; - size_type *xadj; //, nnzA; + size_type *xadj;//, nnzA; scalar_t *values; kk_sparseMatrix_generate_lower_upper_triangle( - uplo, nrows, ncols, nnz, row_size_variance, bandwidth, values, xadj, adj); + uplo, + nrows, ncols, nnz, row_size_variance, bandwidth, + values, xadj, adj); - row_map_view_t rowmap_view("rowmap_view", nrows + 1); + row_map_view_t rowmap_view("rowmap_view", nrows+1); cols_view_t columns_view("colsmap_view", nnz); values_view_t values_view("values_view", nnz); { - typename row_map_view_t::HostMirror hr = - Kokkos::create_mirror_view(rowmap_view); - typename cols_view_t::HostMirror hc = - Kokkos::create_mirror_view(columns_view); - typename values_view_t::HostMirror hv = - Kokkos::create_mirror_view(values_view); - - for (lno_t i = 0; i <= nrows; ++i) { + typename row_map_view_t::HostMirror hr = Kokkos::create_mirror_view (rowmap_view); + typename cols_view_t::HostMirror hc = Kokkos::create_mirror_view (columns_view); + typename values_view_t::HostMirror hv = Kokkos::create_mirror_view (values_view); + + for (lno_t i = 0; i <= nrows; ++i){ hr(i) = xadj[i]; } - for (size_type i = 0; i < nnz; ++i) { + for (size_type i = 0; i < nnz; ++i){ hc(i) = adj[i]; hv(i) = values[i]; } - Kokkos::deep_copy(rowmap_view, hr); - Kokkos::deep_copy(columns_view, hc); - Kokkos::deep_copy(values_view, hv); + Kokkos::deep_copy (rowmap_view , hr); + Kokkos::deep_copy (columns_view , hc); + Kokkos::deep_copy (values_view , hv); Kokkos::fence(); } - graph_t static_graph(columns_view, rowmap_view); + graph_t static_graph (columns_view, rowmap_view); crsMat_t crsmat("CrsMatrix", ncols, values_view, static_graph); - delete[] xadj; - delete[] adj; - delete[] values; + delete [] xadj; delete [] adj; delete [] values; return crsmat; } @@ -394,81 +392,84 @@ crsMat_t kk_generate_sparse_matrix( typename crsMat_t::const_ordinal_type ncols, typename crsMat_t::non_const_size_type &nnz, typename crsMat_t::const_ordinal_type row_size_variance, - typename crsMat_t::const_ordinal_type bandwidth) { + typename crsMat_t::const_ordinal_type bandwidth){ + typedef typename crsMat_t::StaticCrsGraphType graph_t; typedef typename graph_t::row_map_type::non_const_type row_map_view_t; - typedef typename graph_t::entries_type::non_const_type cols_view_t; + typedef typename graph_t::entries_type::non_const_type cols_view_t; typedef typename crsMat_t::values_type::non_const_type values_view_t; + typedef typename row_map_view_t::non_const_value_type size_type; typedef typename cols_view_t::non_const_value_type lno_t; typedef typename values_view_t::non_const_value_type scalar_t; lno_t *adj; - size_type *xadj; //, nnzA; + size_type *xadj;//, nnzA; scalar_t *values; kk_sparseMatrix_generate( - nrows, ncols, nnz, row_size_variance, bandwidth, values, xadj, adj); + nrows, ncols, nnz, row_size_variance, bandwidth, + values, xadj, adj); - row_map_view_t rowmap_view("rowmap_view", nrows + 1); + row_map_view_t rowmap_view("rowmap_view", nrows+1); cols_view_t columns_view("colsmap_view", nnz); values_view_t values_view("values_view", nnz); { - typename row_map_view_t::HostMirror hr = - Kokkos::create_mirror_view(rowmap_view); - typename cols_view_t::HostMirror hc = - Kokkos::create_mirror_view(columns_view); - typename values_view_t::HostMirror hv = - Kokkos::create_mirror_view(values_view); - - for (lno_t i = 0; i <= nrows; ++i) { + typename row_map_view_t::HostMirror hr = Kokkos::create_mirror_view (rowmap_view); + typename cols_view_t::HostMirror hc = Kokkos::create_mirror_view (columns_view); + typename values_view_t::HostMirror hv = Kokkos::create_mirror_view (values_view); + + for (lno_t i = 0; i <= nrows; ++i){ hr(i) = xadj[i]; } - for (size_type i = 0; i < nnz; ++i) { + for (size_type i = 0; i < nnz; ++i){ hc(i) = adj[i]; hv(i) = values[i]; } - Kokkos::deep_copy(rowmap_view, hr); - Kokkos::deep_copy(columns_view, hc); - Kokkos::deep_copy(values_view, hv); + Kokkos::deep_copy (rowmap_view , hr); + Kokkos::deep_copy (columns_view , hc); + Kokkos::deep_copy (values_view , hv); } - graph_t static_graph(columns_view, rowmap_view); + graph_t static_graph (columns_view, rowmap_view); crsMat_t crsmat("CrsMatrix", ncols, values_view, static_graph); - delete[] xadj; - delete[] adj; - delete[] values; + delete [] xadj; delete [] adj; delete [] values; return crsmat; } -// TODO: need to fix the size_type. All over the reading inputs are lno_t. + + + +//TODO: need to fix the size_type. All over the reading inputs are lno_t. template -void md_malloc(stype **arr, size_t n, std::string /*alloc_str*/ = "") { +void md_malloc(stype **arr, size_t n, std::string /*alloc_str*/ = ""){ *arr = new stype[n]; - if (*arr == NULL) { - throw std::runtime_error("Memory Allocation Problem\n"); + if (*arr == NULL){ + throw std::runtime_error ("Memory Allocation Problem\n"); } } template -struct Edge { +struct Edge{ idx src; idx dst; wt ew; - bool operator<(const Edge &a) const { - // return !((this->src < a.src) || (this->src == a.src && this->dst < - // a.dst)); + bool operator<(const Edge & a) const + { + //return !((this->src < a.src) || (this->src == a.src && this->dst < a.dst)); return (this->src < a.src) || (this->src == a.src && this->dst < a.dst); } }; + //////////////////////////////////////////////////////////////////////////////// // From MTGL //////////////////////////////////////////////////////////////////////////////// -inline size_t kk_get_file_size(const char *file) { +inline size_t kk_get_file_size(const char* file) +{ // struct stat stat_buf; #ifdef _WIN32 @@ -483,85 +484,111 @@ inline size_t kk_get_file_size(const char *file) { } template -void buildEdgeListFromBinSrcTarg_undirected(const char *fnameSrc, - const char *fnameTarg, - size_t &numEdges, lno_t **srcs, - lno_t **dst) { - size_t srcFileSize = kk_get_file_size(fnameSrc); - size_t trgFileSize = kk_get_file_size(fnameTarg); - // test these values +void buildEdgeListFromBinSrcTarg_undirected( + const char *fnameSrc, const char*fnameTarg, + size_t &numEdges, + lno_t **srcs, lno_t **dst){ + size_t srcFileSize = kk_get_file_size(fnameSrc); + size_t trgFileSize = kk_get_file_size(fnameTarg); + // test these values - size_t srcSize = srcFileSize / sizeof(lno_t); - size_t trgSize = trgFileSize / sizeof(lno_t); - if (srcSize != trgSize) { - throw std::runtime_error("Src and Target file needs to be the same size"); - } - // Assumption that each edge is listed once - numEdges = srcSize; + size_t srcSize = srcFileSize / sizeof(lno_t); + size_t trgSize = trgFileSize / sizeof(lno_t); + if (srcSize != trgSize) + { + throw std::runtime_error ("Src and Target file needs to be the same size"); + } + //Assumption that each edge is listed once + numEdges = srcSize; - md_malloc(srcs, numEdges); - md_malloc(dst, numEdges); + md_malloc(srcs, numEdges); + md_malloc(dst, numEdges); - //////////////////////////////////////////////////////// - // Read source data into buffer - //////////////////////////////////////////////////////// + //////////////////////////////////////////////////////// + // Read source data into buffer + //////////////////////////////////////////////////////// - std::ifstream myFile(fnameSrc, std::ios::in | std::ios::binary); + std::ifstream myFile (fnameSrc, std::ios::in | std::ios::binary); - myFile.read((char *)*srcs, sizeof(lno_t) * (numEdges)); + myFile.read((char *) *srcs, sizeof(lno_t) * (numEdges)); - myFile.close(); + myFile.close(); - std::ifstream myFile2(fnameTarg, std::ios::in | std::ios::binary); + std::ifstream myFile2 (fnameTarg, std::ios::in | std::ios::binary); - myFile2.read((char *)*dst, sizeof(lno_t) * (numEdges)); + myFile2.read((char *) *dst, sizeof(lno_t) * (numEdges)); + + myFile2.close(); + // - myFile2.close(); - // } + + template -inline void kk_write_1Dview_to_file(idx_array_type view, const char *filename) { +inline void kk_write_1Dview_to_file(idx_array_type view, const char *filename){ + typedef typename idx_array_type::HostMirror host_type; - // typedef typename idx_array_type::size_type idx; - host_type host_view = Kokkos::create_mirror_view(view); - Kokkos::deep_copy(host_view, view); + //typedef typename idx_array_type::size_type idx; + host_type host_view = Kokkos::create_mirror_view (view); + Kokkos::deep_copy (host_view , view); Kokkos::fence(); - std::ofstream myFile(filename, std::ios::out); - for (size_t i = 0; i < view.extent(0); ++i) { - myFile << host_view(i) << std::endl; + std::ofstream myFile (filename, std::ios::out ); + for (size_t i = 0; i < view.extent(0); ++i){ + myFile << host_view(i) << std::endl; } myFile.close(); } template -inline void kk_read_1Dview_from_file(idx_array_type &view, - const char *filename) { +inline void kk_read_1Dview_from_file(idx_array_type &view, const char *filename){ + typedef typename idx_array_type::HostMirror host_type; - // typedef typename idx_array_type::size_type idx; - host_type host_view = Kokkos::create_mirror_view(view); - std::ifstream myFile(filename, std::ios::in); + //typedef typename idx_array_type::size_type idx; + host_type host_view = Kokkos::create_mirror_view (view); + std::ifstream myFile (filename, std::ios::in ); - for (size_t i = 0; i < view.extent(0); ++i) { - myFile >> host_view(i); + for (size_t i = 0; i < view.extent(0); ++i){ + myFile >> host_view(i); } myFile.close(); - Kokkos::deep_copy(view, host_view); + Kokkos::deep_copy (view, host_view); Kokkos::fence(); } + +template +inline void kk_read_MM_vector(idx_array_type &view, const char *filename){ + + typedef typename idx_array_type::HostMirror host_type; + //typedef typename idx_array_type::size_type idx; + host_type host_view = Kokkos::create_mirror_view (view); + std::ifstream myFile (filename, std::ios::in ); + + std::string line; + std::getline(myFile, line); // have line 1 + std::getline(myFile, line); + + for (size_t i = 0; i < view.extent(0); ++i){ + myFile >> host_view(i); + } + myFile.close(); + Kokkos::deep_copy (view, host_view); + Kokkos::fence(); +} + + + template -void convert_crs_to_lower_triangle_edge_list(idx nv, idx *xadj, idx *adj, - idx *lower_triangle_srcs, - idx *lower_triangle_dests) { +void convert_crs_to_lower_triangle_edge_list(idx nv, idx *xadj, idx *adj, idx *lower_triangle_srcs, idx *lower_triangle_dests){ idx ind = 0; - for (idx i = 0; i < nv; ++i) { + for (idx i = 0; i < nv; ++i){ idx xb = xadj[i]; - idx xe = xadj[i + 1]; - for (idx j = xb; j < xe; ++j) { + idx xe = xadj[i+1]; + for (idx j = xb; j < xe; ++j){ idx dst = adj[j]; - if (i < dst) { - lower_triangle_srcs[ind] = i; + if (i < dst){ + lower_triangle_srcs[ind] = i; lower_triangle_dests[ind++] = dst; } } @@ -569,45 +596,50 @@ void convert_crs_to_lower_triangle_edge_list(idx nv, idx *xadj, idx *adj, } template -void convert_crs_to_edge_list(idx nv, idx *xadj, idx *srcs) { - for (idx i = 0; i < nv; ++i) { +void convert_crs_to_edge_list(idx nv, idx *xadj, idx *srcs){ + for (idx i = 0; i < nv; ++i){ idx xb = xadj[i]; - idx xe = xadj[i + 1]; - for (idx j = xb; j < xe; ++j) { + idx xe = xadj[i+1]; + for (idx j = xb; j < xe; ++j){ srcs[j] = i; } } } template -void convert_edge_list_to_csr(lno_t nv, size_type ne, lno_t *srcs, lno_t *dests, - wt *ew, size_type *xadj, lno_t *adj, wt *crs_ew) { - std::vector> edges(ne); - for (size_type i = 0; i < ne; ++i) { +void convert_edge_list_to_csr (lno_t nv, size_type ne, + lno_t *srcs, lno_t *dests, wt *ew, + size_type *xadj, lno_t *adj, wt *crs_ew){ + + std::vector > edges (ne); + for(size_type i = 0; i < ne; ++i){ edges[i].src = srcs[i]; edges[i].dst = dests[i]; - edges[i].ew = ew[i]; + edges[i].ew = ew[i]; } - std::sort(edges.begin(), edges.begin() + ne); + std::sort (edges.begin(), edges.begin() + ne); size_type eind = 0; - for (lno_t i = 0; i < nv; ++i) { + for (lno_t i = 0; i < nv; ++i){ (xadj)[i] = eind; - while (edges[eind].src == i) { - (adj)[eind] = edges[eind].dst; + while (edges[eind].src == i){ + (adj)[eind] = edges[eind].dst; (*crs_ew)[eind] = edges[eind].ew; ++eind; } } xadj[nv] = eind; + } template -void convert_undirected_edge_list_to_csr(lno_t nv, size_type ne, in_lno_t *srcs, - in_lno_t *dests, size_type *xadj, - lno_t *adj) { - std::vector> edges(ne * 2); - for (size_type i = 0; i < ne; ++i) { +void convert_undirected_edge_list_to_csr ( + lno_t nv, size_type ne, + in_lno_t *srcs, in_lno_t *dests, + size_type *xadj, lno_t *adj){ + + std::vector > edges (ne * 2); + for(size_type i = 0; i < ne; ++i){ edges[i * 2].src = srcs[i]; edges[i * 2].dst = dests[i]; @@ -615,21 +647,22 @@ void convert_undirected_edge_list_to_csr(lno_t nv, size_type ne, in_lno_t *srcs, edges[i * 2 + 1].dst = srcs[i]; } #ifdef KOKKOSKERNELS_HAVE_OUTER -#include -#include -#include -#include - __gnu_parallel::parallel_sort_mwms *>( - &(edges[0]), &(edges[0]) + ne * 2, - std::less>(), 64); +#include +#include +#include +#include + __gnu_parallel::parallel_sort_mwms *> + (&(edges[0]), &(edges[0])+ne*2, + std::less >(), 64); #else - std::sort(edges.begin(), edges.begin() + ne * 2); + std::sort (edges.begin(), edges.begin() + ne * 2); #endif + size_type eind = 0; - for (lno_t i = 0; i < nv; ++i) { + for (lno_t i = 0; i < nv; ++i){ (xadj)[i] = eind; - while (edges[eind].src == i) { + while (edges[eind].src == i){ (adj)[eind] = edges[eind].dst; //(*crs_ew)[eind] = edges[eind].ew; ++eind; @@ -674,188 +707,216 @@ void read_graph_src_dst_bin( } */ + template -void write_edgelist_bin(size_t ne, const idx *edge_begins, const idx *edge_ends, - const wt *ew, const char *filename) { - std::ofstream myFile(filename, std::ios::out | std::ios::binary); - myFile.write((char *)&ne, sizeof(idx)); - myFile.write((char *)edge_begins, sizeof(idx) * (ne)); - myFile.write((char *)edge_ends, sizeof(idx) * (ne)); - myFile.write((char *)ew, sizeof(wt) * (ne)); +void write_edgelist_bin( + size_t ne, + const idx *edge_begins, + const idx *edge_ends, + const wt *ew, + const char *filename){ + std::ofstream myFile (filename, std::ios::out | std::ios::binary); + myFile.write((char *) &ne, sizeof(idx)); + myFile.write((char *) edge_begins, sizeof(idx) * (ne)); + myFile.write((char *) edge_ends, sizeof(idx) * (ne)); + myFile.write((char *) ew, sizeof(wt) * (ne)); myFile.close(); } template -void read_edgelist_bin(idx *ne, idx **edge_begins, idx **edge_ends, wt **ew, - const char *filename) { - std::ifstream myFile(filename, std::ios::in | std::ios::binary); +void read_edgelist_bin( + idx *ne, + idx **edge_begins, + idx **edge_ends, + wt **ew, + const char *filename){ + + std::ifstream myFile (filename, std::ios::in | std::ios::binary); - myFile.read((char *)ne, sizeof(idx)); + + myFile.read((char *) ne, sizeof(idx)); md_malloc(edge_begins, *ne); md_malloc(edge_ends, *ne); - md_malloc(ew, *ne); - myFile.read((char *)*edge_begins, sizeof(idx) * (*ne)); - myFile.read((char *)*edge_ends, sizeof(idx) * (*ne)); - myFile.read((char *)*ew, sizeof(wt) * (*ne)); + md_malloc (ew, *ne); + myFile.read((char *) *edge_begins, sizeof(idx) * (*ne)); + myFile.read((char *) *edge_ends, sizeof(idx) * (*ne)); + myFile.read((char *) *ew, sizeof(wt) * (*ne)); myFile.close(); } + + + + template -void write_graph_bin(lno_t nv, size_type ne, const size_type *xadj, - const lno_t *adj, const scalar_t *ew, - const char *filename) { - std::ofstream myFile(filename, std::ios::out | std::ios::binary); - myFile.write((char *)&nv, sizeof(lno_t)); - myFile.write((char *)&ne, sizeof(size_type)); - myFile.write((char *)xadj, sizeof(size_type) * (nv + 1)); +void write_graph_bin(lno_t nv, size_type ne,const size_type *xadj,const lno_t *adj,const scalar_t *ew,const char *filename){ + std::ofstream myFile (filename, std::ios::out | std::ios::binary); + myFile.write((char *) &nv, sizeof(lno_t)); + myFile.write((char *) &ne, sizeof(size_type)); + myFile.write((char *) xadj, sizeof(size_type) * (nv + 1)); - myFile.write((char *)adj, sizeof(lno_t) * (ne)); + myFile.write((char *) adj, sizeof(lno_t) * (ne)); - myFile.write((char *)ew, sizeof(scalar_t) * (ne)); + myFile.write((char *) ew, sizeof(scalar_t) * (ne)); myFile.close(); } template -void write_graph_crs(lno_t nv, size_type ne, const size_type *xadj, - const lno_t *adj, const scalar_t *ew, - const char *filename) { - std::ofstream myFile(filename, std::ios::out); +void write_graph_crs(lno_t nv, size_type ne,const size_type *xadj,const lno_t *adj,const scalar_t *ew,const char *filename){ + std::ofstream myFile (filename, std::ios::out ); myFile << nv << " " << ne << std::endl; - for (lno_t i = 0; i <= nv; ++i) { - myFile << xadj[i] << " "; + for (lno_t i = 0; i <= nv; ++i){ + myFile << xadj[i] << " "; } - myFile << std::endl; + myFile << std::endl; - for (lno_t i = 0; i < nv; ++i) { + for (lno_t i = 0; i < nv; ++i){ size_type b = xadj[i]; size_type e = xadj[i + 1]; - for (size_type j = b; j < e; ++j) { - myFile << adj[j] << " "; + for (size_type j = b; j < e; ++j){ + myFile << adj[j] << " "; } - myFile << std::endl; + myFile << std::endl; } - for (size_type i = 0; i < ne; ++i) { - myFile << ew[i] << " "; + for (size_type i = 0; i < ne; ++i){ + myFile << ew[i] << " "; } - myFile << std::endl; + myFile << std::endl; myFile.close(); } template -void write_graph_ligra(lno_t nv, size_type ne, const size_type *xadj, - const lno_t *adj, const scalar_t * /*ew*/, - const char *filename) { - std::ofstream ff(filename); +void write_graph_ligra(lno_t nv, size_type ne,const size_type *xadj,const lno_t *adj,const scalar_t * /*ew*/,const char *filename){ + + std::ofstream ff (filename); ff << "AdjacencyGraph" << std::endl; ff << nv << std::endl << ne << std::endl; - for (lno_t i = 0; i < nv; ++i) { + for (lno_t i = 0; i < nv; ++i){ ff << xadj[i] << std::endl; } - for (size_type i = 0; i < ne; ++i) { + for (size_type i = 0; i < ne; ++i){ ff << adj[i] << std::endl; } ff.close(); } -// MM: types and utility functions for parsing the MatrixMarket format -namespace MM { -enum MtxObject { UNDEFINED_OBJECT, MATRIX, VECTOR }; -enum MtxFormat { UNDEFINED_FORMAT, COORDINATE, ARRAY }; -enum MtxField { - UNDEFINED_FIELD, - REAL, // includes both float and double - COMPLEX, // includes complex and complex - INTEGER, // includes all integer types - PATTERN // not a type, but means the value for every entry is 1 -}; -enum MtxSym { - UNDEFINED_SYMMETRY, - GENERAL, - SYMMETRIC, // A(i, j) = A(j, i) - SKEW_SYMMETRIC, // A(i, j) = -A(j, i) - HERMITIAN // A(i, j) = a + bi; A(j, i) = a - bi -}; - -// readScalar/writeScalar: read and write a scalar in the form that it appears -// in an .mtx file. The >> and << operators won't work, because complex appears -// as "real imag", not "(real, imag)" -template -scalar_t readScalar(std::istream &is) { - scalar_t val; - is >> val; - return val; -} +//MM: types and utility functions for parsing the MatrixMarket format +namespace MM +{ + enum MtxObject + { + UNDEFINED_OBJECT, + MATRIX, + VECTOR + }; + enum MtxFormat + { + UNDEFINED_FORMAT, + COORDINATE, + ARRAY + }; + enum MtxField + { + UNDEFINED_FIELD, + REAL, //includes both float and double + COMPLEX, //includes complex and complex + INTEGER, //includes all integer types + PATTERN //not a type, but means the value for every entry is 1 + }; + enum MtxSym + { + UNDEFINED_SYMMETRY, + GENERAL, + SYMMETRIC, //A(i, j) = A(j, i) + SKEW_SYMMETRIC, //A(i, j) = -A(j, i) + HERMITIAN //A(i, j) = a + bi; A(j, i) = a - bi + }; + + //readScalar/writeScalar: read and write a scalar in the form that it appears in an .mtx file. + //The >> and << operators won't work, because complex appears as "real imag", not "(real, imag)" + template + scalar_t readScalar(std::istream& is) + { + scalar_t val; + is >> val; + return val; + } -template <> -inline Kokkos::complex readScalar(std::istream &is) { - float r, i; - is >> r; - is >> i; - return Kokkos::complex(r, i); -} + template<> + inline Kokkos::complex readScalar(std::istream& is) + { + float r, i; + is >> r; + is >> i; + return Kokkos::complex(r, i); + } -template <> -inline Kokkos::complex readScalar(std::istream &is) { - double r, i; - is >> r; - is >> i; - return Kokkos::complex(r, i); -} + template<> + inline Kokkos::complex readScalar(std::istream& is) + { + double r, i; + is >> r; + is >> i; + return Kokkos::complex(r, i); + } -template -void writeScalar(std::ostream &os, scalar_t val) { - os << val; -} + template + void writeScalar(std::ostream& os, scalar_t val) + { + os << val; + } -template <> -inline void writeScalar(std::ostream &os, Kokkos::complex val) { - os << val.real() << ' ' << val.imag(); -} + template<> + inline void writeScalar(std::ostream& os, Kokkos::complex val) + { + os << val.real() << ' ' << val.imag(); + } -template <> -inline void writeScalar(std::ostream &os, Kokkos::complex val) { - os << val.real() << ' ' << val.imag(); -} + template<> + inline void writeScalar(std::ostream& os, Kokkos::complex val) + { + os << val.real() << ' ' << val.imag(); + } -// symmetryFlip: given a value for A(i, j), return the value that -// should be inserted at A(j, i) (if any) -template -scalar_t symmetryFlip(scalar_t val, MtxSym symFlag) { - if (symFlag == SKEW_SYMMETRIC) return -val; - return val; -} + //symmetryFlip: given a value for A(i, j), return the value that + //should be inserted at A(j, i) (if any) + template + scalar_t symmetryFlip(scalar_t val, MtxSym symFlag) + { + if(symFlag == SKEW_SYMMETRIC) + return -val; + return val; + } -template <> -inline Kokkos::complex symmetryFlip(Kokkos::complex val, - MtxSym symFlag) { - if (symFlag == HERMITIAN) - return Kokkos::conj(val); - else if (symFlag == SKEW_SYMMETRIC) - return -val; - return val; -} + template<> + inline Kokkos::complex symmetryFlip(Kokkos::complex val, MtxSym symFlag) + { + if(symFlag == HERMITIAN) + return Kokkos::conj(val); + else if(symFlag == SKEW_SYMMETRIC) + return -val; + return val; + } -template <> -inline Kokkos::complex symmetryFlip(Kokkos::complex val, - MtxSym symFlag) { - if (symFlag == HERMITIAN) - return Kokkos::conj(val); - else if (symFlag == SKEW_SYMMETRIC) - return -val; - return val; + template<> + inline Kokkos::complex symmetryFlip(Kokkos::complex val, MtxSym symFlag) + { + if(symFlag == HERMITIAN) + return Kokkos::conj(val); + else if(symFlag == SKEW_SYMMETRIC) + return -val; + return val; + } } -} // namespace MM template -void write_matrix_mtx(lno_t nrows, lno_t ncols, size_type nentries, - const size_type *xadj, const lno_t *adj, - const scalar_t *vals, const char *filename) { - std::ofstream myFile(filename); - myFile << "%%MatrixMarket matrix coordinate "; - if (std::is_same>::value || +void write_matrix_mtx(lno_t nrows, lno_t ncols, size_type nentries, const size_type *xadj, const lno_t *adj, const scalar_t *vals, const char *filename) { + std::ofstream myFile (filename); + myFile << "%%MatrixMarket matrix coordinate "; + if(std::is_same>::value || std::is_same>::value) myFile << "complex"; else @@ -867,7 +928,7 @@ void write_matrix_mtx(lno_t nrows, lno_t ncols, size_type nentries, size_type b = xadj[i]; size_type e = xadj[i + 1]; for (size_type j = b; j < e; ++j) { - myFile << i + 1 << " " << adj[j] + 1 << " "; + myFile << i + 1 << " " << adj[j] + 1 << " "; MM::writeScalar(myFile, vals[j]); myFile << '\n'; } @@ -876,24 +937,23 @@ void write_matrix_mtx(lno_t nrows, lno_t ncols, size_type nentries, } template -void write_graph_mtx(lno_t nv, size_type ne, const size_type *xadj, - const lno_t *adj, const scalar_t *ew, - const char *filename) { - std::ofstream myFile(filename); - myFile << "%%MatrixMarket matrix coordinate "; - if (std::is_same>::value || +void write_graph_mtx(lno_t nv, size_type ne,const size_type *xadj,const lno_t *adj,const scalar_t *ew,const char *filename){ + + std::ofstream myFile (filename); + myFile << "%%MatrixMarket matrix coordinate "; + if(std::is_same>::value || std::is_same>::value) myFile << "complex"; else myFile << "real"; myFile << " general\n"; - myFile << nv << " " << nv << " " << ne << '\n'; + myFile << nv << " " << nv << " " << ne << '\n'; myFile << std::setprecision(8) << std::scientific; - for (lno_t i = 0; i < nv; ++i) { + for (lno_t i = 0; i < nv; ++i){ size_type b = xadj[i]; size_type e = xadj[i + 1]; - for (size_type j = b; j < e; ++j) { - myFile << i + 1 << " " << (adj)[j] + 1 << " "; + for (size_type j = b; j < e; ++j){ + myFile << i + 1 << " " << (adj)[j] + 1 << " "; MM::writeScalar(myFile, ew[j]); myFile << '\n'; } @@ -902,344 +962,353 @@ void write_graph_mtx(lno_t nv, size_type ne, const size_type *xadj, myFile.close(); } + + template -void read_graph_bin(lno_t *nv, size_type *ne, size_type **xadj, lno_t **adj, - scalar_t **ew, const char *filename) { - std::ifstream myFile(filename, std::ios::in | std::ios::binary); +void read_graph_bin(lno_t *nv, size_type *ne,size_type **xadj, lno_t **adj, scalar_t **ew, const char *filename){ + + std::ifstream myFile (filename, std::ios::in | std::ios::binary); - myFile.read((char *)nv, sizeof(lno_t)); - myFile.read((char *)ne, sizeof(size_type)); - md_malloc(xadj, *nv + 1); + myFile.read((char *) nv, sizeof(lno_t)); + myFile.read((char *) ne, sizeof(size_type)); + md_malloc(xadj, *nv+1); md_malloc(adj, *ne); - md_malloc(ew, *ne); - myFile.read((char *)*xadj, sizeof(size_type) * (*nv + 1)); - myFile.read((char *)*adj, sizeof(lno_t) * (*ne)); - myFile.read((char *)*ew, sizeof(scalar_t) * (*ne)); + md_malloc (ew, *ne); + myFile.read((char *) *xadj, sizeof(size_type) * (*nv + 1)); + myFile.read((char *) *adj, sizeof(lno_t) * (*ne)); + myFile.read((char *) *ew, sizeof(scalar_t) * (*ne)); myFile.close(); } -// When Kokkos issue #2313 is resolved, can delete -// parseScalar and just use operator>> -template -scalar_t parseScalar(std::istream &is) { +//When Kokkos issue #2313 is resolved, can delete +//parseScalar and just use operator>> +template +scalar_t parseScalar(std::istream& is) +{ scalar_t val; is >> val; return val; } -template <> -inline Kokkos::complex parseScalar(std::istream &is) { +template<> +inline Kokkos::complex parseScalar(std::istream& is) +{ std::complex val; is >> val; return Kokkos::complex(val); } -template <> -inline Kokkos::complex parseScalar(std::istream &is) { +template<> +inline Kokkos::complex parseScalar(std::istream& is) +{ std::complex val; is >> val; return Kokkos::complex(val); } template -void read_graph_crs(lno_t *nv, size_type *ne, size_type **xadj, lno_t **adj, - scalar_t **ew, const char *filename) { - std::ifstream myFile(filename, std::ios::in); +void read_graph_crs(lno_t *nv, size_type *ne,size_type **xadj, lno_t **adj, scalar_t **ew, const char *filename){ + + std::ifstream myFile (filename, std::ios::in ); myFile >> *nv >> *ne; - md_malloc(xadj, *nv + 1); + md_malloc(xadj, *nv+1); md_malloc(adj, *ne); - md_malloc(ew, *ne); + md_malloc (ew, *ne); - for (lno_t i = 0; i <= *nv; ++i) { - myFile >> (*xadj)[i]; + for (lno_t i = 0; i <= *nv; ++i){ + myFile >> (*xadj)[i]; } - for (size_type i = 0; i < *ne; ++i) { - myFile >> (*adj)[i]; + for (size_type i = 0; i < *ne; ++i){ + myFile >> (*adj)[i]; } - for (size_type i = 0; i < *ne; ++i) { + for (size_type i = 0; i < *ne; ++i){ (*ew)[i] = parseScalar(myFile); } myFile.close(); } -inline bool endswith(std::string const &fullString, std::string const &ending) { - if (fullString.length() >= ending.length()) { - return (0 == fullString.compare(fullString.length() - ending.length(), - ending.length(), ending)); - } else { - return false; - } + + + +inline bool endswith (std::string const &fullString, std::string const &ending) { + if (fullString.length() >= ending.length()) { + return (0 == fullString.compare (fullString.length() - ending.length(), ending.length(), ending)); + } else { + return false; + } } + template -void write_kokkos_crst_matrix(crs_matrix_t a_crsmat, const char *filename) { +void write_kokkos_crst_matrix(crs_matrix_t a_crsmat,const char *filename){ typedef typename crs_matrix_t::StaticCrsGraphType graph_t; - typedef typename graph_t::row_map_type::non_const_type row_map_view_t; - typedef typename graph_t::entries_type::non_const_type cols_view_t; + typedef typename graph_t::row_map_type::non_const_type row_map_view_t; + typedef typename graph_t::entries_type::non_const_type cols_view_t; typedef typename crs_matrix_t::values_type::non_const_type values_view_t; typedef typename row_map_view_t::value_type offset_t; - typedef typename cols_view_t::value_type lno_t; - typedef typename values_view_t::value_type scalar_t; - typedef typename values_view_t::size_type size_type; + typedef typename cols_view_t::value_type lno_t; + typedef typename values_view_t::value_type scalar_t; + typedef typename values_view_t::size_type size_type; size_type nnz = a_crsmat.nnz(); - auto a_rowmap_view = Kokkos::create_mirror_view_and_copy( - Kokkos::HostSpace(), a_crsmat.graph.row_map); - auto a_entries_view = Kokkos::create_mirror_view_and_copy( - Kokkos::HostSpace(), a_crsmat.graph.entries); - auto a_values_view = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), a_crsmat.values); - offset_t *a_rowmap = const_cast(a_rowmap_view.data()); - lno_t *a_entries = a_entries_view.data(); - scalar_t *a_values = a_values_view.data(); + auto a_rowmap_view = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), a_crsmat.graph.row_map); + auto a_entries_view = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), a_crsmat.graph.entries); + auto a_values_view = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), a_crsmat.values); + offset_t* a_rowmap = const_cast(a_rowmap_view.data()); + lno_t* a_entries = a_entries_view.data(); + scalar_t* a_values = a_values_view.data(); std::string strfilename(filename); - if (endswith(strfilename, ".mtx") || endswith(strfilename, ".mm")) { + if (endswith(strfilename, ".mtx") || endswith(strfilename, ".mm")){ write_matrix_mtx( - a_crsmat.numRows(), a_crsmat.numCols(), a_crsmat.nnz(), a_rowmap, - a_entries, a_values, filename); + a_crsmat.numRows(), a_crsmat.numCols(), a_crsmat.nnz(), + a_rowmap, a_entries, a_values, filename); return; - } else if (a_crsmat.numRows() != a_crsmat.numCols()) { - throw std::runtime_error( - "For formats other than MatrixMarket (suffix .mm or .mtx),\n" + } + else if(a_crsmat.numRows() != a_crsmat.numCols()) + { + throw std::runtime_error("For formats other than MatrixMarket (suffix .mm or .mtx),\n" "write_kokkos_crst_matrix only supports square matrices"); } - if (endswith(strfilename, ".bin")) { - write_graph_bin( - a_crsmat.numRows(), nnz, a_rowmap, a_entries, a_values, filename); - } else if (endswith(strfilename, ".ligra")) { - write_graph_ligra( - a_crsmat.numRows(), nnz, a_rowmap, a_entries, a_values, filename); - } else if (endswith(strfilename, ".crs")) { - write_graph_crs( - a_crsmat.numRows(), nnz, a_rowmap, a_entries, a_values, filename); - } else { - std::string errMsg = - std::string("write_kokkos_crst_matrix: File extension on ") + filename + - " does not correspond to a known format"; + if (endswith(strfilename, ".bin")){ + write_graph_bin(a_crsmat.numRows(), + nnz, a_rowmap, a_entries, a_values, filename); + } + else if (endswith(strfilename, ".ligra")){ + write_graph_ligra(a_crsmat.numRows(), + nnz, a_rowmap, a_entries, a_values, filename); + } + else if (endswith(strfilename, ".crs")){ + write_graph_crs(a_crsmat.numRows(), + nnz, a_rowmap, a_entries, a_values, filename); + } + else { + std::string errMsg = std::string("write_kokkos_crst_matrix: File extension on ") + filename + " does not correspond to a known format"; throw std::runtime_error(errMsg); } } template -int read_mtx(const char *fileName, lno_t *nrows, lno_t *ncols, size_type *ne, - size_type **xadj, lno_t **adj, scalar_t **ew, - bool symmetrize = false, bool remove_diagonal = true, - bool transpose = false) { +int read_mtx ( + const char *fileName, + lno_t *nrows, lno_t* ncols, size_type *ne, + size_type **xadj, lno_t **adj, scalar_t **ew, + bool symmetrize = false, bool remove_diagonal = true, + bool transpose = false) +{ using namespace MM; - std::ifstream mmf(fileName, std::ifstream::in); + std::ifstream mmf (fileName, std::ifstream::in); if (!mmf.is_open()) { - throw std::runtime_error("File cannot be opened\n"); + throw std::runtime_error ("File cannot be opened\n"); } std::string fline = ""; getline(mmf, fline); - if (fline.size() < 2 || fline[0] != '%' || fline[1] != '%') { - throw std::runtime_error("Invalid MM file. Line-1\n"); + if (fline.size() < 2 || fline[0] != '%' || fline[1] != '%'){ + throw std::runtime_error ("Invalid MM file. Line-1\n"); } - // make sure every required field is in the file, by initializing them to - // UNDEFINED_* + //make sure every required field is in the file, by initializing them to UNDEFINED_* MtxObject mtx_object = UNDEFINED_OBJECT; MtxFormat mtx_format = UNDEFINED_FORMAT; - MtxField mtx_field = UNDEFINED_FIELD; - MtxSym mtx_sym = UNDEFINED_SYMMETRY; + MtxField mtx_field = UNDEFINED_FIELD; + MtxSym mtx_sym = UNDEFINED_SYMMETRY; - if (fline.find("matrix") != std::string::npos) { + if (fline.find("matrix") != std::string::npos){ mtx_object = MATRIX; - } else if (fline.find("vector") != std::string::npos) { + } else if (fline.find("vector") != std::string::npos){ mtx_object = VECTOR; - throw std::runtime_error( - "MatrixMarket \"vector\" is not supported by KokkosKernels read_mtx()"); + throw std::runtime_error("MatrixMarket \"vector\" is not supported by KokkosKernels read_mtx()"); } - if (fline.find("coordinate") != std::string::npos) { - // sparse + if (fline.find("coordinate") != std::string::npos){ + //sparse mtx_format = COORDINATE; - } else if (fline.find("array") != std::string::npos) { - // dense + } + else if (fline.find("array") != std::string::npos){ + //dense mtx_format = ARRAY; } - if (fline.find("real") != std::string::npos || - fline.find("double") != std::string::npos) { - if (std::is_same::value) + if(fline.find("real") != std::string::npos || + fline.find("double") != std::string::npos) + { + if (std::is_same::value) mtx_field = REAL; else { if (!std::is_floating_point::value) throw std::runtime_error( - "scalar_t in read_mtx() incompatible with float or double typed " - "MatrixMarket file."); + "scalar_t in read_mtx() incompatible with float or double typed MatrixMarket file."); else mtx_field = REAL; } - } else if (fline.find("complex") != std::string::npos) { - if (!(std::is_same>::value || + } + else if (fline.find("complex") != std::string::npos){ + if(!(std::is_same>::value || std::is_same>::value)) - throw std::runtime_error( - "scalar_t in read_mtx() incompatible with complex-typed MatrixMarket " - "file."); + throw std::runtime_error("scalar_t in read_mtx() incompatible with complex-typed MatrixMarket file."); else mtx_field = COMPLEX; - } else if (fline.find("integer") != std::string::npos) { - if (std::is_integral::value || - std::is_floating_point::value || - std::is_same::value) + } + else if (fline.find("integer") != std::string::npos){ + if(std::is_integral::value) mtx_field = INTEGER; else - throw std::runtime_error( - "scalar_t in read_mtx() incompatible with integer-typed MatrixMarket " - "file."); - } else if (fline.find("pattern") != std::string::npos) { + throw std::runtime_error("scalar_t in read_mtx() incompatible with integer-typed MatrixMarket file."); + } + else if (fline.find("pattern") != std::string::npos){ mtx_field = PATTERN; - // any reasonable choice for scalar_t can represent "1" or "1.0 + 0i", so - // nothing to check here + //any reasonable choice for scalar_t can represent "1" or "1.0 + 0i", so nothing to check here } - if (fline.find("general") != std::string::npos) { + if (fline.find("general") != std::string::npos){ mtx_sym = GENERAL; - } else if (fline.find("skew-symmetric") != std::string::npos) { + } + else if (fline.find("skew-symmetric") != std::string::npos){ mtx_sym = SKEW_SYMMETRIC; - } else if (fline.find("symmetric") != std::string::npos) { - // checking for "symmetric" after "skew-symmetric" because it's a substring + } + else if (fline.find("symmetric") != std::string::npos){ + //checking for "symmetric" after "skew-symmetric" because it's a substring mtx_sym = SYMMETRIC; - } else if (fline.find("hermitian") != std::string::npos || - fline.find("Hermitian") != std::string::npos) { + } + else if (fline.find("hermitian") != std::string::npos || + fline.find("Hermitian") != std::string::npos){ mtx_sym = HERMITIAN; } - // Validate the matrix attributes - if (mtx_format == ARRAY) { - if (mtx_sym == UNDEFINED_SYMMETRY) mtx_sym = GENERAL; - if (mtx_sym != GENERAL) - throw std::runtime_error( - "array format MatrixMarket file must have general symmetry (optional " - "to include \"general\")"); - } - if (mtx_object == UNDEFINED_OBJECT) - throw std::runtime_error( - "MatrixMarket file header is missing the object type."); - if (mtx_format == UNDEFINED_FORMAT) + //Validate the matrix attributes + if(mtx_format == ARRAY) + { + if(mtx_sym == UNDEFINED_SYMMETRY) + mtx_sym = GENERAL; + if(mtx_sym != GENERAL) + throw std::runtime_error("array format MatrixMarket file must have general symmetry (optional to include \"general\")"); + } + if(mtx_object == UNDEFINED_OBJECT) + throw std::runtime_error("MatrixMarket file header is missing the object type."); + if(mtx_format == UNDEFINED_FORMAT) throw std::runtime_error("MatrixMarket file header is missing the format."); - if (mtx_field == UNDEFINED_FIELD) - throw std::runtime_error( - "MatrixMarket file header is missing the field type."); - if (mtx_sym == UNDEFINED_SYMMETRY) - throw std::runtime_error( - "MatrixMarket file header is missing the symmetry type."); - - while (1) { + if(mtx_field == UNDEFINED_FIELD) + throw std::runtime_error("MatrixMarket file header is missing the field type."); + if(mtx_sym == UNDEFINED_SYMMETRY) + throw std::runtime_error("MatrixMarket file header is missing the symmetry type."); + + while(1){ getline(mmf, fline); - if (fline[0] != '%') break; + if(fline[0] != '%') break; } - std::stringstream ss(fline); + std::stringstream ss (fline); lno_t nr = 0, nc = 0; size_type nnz = 0; ss >> nr >> nc; - if (mtx_format == COORDINATE) + if(mtx_format == COORDINATE) ss >> nnz; else nnz = nr * nc; size_type numEdges = nnz; - symmetrize = symmetrize || mtx_sym != GENERAL; - if (symmetrize && nr != nc) { + symmetrize = symmetrize || mtx_sym != GENERAL; + if(symmetrize && nr != nc) + { throw std::runtime_error("A non-square matrix cannot be symmetrized."); } - if (mtx_format == ARRAY) { - // Array format only supports general symmetry and non-pattern - if (symmetrize) - throw std::runtime_error( - "array format MatrixMarket file cannot be symmetrized."); - if (mtx_field == PATTERN) - throw std::runtime_error( - "array format MatrixMarket file can't have \"pattern\" field type."); + if(mtx_format == ARRAY) + { + //Array format only supports general symmetry and non-pattern + if(symmetrize) + throw std::runtime_error("array format MatrixMarket file cannot be symmetrized."); + if(mtx_field == PATTERN) + throw std::runtime_error("array format MatrixMarket file can't have \"pattern\" field type."); } - if (symmetrize) { + if(symmetrize) + { numEdges = 2 * nnz; } - // numEdges is only an upper bound (diagonal entries may be removed) - std::vector> edges(numEdges); - size_type nE = 0; + //numEdges is only an upper bound (diagonal entries may be removed) + std::vector > edges (numEdges); + size_type nE = 0; lno_t numDiagonal = 0; - for (size_type i = 0; i < nnz; ++i) { + for (size_type i = 0; i < nnz; ++i){ getline(mmf, fline); - std::stringstream ss2(fline); + std::stringstream ss2 (fline); struct Edge tmp; - // read source, dest (edge) and weight (value) - lno_t s, d; + //read source, dest (edge) and weight (value) + lno_t s,d; scalar_t w; - if (mtx_format == ARRAY) { - // In array format, entries are listed in column major order, - // so the row and column can be determined just from the index i + if(mtx_format == ARRAY) + { + //In array format, entries are listed in column major order, + //so the row and column can be determined just from the index i //(but make them 1-based indices, to match the way coordinate works) - s = i % nr + 1; // row - d = i / nr + 1; // col - } else { - // In coordinate format, row and col of each entry is read from file + s = i % nr + 1; //row + d = i / nr + 1; //col + } + else + { + //In coordinate format, row and col of each entry is read from file ss2 >> s >> d; } - if (mtx_field == PATTERN) + if(mtx_field == PATTERN) w = 1; else w = readScalar(ss2); - if (!transpose) { + if (!transpose){ tmp.src = s - 1; tmp.dst = d - 1; - tmp.ew = w; - } else { + tmp.ew = w; + } + else { tmp.src = d - 1; tmp.dst = s - 1; - tmp.ew = w; + tmp.ew = w; } - if (tmp.src == tmp.dst) { + if (tmp.src == tmp.dst){ numDiagonal++; - if (!remove_diagonal) { + if (!remove_diagonal){ edges[nE++] = tmp; } continue; } edges[nE++] = tmp; - if (symmetrize) { + if (symmetrize){ struct Edge tmp2; tmp2.src = tmp.dst; tmp2.dst = tmp.src; - // the symmetrized value is w, -w or conj(w) if mtx_sym is - // SYMMETRIC, SKEW_SYMMETRIC or HERMITIAN, respectively. - tmp2.ew = symmetryFlip(tmp.ew, mtx_sym); + //the symmetrized value is w, -w or conj(w) if mtx_sym is + //SYMMETRIC, SKEW_SYMMETRIC or HERMITIAN, respectively. + tmp2.ew = symmetryFlip(tmp.ew, mtx_sym); edges[nE++] = tmp2; } } mmf.close(); - std::sort(edges.begin(), edges.begin() + nE); - if (transpose) { + std::sort (edges.begin(), edges.begin() + nE); + if (transpose){ lno_t tmp = nr; - nr = nc; - nc = tmp; + nr = nc; + nc = tmp; } - // idx *nv, idx *ne, idx **xadj, idx **adj, wt **wt + //idx *nv, idx *ne, idx **xadj, idx **adj, wt **wt *nrows = nr; *ncols = nc; - *ne = nE; + *ne = nE; //*xadj = new idx[nr + 1]; - md_malloc(xadj, nr + 1); + md_malloc(xadj, nr+1); //*adj = new idx[nE]; md_malloc(adj, nE); //*ew = new wt[nE]; md_malloc(ew, nE); - size_type eind = 0; + size_type eind = 0; size_type actual = 0; - for (lno_t i = 0; i < nr; ++i) { - (*xadj)[i] = actual; + for (lno_t i = 0; i < nr; ++i){ + (*xadj)[i] = actual; bool is_first = true; - while (eind < nE && edges[eind].src == i) { - if (is_first || !symmetrize || eind == 0 || - (eind > 0 && edges[eind - 1].dst != edges[eind].dst)) { + while (eind < nE && edges[eind].src == i){ + if (is_first || !symmetrize || eind == 0 || (eind > 0 && edges[eind - 1].dst != edges[eind].dst)){ (*adj)[actual] = edges[eind].dst; - (*ew)[actual] = edges[eind].ew; + (*ew)[actual] = edges[eind].ew; ++actual; } is_first = false; @@ -1247,71 +1316,78 @@ int read_mtx(const char *fileName, lno_t *nrows, lno_t *ncols, size_type *ne, } } (*xadj)[nr] = actual; - *ne = actual; + *ne = actual; return 0; } -// Version of read_mtx which does not capture the number of columns. -// This is the old interface; it's kept for backwards compatibility. +//Version of read_mtx which does not capture the number of columns. +//This is the old interface; it's kept for backwards compatibility. template -int read_mtx(const char *fileName, lno_t *nv, size_type *ne, size_type **xadj, - lno_t **adj, scalar_t **ew, bool symmetrize = false, - bool remove_diagonal = true, bool transpose = false) { - lno_t ncol; // will discard - return read_mtx(fileName, nv, &ncol, ne, xadj, - adj, ew, symmetrize, - remove_diagonal, transpose); +int read_mtx ( + const char *fileName, + lno_t *nv, size_type *ne, + size_type **xadj, lno_t **adj, scalar_t **ew, + bool symmetrize = false, bool remove_diagonal = true, + bool transpose = false) +{ + lno_t ncol; //will discard + return read_mtx(fileName, nv, &ncol, ne, xadj, adj, ew, symmetrize, remove_diagonal, transpose); } template -void read_matrix(lno_t *nv, size_type *ne, size_type **xadj, lno_t **adj, - scalar_t **ew, const char *filename) { +void read_matrix(lno_t *nv, size_type *ne,size_type **xadj, lno_t **adj, scalar_t **ew, const char *filename){ + std::string strfilename(filename); - if (endswith(strfilename, ".mtx") || endswith(strfilename, ".mm")) { - read_mtx(filename, nv, ne, xadj, adj, ew, false, false, false); + if (endswith(strfilename, ".mtx") || endswith(strfilename, ".mm")){ + read_mtx ( + filename, + nv, ne, + xadj, adj, ew,false,false,false); } - else if (endswith(strfilename, ".bin")) { - read_graph_bin(nv, ne, xadj, adj, ew, filename); + else if (endswith(strfilename, ".bin")){ + read_graph_bin(nv, ne,xadj, adj, ew, filename); } - else if (endswith(strfilename, ".crs")) { - read_graph_crs(nv, ne, xadj, adj, ew, filename); + else if (endswith(strfilename, ".crs")){ + + read_graph_crs(nv, ne,xadj, adj, ew, filename); } else { - throw std::runtime_error("Reader is not available\n"); + throw std::runtime_error ("Reader is not available\n"); } } template -crsMat_t read_kokkos_crst_matrix(const char *filename_) { +crsMat_t read_kokkos_crst_matrix(const char * filename_){ std::string strfilename(filename_); - bool isMatrixMarket = - endswith(strfilename, ".mtx") || endswith(strfilename, ".mm"); + bool isMatrixMarket = endswith(strfilename, ".mtx") || endswith(strfilename, ".mm"); typedef typename crsMat_t::StaticCrsGraphType graph_t; typedef typename graph_t::row_map_type::non_const_type row_map_view_t; - typedef typename graph_t::entries_type::non_const_type cols_view_t; + typedef typename graph_t::entries_type::non_const_type cols_view_t; typedef typename crsMat_t::values_type::non_const_type values_view_t; typedef typename row_map_view_t::value_type size_type; - typedef typename cols_view_t::value_type lno_t; + typedef typename cols_view_t::value_type lno_t; typedef typename values_view_t::value_type scalar_t; lno_t nr, nc, *adj; size_type *xadj, nnzA; scalar_t *values; - if (isMatrixMarket) { - // MatrixMarket file contains the exact number of columns - read_mtx(filename_, &nr, &nc, &nnzA, &xadj, - &adj, &values, false, false, false); - } else { - //.crs and .bin files don't contain #cols, so will compute it later based on - // the entries - read_matrix(&nr, &nnzA, &xadj, &adj, &values, - filename_); + if(isMatrixMarket) + { + //MatrixMarket file contains the exact number of columns + read_mtx ( + filename_, &nr, &nc, &nnzA, &xadj, &adj, &values, false, false, false); + } + else + { + //.crs and .bin files don't contain #cols, so will compute it later based on the entries + read_matrix( + &nr, &nnzA, &xadj, &adj, &values, filename_); } row_map_view_t rowmap_view("rowmap_view", nr + 1); @@ -1319,136 +1395,134 @@ crsMat_t read_kokkos_crst_matrix(const char *filename_) { values_view_t values_view("values_view", nnzA); { - Kokkos::View> - hr(xadj, nr + 1); - Kokkos::View> - hc(adj, nnzA); - Kokkos::View> - hv(values, nnzA); - Kokkos::deep_copy(rowmap_view, hr); - Kokkos::deep_copy(columns_view, hc); - Kokkos::deep_copy(values_view, hv); - } - - if (!isMatrixMarket) { - KokkosKernels::Impl::kk_view_reduce_max( - nnzA, columns_view, nc); - nc++; + Kokkos::View> hr(xadj, nr + 1); + Kokkos::View> hc(adj, nnzA); + Kokkos::View> hv(values, nnzA); + Kokkos::deep_copy (rowmap_view , hr); + Kokkos::deep_copy (columns_view , hc); + Kokkos::deep_copy (values_view , hv); } - graph_t static_graph(columns_view, rowmap_view); + if(!isMatrixMarket) + { + KokkosKernels::Impl::kk_view_reduce_max + (nnzA, columns_view, nc); + nc++; + } + + graph_t static_graph (columns_view, rowmap_view); crsMat_t crsmat("CrsMatrix", nc, values_view, static_graph); - delete[] xadj; - delete[] adj; - delete[] values; + delete [] xadj; delete [] adj; delete [] values; return crsmat; } + template -crsGraph_t read_kokkos_crst_graph(const char *filename_) { +crsGraph_t read_kokkos_crst_graph(const char * filename_){ + typedef typename crsGraph_t::row_map_type::non_const_type row_map_view_t; - typedef typename crsGraph_t::entries_type::non_const_type cols_view_t; + typedef typename crsGraph_t::entries_type::non_const_type cols_view_t; typedef typename row_map_view_t::value_type size_type; - typedef typename cols_view_t::value_type lno_t; - typedef double scalar_t; + typedef typename cols_view_t::value_type lno_t; + typedef double scalar_t ; lno_t nv, *adj; size_type *xadj, nnzA; scalar_t *values; - read_matrix(&nv, &nnzA, &xadj, &adj, &values, - filename_); + read_matrix( + &nv, &nnzA, &xadj, &adj, &values, filename_); - row_map_view_t rowmap_view("rowmap_view", nv + 1); + row_map_view_t rowmap_view("rowmap_view", nv+1); cols_view_t columns_view("colsmap_view", nnzA); + + { - typename row_map_view_t::HostMirror hr = - Kokkos::create_mirror_view(rowmap_view); - typename cols_view_t::HostMirror hc = - Kokkos::create_mirror_view(columns_view); + typename row_map_view_t::HostMirror hr = Kokkos::create_mirror_view (rowmap_view); + typename cols_view_t::HostMirror hc = Kokkos::create_mirror_view (columns_view); - for (lno_t i = 0; i <= nv; ++i) { + for (lno_t i = 0; i <= nv; ++i){ hr(i) = xadj[i]; } - for (size_type i = 0; i < nnzA; ++i) { + for (size_type i = 0; i < nnzA; ++i){ hc(i) = adj[i]; } - Kokkos::deep_copy(rowmap_view, hr); - Kokkos::deep_copy(columns_view, hc); + Kokkos::deep_copy (rowmap_view , hr); + Kokkos::deep_copy (columns_view , hc); } lno_t ncols = 0; - KokkosKernels::Impl::kk_view_reduce_max( - nnzA, columns_view, ncols); + KokkosKernels::Impl::kk_view_reduce_max + (nnzA, columns_view, ncols); ncols += 1; - crsGraph_t static_graph(columns_view, rowmap_view, ncols); - delete[] xadj; - delete[] adj; - delete[] values; + crsGraph_t static_graph (columns_view, rowmap_view, ncols); + delete [] xadj; delete [] adj; delete [] values; return static_graph; } + + template inline void kk_sequential_create_incidence_matrix( - nnz_lno_t num_rows, const size_type *xadj, const nnz_lno_t *adj, - size_type *i_adj // output. preallocated -) { + nnz_lno_t num_rows, + const size_type *xadj, + const nnz_lno_t *adj, + size_type *i_adj //output. preallocated + ){ + std::vector c_xadj(num_rows); - for (nnz_lno_t i = 0; i < num_rows; i++) { + for (nnz_lno_t i = 0; i < num_rows; i++){ c_xadj[i] = xadj[i]; } - int eCnt = 0; - for (nnz_lno_t i = 0; i < num_rows; i++) { - size_type begin = xadj[i]; - size_type end = xadj[i + 1]; + int eCnt=0; + for (nnz_lno_t i = 0; i < num_rows; i++){ + size_type begin = xadj[i]; + size_type end = xadj[i + 1]; nnz_lno_t adjsize = end - begin; - for (nnz_lno_t j = 0; j < adjsize; j++) { + for (nnz_lno_t j = 0; j < adjsize; j++){ size_type aind = j + begin; - nnz_lno_t col = adj[aind]; - if (i < col) { - i_adj[c_xadj[i]++] = eCnt; + nnz_lno_t col = adj[aind]; + if (i < col){ + i_adj[c_xadj[i]++] = eCnt; i_adj[c_xadj[col]++] = eCnt++; } } } - for (nnz_lno_t i = 0; i < num_rows; i++) { - if (c_xadj[i] != xadj[i + 1]) { - std::cout << "i:" << i << " c_xadj[i]:" << c_xadj[i] - << " xadj[i+1]:" << xadj[i + 1] << std::endl; + for (nnz_lno_t i = 0; i < num_rows; i++){ + if (c_xadj[i] != xadj[i+1]){ + std::cout << "i:" << i << " c_xadj[i]:" << c_xadj[i] << " xadj[i+1]:" << xadj[i+1] << std::endl; } } } template inline void kk_sequential_create_incidence_matrix_transpose( - const nnz_lno_t num_rows, const size_type num_edges, const size_type *xadj, + const nnz_lno_t num_rows, + const size_type num_edges, + const size_type *xadj, const nnz_lno_t *adj, - size_type *i_xadj, // output. preallocated - nnz_lno_t *i_adj // output. preallocated -) { - for (nnz_lno_t i = 0; i < num_edges / 2 + 1; i++) { + size_type *i_xadj, //output. preallocated + nnz_lno_t *i_adj //output. preallocated + ){ + + for (nnz_lno_t i = 0; i < num_edges/2 + 1; i++){ i_xadj[i] = i * 2; } - int eCnt = 0; - for (nnz_lno_t i = 0; i < num_rows; i++) { - size_type begin = xadj[i]; - size_type end = xadj[i + 1]; + int eCnt=0; + for (nnz_lno_t i = 0; i < num_rows; i++){ + size_type begin = xadj[i]; + size_type end = xadj[i + 1]; nnz_lno_t adjsize = end - begin; - for (nnz_lno_t j = 0; j < adjsize; j++) { + for (nnz_lno_t j = 0; j < adjsize; j++){ size_type aind = j + begin; - nnz_lno_t col = adj[aind]; - if (i < col) { + nnz_lno_t col = adj[aind]; + if (i < col){ i_adj[eCnt++] = i; i_adj[eCnt++] = col; } @@ -1456,7 +1530,8 @@ inline void kk_sequential_create_incidence_matrix_transpose( } } -} // namespace Impl -} // namespace KokkosKernels + +} +} #endif