alternate coo2csr implementation

resolve/matrix/CMakeLists.txt

@@ -15,6 +15,7 @@ set(Matrix_SRC
  Coo.cpp
  MatrixHandler.cpp
  MatrixHandlerCpu.cpp
+ Utilities.cpp
 )
 
 # C++ code that depends on CUDA SDK libraries
@@ -35,6 +36,7 @@ set(Matrix_HEADER_INSTALL
  Csr.hpp
  Csc.hpp
  MatrixHandler.hpp
+ Utilities.hpp
 )
 
 # Build shared library ReSolve::matrix

resolve/matrix/Csr.cpp

@@ -7,6 +7,7 @@
 #include "Coo.hpp"
 #include <resolve/utilities/misc/IndexValuePair.hpp>
 #include <resolve/utilities/logger/Logger.hpp>
+#include <resolve/matrix/Utilities.hpp>
 
 namespace ReSolve 
 {
@@ -35,7 +36,7 @@ namespace ReSolve
  A_coo->symmetric(),
  A_coo->expanded())
  {
- coo2csr(A_coo, memspace);
+ matrix::coo2csr(A_coo, this, memspace);
  }
 
  /**
@@ -374,149 +375,10 @@ namespace ReSolve
  assert(nnz_ == A_coo->getNnz());
  assert(is_symmetric_ == A_coo->symmetric()); // <- Do we need to check for this?
 
- return coo2csr(A_coo, memspaceOut);
+ return matrix::coo2csr(A_coo, this, memspaceOut);
  }
 
 
- int matrix::Csr::coo2csr(matrix::Coo* A_coo, memory::MemorySpace memspace)
- {
- //count nnzs first
- index_type nnz_unpacked = 0;
- index_type nnz = A_coo->getNnz();
- index_type n = A_coo->getNumRows();
- bool symmetric = A_coo->symmetric();
- bool expanded = A_coo->expanded();
-
- index_type* nnz_counts = new index_type[n];
- std::fill_n(nnz_counts, n, 0);
- index_type* coo_rows = A_coo->getRowData(memory::HOST);
- index_type* coo_cols = A_coo->getColData(memory::HOST);
- real_type* coo_vals = A_coo->getValues( memory::HOST);
-
- index_type* diag_control = new index_type[n]; //for DEDUPLICATION of the diagonal
- std::fill_n(diag_control, n, 0);
- index_type nnz_unpacked_no_duplicates = 0;
- index_type nnz_no_duplicates = nnz;
-
-
- //maybe check if they exist?
- for (index_type i = 0; i < nnz; ++i)
- {
- nnz_counts[coo_rows[i]]++;
- nnz_unpacked++;
- nnz_unpacked_no_duplicates++;
- if ((coo_rows[i] != coo_cols[i])&& (symmetric) && (!expanded))
- {
- nnz_counts[coo_cols[i]]++;
- nnz_unpacked++;
- nnz_unpacked_no_duplicates++;
- }
- if (coo_rows[i] == coo_cols[i]){
- if (diag_control[coo_rows[i]] > 0) {
- //duplicate
- nnz_unpacked_no_duplicates--;
- nnz_no_duplicates--;
- }
- diag_control[coo_rows[i]]++;
- }
- }
- this->setExpanded(true);
- this->setNnzExpanded(nnz_unpacked_no_duplicates);
- index_type* csr_ia = new index_type[n+1];
- std::fill_n(csr_ia, n + 1, 0);
- index_type* csr_ja = new index_type[nnz_unpacked];
- real_type* csr_a = new real_type[nnz_unpacked];
- index_type* nnz_shifts = new index_type[n];
- std::fill_n(nnz_shifts, n , 0);
-
- IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; 
-
- csr_ia[0] = 0;
-
- for (index_type i = 1; i < n + 1; ++i){
- csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1);
- }
-
- int r, start;
-
-
- for (index_type i = 0; i < nnz; ++i){
- //which row
- r = coo_rows[i];
- start = csr_ia[r];
-
- if ((start + nnz_shifts[r]) > nnz_unpacked) {
- out::warning() << "index out of bounds (case 1) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl;
- }
- if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates
- bool already_there = false; 
- for (index_type j = start; j < start + nnz_shifts[r]; ++j)
- {
- index_type c = tmp[j].getIdx();
- if (c == r) {
- real_type val = tmp[j].getValue();
- val += coo_vals[i];
- tmp[j].setValue(val);
- already_there = true;
- out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl;
- } 
- } 
- if (!already_there){ // first time this duplicates appears
-
- tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]);
- tmp[start + nnz_shifts[r]].setValue(coo_vals[i]);
-
- nnz_shifts[r]++;
- }
- } else {//not diagonal
- tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]);
- tmp[start + nnz_shifts[r]].setValue(coo_vals[i]);
- nnz_shifts[r]++;
-
- if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1))
- {
- r = coo_cols[i];
- start = csr_ia[r];
-
- if ((start + nnz_shifts[r]) > nnz_unpacked)
- out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl;
- tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]);
- tmp[start + nnz_shifts[r]].setValue(coo_vals[i]);
- nnz_shifts[r]++;
- }
- }
- }
- //now sort whatever is inside rows
-
- for (int i = 0; i < n; ++i)
- {
- //now sorting (and adding 1)
- int colStart = csr_ia[i];
- int colEnd = csr_ia[i + 1];
- int length = colEnd - colStart;
- std::sort(&tmp[colStart],&tmp[colStart] + length);
- }
-
- for (index_type i = 0; i < nnz_unpacked; ++i)
- {
- csr_ja[i] = tmp[i].getIdx();
- csr_a[i] = tmp[i].getValue();
- }
-
- this->setNnz(nnz_no_duplicates);
- this->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace);
-
- delete [] nnz_counts;
- delete [] tmp;
- delete [] nnz_shifts;
- delete [] csr_ia;
- delete [] csr_ja;
- delete [] csr_a;
- delete [] diag_control; 
-
- return 0;
- }
-
  /**
  * @brief Prints matrix data.
  * 

resolve/matrix/Csr.hpp

@@ -49,8 +49,6 @@ namespace ReSolve { namespace matrix {
 
  int updateFromCoo(matrix::Coo* mat, memory::MemorySpace memspaceOut);
 
- private:
- int coo2csr(matrix::Coo* mat, memory::MemorySpace memspace);
  };
 
 }} // namespace ReSolve::matrix

resolve/matrix/MatrixHandler.cpp

@@ -6,7 +6,7 @@
 #include <resolve/matrix/Csc.hpp>
 #include <resolve/matrix/Csr.hpp>
 #include <resolve/workspace/LinAlgWorkspace.hpp>
-#include <resolve/utilities/misc/IndexValuePair.hpp>
+#include <resolve/matrix/Utilities.hpp>
 #include "MatrixHandler.hpp"
 #include "MatrixHandlerCpu.hpp"
 
@@ -116,150 +116,7 @@ namespace ReSolve {
  */
  int MatrixHandler::coo2csr(matrix::Coo* A_coo, matrix::Csr* A_csr, memory::MemorySpace memspace)
  {
- //count nnzs first
- index_type nnz_unpacked = 0;
- index_type nnz = A_coo->getNnz();
- index_type n = A_coo->getNumRows();
- bool symmetric = A_coo->symmetric();
- bool expanded = A_coo->expanded();
-
- index_type* nnz_counts = new index_type[n];
- std::fill_n(nnz_counts, n, 0);
- index_type* coo_rows = A_coo->getRowData(memory::HOST);
- index_type* coo_cols = A_coo->getColData(memory::HOST);
- real_type* coo_vals = A_coo->getValues( memory::HOST);
-
- index_type* diag_control = new index_type[n]; //for DEDUPLICATION of the diagonal
- std::fill_n(diag_control, n, 0);
- index_type nnz_unpacked_no_duplicates = 0;
- index_type nnz_no_duplicates = nnz;
-
-
- //maybe check if they exist?
- for (index_type i = 0; i < nnz; ++i)
- {
- nnz_counts[coo_rows[i]]++;
- nnz_unpacked++;
- nnz_unpacked_no_duplicates++;
- if ((coo_rows[i] != coo_cols[i])&& (symmetric) && (!expanded))
- {
- nnz_counts[coo_cols[i]]++;
- nnz_unpacked++;
- nnz_unpacked_no_duplicates++;
- }
- if (coo_rows[i] == coo_cols[i]){
- if (diag_control[coo_rows[i]] > 0) {
- //duplicate
- nnz_unpacked_no_duplicates--;
- nnz_no_duplicates--;
- }
- diag_control[coo_rows[i]]++;
- }
- }
- A_csr->setExpanded(true);
- A_csr->setNnzExpanded(nnz_unpacked_no_duplicates);
- index_type* csr_ia = new index_type[n+1];
- std::fill_n(csr_ia, n + 1, 0);
- index_type* csr_ja = new index_type[nnz_unpacked];
- real_type* csr_a = new real_type[nnz_unpacked];
- index_type* nnz_shifts = new index_type[n];
- std::fill_n(nnz_shifts, n , 0);
-
- IndexValuePair* tmp = new IndexValuePair[nnz_unpacked]; 
-
- csr_ia[0] = 0;
-
- for (index_type i = 1; i < n + 1; ++i){
- csr_ia[i] = csr_ia[i - 1] + nnz_counts[i - 1] - (diag_control[i-1] - 1);
- }
-
- int r, start;
-
-
- for (index_type i = 0; i < nnz; ++i){
- //which row
- r = coo_rows[i];
- start = csr_ia[r];
-
- if ((start + nnz_shifts[r]) > nnz_unpacked) {
- out::warning() << "index out of bounds (case 1) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl;
- }
- if ((r == coo_cols[i]) && (diag_control[r] > 1)) {//diagonal, and there are duplicates
- bool already_there = false; 
- for (index_type j = start; j < start + nnz_shifts[r]; ++j)
- {
- index_type c = tmp[j].getIdx();
- if (c == r) {
- real_type val = tmp[j].getValue();
- val += coo_vals[i];
- tmp[j].setValue(val);
- already_there = true;
- out::warning() << " duplicate found, row " << c << " adding in place " << j << " current value: " << val << std::endl;
- } 
- } 
- if (!already_there){ // first time this duplicates appears
-
- tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]);
- tmp[start + nnz_shifts[r]].setValue(coo_vals[i]);
-
- nnz_shifts[r]++;
- }
- } else {//not diagonal
- tmp[start + nnz_shifts[r]].setIdx(coo_cols[i]);
- tmp[start + nnz_shifts[r]].setValue(coo_vals[i]);
- nnz_shifts[r]++;
-
- if ((coo_rows[i] != coo_cols[i]) && (symmetric == 1))
- {
- r = coo_cols[i];
- start = csr_ia[r];
-
- if ((start + nnz_shifts[r]) > nnz_unpacked)
- out::warning() << "index out of bounds (case 2) start: " << start << "nnz_shifts[" << r << "] = " << nnz_shifts[r] << std::endl;
- tmp[start + nnz_shifts[r]].setIdx(coo_rows[i]);
- tmp[start + nnz_shifts[r]].setValue(coo_vals[i]);
- nnz_shifts[r]++;
- }
- }
- }
- //now sort whatever is inside rows
-
- for (int i = 0; i < n; ++i)
- {
-
- //now sorting (and adding 1)
- int colStart = csr_ia[i];
- int colEnd = csr_ia[i + 1];
- int length = colEnd - colStart;
- std::sort(&tmp[colStart],&tmp[colStart] + length);
- }
-
- for (index_type i = 0; i < nnz_unpacked; ++i)
- {
- csr_ja[i] = tmp[i].getIdx();
- csr_a[i] = tmp[i].getValue();
- }
-#if 0
- for (int i = 0; i<n; ++i){
- printf("Row: %d \n", i);
- for (int j = csr_ia[i]; j<csr_ia[i+1]; ++j){
- printf("(%d %16.16f) ", csr_ja[j], csr_a[j]);
- }
- printf("\n");
- }
-#endif
- A_csr->setNnz(nnz_no_duplicates);
- A_csr->updateData(csr_ia, csr_ja, csr_a, memory::HOST, memspace);
-
- delete [] nnz_counts;
- delete [] tmp;
- delete [] nnz_shifts;
- delete [] csr_ia;
- delete [] csr_ja;
- delete [] csr_a;
- delete [] diag_control; 
-
- return 0;
+ return matrix::coo2csr(A_coo, A_csr, memspace);
  }
 
  /**