From ddd4e1f42a62fc075f20c1c337a28fe367e2bc34 Mon Sep 17 00:00:00 2001 From: vasiloglou Date: Wed, 12 Dec 2007 20:40:53 +0000 Subject: [PATCH] sparse matrix works now, unit tests compile --- .../u/nvasil/sparse_matrix/sparse_matrix.h | 24 ++++++--- .../nvasil/sparse_matrix/sparse_matrix_impl.h | 49 ++++++++++++++----- .../sparse_matrix/sparse_matrix_test.cc | 39 +++++++++++---- 3 files changed, 83 insertions(+), 29 deletions(-) diff --git a/fastlib/u/nvasil/sparse_matrix/sparse_matrix.h b/fastlib/u/nvasil/sparse_matrix/sparse_matrix.h index 67e1f58ac0..b81964d72c 100644 --- a/fastlib/u/nvasil/sparse_matrix/sparse_matrix.h +++ b/fastlib/u/nvasil/sparse_matrix/sparse_matrix.h @@ -26,9 +26,15 @@ #include #include #include +#include #include #include "fastlib/fastlib.h" #include "la/matrix.h" +// you need this because trillinos redifines it. It's ok +// if you don't have it, but you will get an annoying warning +#ifdef F77_FUNC +#undef F77_FUNC +#endif #include "Epetra_CrsMatrix.h" #include "Epetra_SerialComm.h" #include "Epetra_Map.h" @@ -36,12 +42,10 @@ class SparseMatrix { public: SparseMatrix() ; - SparseMatrix(index_t num_of_rows, index_t num_of_columns) { - num_of_rows_ = num_of_rows; - num_of_columns_ = num_of_columns; - } // This constructor creates square matrices - SparseMatrix(index_t num_of_rows); + SparseMatrix(const index_t num_of_rows, + const index_t num_of_cols, + const index_t nnz_per_row); // Copy constructor SparseMatrix(const SparseMatrix &other); ~SparseMatrix() { @@ -49,8 +53,9 @@ class SparseMatrix { } void Destruct(); // Initializer: nnz_per_row is the estimated non zero elements per row - void Init(index_t num_of_rows, index_t num_of_columns, - index_t nnz_per_row); + void Init(const index_t num_of_rows, + const index_t num_of_columns, + const index_t nnz_per_row); // Initializer nnz_per_row has the estimated non zero elements per row // notice that it is not the same for every row; void Init(index_t num_of_rows, index_t num_of_columns, index_t *nnz_per_row); @@ -92,6 +97,11 @@ void Init(const std::vector &row_indices, void CopyValues(const SparseMatrix& other); double get(index_t r, index_t c); void set(index_t r, index_t c, double v); + std::string Print() { + std::ostringstream s1; + matrix_->Print(s1); + return s1.str(); + } index_t get_num_of_rows() { return num_of_rows_; } diff --git a/fastlib/u/nvasil/sparse_matrix/sparse_matrix_impl.h b/fastlib/u/nvasil/sparse_matrix/sparse_matrix_impl.h index 7cc14fb7b5..d561ffdbfa 100644 --- a/fastlib/u/nvasil/sparse_matrix/sparse_matrix_impl.h +++ b/fastlib/u/nvasil/sparse_matrix/sparse_matrix_impl.h @@ -18,20 +18,32 @@ */ SparseMatrix::SparseMatrix() { - + map_ = NULL; + matrix_ = NULL; } -void SparseMatrix::Init(index_t num_of_rows, - index_t num_of_columns, - index_t nnz_per_row) { +SparseMatrix::SparseMatrix(const index_t num_of_rows, + const index_t num_of_columns, + const index_t nnz_per_row) { + Init(num_of_rows, num_of_columns, nnz_per_row); +} +void SparseMatrix::Init(const index_t num_of_rows, + const index_t num_of_columns, + const index_t nnz_per_row) { + if likely(num_of_rows < num_of_columns) { + FATAL("Num of rows %i should be greater than the num of columns %i\n", + num_of_rows, num_of_columns); + } num_of_rows_=num_of_rows; num_of_columns_=num_of_columns; + dimension_ = num_of_rows_; if (num_of_rows_ < num_of_columns_) { FATAL("Num of rows %i is less than the number or columns %i", num_of_rows_, num_of_columns_); } map_ = new Epetra_Map(num_of_rows_, 0, comm_); matrix_ = new Epetra_CrsMatrix((Epetra_DataAccess)0, *map_, nnz_per_row); + StartLoadingRows(); } void SparseMatrix::Init(index_t num_of_rows, @@ -39,12 +51,14 @@ void SparseMatrix::Init(index_t num_of_rows, index_t *nnz_per_row) { num_of_rows_=num_of_rows; num_of_columns_=num_of_columns; + dimension_ = num_of_rows; if (num_of_rows_ < num_of_columns_) { FATAL("Num of rows %i is less than the number or columns %i", num_of_rows_, num_of_columns_); } map_ = new Epetra_Map(num_of_rows_, 0, comm_); matrix_ = new Epetra_CrsMatrix((Epetra_DataAccess)0 , *map_, nnz_per_row); + StartLoadingRows(); } @@ -69,7 +83,9 @@ void SparseMatrix::Init(const std::vector &rows, num_of_columns_ = columns[i]; } } - if ((index_t)frequencies.size()-1!=num_of_rows_) { + num_of_columns_++; + num_of_rows_++; + if ((index_t)frequencies.size()!=num_of_rows_) { NONFATAL("Some of the rows are zeros only!"); } index_t *nnz= new index_t[dimension_]; @@ -77,9 +93,9 @@ void SparseMatrix::Init(const std::vector &rows, nnz[i]=frequencies[i]; } Init(num_of_rows_, num_of_columns_, nnz); - delete []nnz; - Load(rows, columns, values); + //delete []nnz; } + Load(rows, columns, values); } void SparseMatrix::Init(const std::vector &rows, @@ -112,6 +128,9 @@ void SparseMatrix::Init(std::string filename) { vals.push_back(v); } fclose(fp); + /*for(index_t i=0; i< (index_t)rows.size(); i++) { + printf("%i %i %lg\n", rows[i], cols[i], vals[i]); + }*/ Init(rows, cols, vals, -1, -1); } @@ -176,6 +195,8 @@ void SparseMatrix::EndLoading() { void SparseMatrix::Load(const std::vector &rows, const std::vector &columns, const Vector &values) { + DEBUG_ASSERT(rows.size() ==columns.size()); + DEBUG_ASSERT((index_t)columns.size() == values.length()); my_global_elements_ = map_->MyGlobalElements(); index_t i=0; index_t cur_row = rows[i]; @@ -183,18 +204,22 @@ void SparseMatrix::Load(const std::vector &rows, std::vector indices; std::vector row_values; while (true) { - while (rows[cur_row]==rows[prev_row]) { + indices.clear(); + row_values.clear(); + while (likely((rows[cur_row]==rows[prev_row]) && + (i < (index_t)rows.size()))) { indices.push_back(columns[i]); row_values.push_back(values[i]); i++; prev_row=i-1; cur_row=i; } - matrix_->InsertGlobalValues(my_global_elements_[cur_row], - values.length(), + matrix_->InsertGlobalValues(my_global_elements_[rows[prev_row]], + row_values.size(), &row_values[0], &indices[0]); - if (i==dimension_) { + prev_row=cur_row; + if (i >= (index_t)rows.size()) { break; } } @@ -212,7 +237,7 @@ double SparseMatrix::get(index_t r, index_t c) { if (pos==indices+num_of_entries) { return 0; } - return values[*pos]; + return values[(ptrdiff_t)(pos-indices)]; } void SparseMatrix::set(index_t r, index_t c, double v) { diff --git a/fastlib/u/nvasil/sparse_matrix/sparse_matrix_test.cc b/fastlib/u/nvasil/sparse_matrix/sparse_matrix_test.cc index f5e4d5fa4a..4c8e81e745 100644 --- a/fastlib/u/nvasil/sparse_matrix/sparse_matrix_test.cc +++ b/fastlib/u/nvasil/sparse_matrix/sparse_matrix_test.cc @@ -44,7 +44,9 @@ class SparseMatrixTest { delete smat_; } void TestInit1() { - smat_ = new SparseMatrix(num_of_rows_, num_of_cols_); + smat_ = new SparseMatrix(num_of_rows_, + num_of_cols_, + num_of_nnz_); smat_->StartLoadingRows(); std::vector ind; std::vector val; @@ -53,12 +55,13 @@ class SparseMatrixTest { val.clear(); for(index_t j=0; jLoadRow(i, ind, val); } + for(index_t i=0; iget(i,j), @@ -66,30 +69,40 @@ class SparseMatrixTest { std::numeric_limits::epsilon()); } } + NONFATAL("TestInit1 sucess!!\n"); } void TestInit2() { std::vector rows; std::vector cols; std::vector vals; std::vector nnz(num_of_rows_); - for(index_t i=0; iInit(rows, cols, vals, *(std::max_element(nnz.begin(), nnz.end())), num_of_rows_); + // printf("%s\n", smat_->Print().c_str()); for(index_t i=0; iget(i,j), mat_[i][j], std::numeric_limits::epsilon()); } } + NONFATAL("TestInit2 success!!\n"); } void TestInit3() { FILE *fp = fopen("temp.txt", "w"); @@ -99,11 +112,13 @@ class SparseMatrixTest { for(index_t i=0; iInit("temp.txt"); + fclose(fp); + smat_ = new SparseMatrix(); + smat_->Init("temp.txt"); unlink("temp.txt"); for(index_t i=0; i::epsilon()); } } + NONFATAL("TestInit3 success!!"); } void TestCopyConstructor() { - } + smat_ = new SparseMatrix(); + NONFATAL("TestCopyConstructor success!!\n"); + } void TestMakeSymmetric() { TestInit1(); smat_->set(2, 3, 1.44); @@ -129,6 +147,7 @@ class SparseMatrixTest { std::numeric_limits::epsilon()); } } + NONFATAL("Test MakeSymmetric success!!\n"); } void TestAll(){ Init(); @@ -150,12 +169,12 @@ class SparseMatrixTest { private: SparseMatrix *smat_; - static const index_t num_of_cols_ = 40; - static const index_t num_of_rows_ = 40; + static const index_t num_of_cols_ = 10; + static const index_t num_of_rows_ = 10; + static const index_t num_of_nnz_ = 2; double mat_[num_of_rows_][num_of_cols_]; std::vector indices_; std::vector rows_; - Vector values_; }; int main() {