sparse matrix works now, unit tests compile

This commit is contained in:
vasiloglou
2007-12-12 20:40:53 +00:00
parent 042146937f
commit ddd4e1f42a
3 changed files with 83 additions and 29 deletions
+17 -7
View File
@@ -26,9 +26,15 @@
#include <string>
#include <map>
#include <vector>
#include <sstream>
#include <algorithm>
#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<index_t> &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_;
}
@@ -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<index_t> &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<index_t> &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<index_t> &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<index_t> &rows,
const std::vector<index_t> &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<index_t> &rows,
std::vector<index_t> indices;
std::vector<double> 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) {
@@ -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<index_t> ind;
std::vector<double> val;
@@ -53,12 +55,13 @@ class SparseMatrixTest {
val.clear();
for(index_t j=0; j<num_of_cols_; j++) {
if (mat_[i][j] != 0) {
ind.push_back(i);
ind.push_back(j);
val.push_back(mat_[i][j]);
}
}
smat_->LoadRow(i, ind, val);
}
for(index_t i=0; i<num_of_rows_; i++) {
for(index_t j=0; j<num_of_cols_; j++) {
TEST_DOUBLE_APPROX(smat_->get(i,j),
@@ -66,30 +69,40 @@ class SparseMatrixTest {
std::numeric_limits<double>::epsilon());
}
}
NONFATAL("TestInit1 sucess!!\n");
}
void TestInit2() {
std::vector<index_t> rows;
std::vector<index_t> cols;
std::vector<double> vals;
std::vector<index_t> nnz(num_of_rows_);
for(index_t i=0; i<num_of_cols_; i++) {
for(index_t i=0; i<num_of_rows_; i++) {
for(index_t j=0; j<num_of_cols_; j++) {
if (mat_[i][j] != 0) {
rows.push_back(i);
rows.push_back(j);
cols.push_back(j);
vals.push_back(mat_[i][j]);
nnz[i]++;
}
}
}
/*for(index_t i=0; i<(index_t)rows.size(); i++) {
printf("%i %i %lg\n",
rows[i],
cols[i],
vals[i]);
}*/
smat_ = new SparseMatrix();
smat_->Init(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; i<num_of_rows_; i++) {
for(index_t j=0; j<num_of_cols_; j++) {
TEST_DOUBLE_APPROX(smat_->get(i,j), mat_[i][j],
std::numeric_limits<double>::epsilon());
}
}
NONFATAL("TestInit2 success!!\n");
}
void TestInit3() {
FILE *fp = fopen("temp.txt", "w");
@@ -99,11 +112,13 @@ class SparseMatrixTest {
for(index_t i=0; i<num_of_cols_; i++) {
for(index_t j=0; j<num_of_cols_; j++) {
if (mat_[i][j] != 0) {
fprintf(fp, "%i %i %g", i, j, mat_[i][j]);
fprintf(fp, "%i %i %g\n", i, j, mat_[i][j]);
}
}
}
smat_->Init("temp.txt");
fclose(fp);
smat_ = new SparseMatrix();
smat_->Init("temp.txt");
unlink("temp.txt");
for(index_t i=0; i<num_of_rows_; i++) {
for(index_t j=0; j<num_of_cols_; j++) {
@@ -112,9 +127,12 @@ class SparseMatrixTest {
std::numeric_limits<double>::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<double>::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<index_t> indices_;
std::vector<index_t> rows_;
Vector values_;
};
int main() {