From 3c1a9263cbb103c1f6ff4bf7fd50cf449db0f8e0 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Tue, 8 Nov 2011 16:55:51 +0100
Subject: [PATCH] implementation of method eraseEntries() in class template
CRSMatrix - backup of work
---
MathLib/LinAlg/Sparse/CRSMatrix.h | 206 +++++++++++++++++-
MathLib/LinAlg/Sparse/CRSTranspose.h | 8 +-
MathLib/sparse.h | 71 +++---
SimpleTests/MatrixTests/CMakeLists.txt | 16 ++
.../MatrixTests/MatTestRemoveRowsCols.cpp | 71 ++++++
5 files changed, 322 insertions(+), 50 deletions(-)
create mode 100644 SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
diff --git a/MathLib/LinAlg/Sparse/CRSMatrix.h b/MathLib/LinAlg/Sparse/CRSMatrix.h
index ba7042f5122..49737ebe79e 100644
--- a/MathLib/LinAlg/Sparse/CRSMatrix.h
+++ b/MathLib/LinAlg/Sparse/CRSMatrix.h
@@ -12,6 +12,11 @@
#include <fstream>
#include <iostream>
#include <cassert>
+
+// Base
+#include "swap.h"
+
+// MathLib
#include "SparseMatrixBase.h"
#include "sparse.h"
#include "amuxCRS.h"
@@ -113,23 +118,206 @@ public:
}
}
- void eraseEntries (IDX_TYPE n_entries, IDX_TYPE* rows, IDX_TYPE *cols)
+ /**
+ * get const access to the row pointer array of CRS matrix
+ * @return the index array _row_ptr
+ */
+ IDX_TYPE const* getRowPtrArray() const { return _row_ptr; }
+
+ /**
+ * get const access to the column index array of CRS matrix
+ * @return the index array _col_idx
+ */
+ IDX_TYPE const* getColIdxArray ()const { return _col_idx; }
+
+ /**
+ * get the matrix entries within an array of CRS matrix
+ * @return
+ */
+ FP_TYPE const* getEntryArray() const { return _data; }
+
+ /**
+ * erase rows and columns from sparse matrix
+ * @param n_rows_cols number of rows / columns to remove
+ * @param rows sorted list of rows that should be removed
+ * @param cols sorted list of columns that should be removed
+ */
+ void eraseEntries (IDX_TYPE n_rows_cols,
+ IDX_TYPE const*const rows, IDX_TYPE const*const cols)
{
- // ToDo
- //*** determine the number of new entries
- unsigned cnt_new_entries(_row_ptr[MatrixBase::_n_rows]);
+ IDX_TYPE n_cols(MatrixBase::_n_rows);
+ //*** remove the rows
+ removeRows(n_rows_cols, rows);
+ //*** transpose
+ transpose(n_cols);
+ //*** remove columns in original means removing rows in the transposed
+ removeRows(n_rows_cols, cols);
+ //*** transpose again
+ transpose(MatrixBase::_n_rows);
+ //*** apply changes to the column indices
+ correctColumnIndices(n_rows_cols, cols);
+ }
+
+ /**
+ * get the j-th column of the sparse matrix
+ * @param j the column number that should be returned
+ * @param column_entries the column entries (have to be allocated
+ */
+ void getColumn(IDX_TYPE j, FP_TYPE* column_entries) const
+ {
+ for (IDX_TYPE k(0); k<MatrixBase::_n_rows; k++) {
+ const IDX_TYPE end_row(_row_ptr[k+1]);
+ IDX_TYPE i(_row_ptr[k+1]);
+ while (i<end_row && _col_idx[i] != j) {
+ i++;
+ }
+ if (i==end_row) {
+ column_entries[k] = 0.0;
+ } else {
+ column_entries[k] = _data[i];
+ }
+ }
+ }
+
+protected:
+ void removeRows (IDX_TYPE n_rows_cols, IDX_TYPE const*const rows)
+ {
+ //*** determine the number of new rows and the number of entries without the rows
+ const IDX_TYPE n_new_rows(MatrixBase::_n_rows - n_rows_cols);
+ IDX_TYPE *row_ptr_new(new IDX_TYPE[n_new_rows+1]);
+ row_ptr_new[0] = 0;
+ IDX_TYPE row_cnt (1), erase_row_cnt(0);
+ for (unsigned k(0); k<MatrixBase::_n_rows; k++) {
+ if (k != rows[erase_row_cnt]) {
+ row_ptr_new[row_cnt] = _row_ptr[k+1] - _row_ptr[k];
+ row_cnt++;
+ } else {
+ erase_row_cnt++;
+ }
+ }
+
+ //*** sum up the entries
+ for (IDX_TYPE k(0); k<n_new_rows; k++) {
+ row_ptr_new[k+1] = row_ptr_new[k+1] + row_ptr_new[k];
+ }
+
+ //*** create new memory for col_idx and data
+ IDX_TYPE nnz_new(row_ptr_new[n_new_rows]);
+ IDX_TYPE *col_idx_new (new IDX_TYPE[nnz_new]);
+ FP_TYPE *data_new (new FP_TYPE[nnz_new]);
- //*** create new memory
//*** copy the entries
- //*** apply changes to attributes
+ // initialization
+ IDX_TYPE *row_ptr_new_tmp(new IDX_TYPE[n_new_rows+1]);
+ for (unsigned k(0); k<=n_new_rows; k++) {
+ row_ptr_new_tmp[k] = row_ptr_new[k];
+ }
+ erase_row_cnt = 0;
+ row_cnt = 0;
+ // copy column index and data entries
+ for (IDX_TYPE k(0); k<MatrixBase::_n_rows; k++) {
+ if (k != rows[erase_row_cnt]) {
+ const IDX_TYPE end (_row_ptr[k+1]);
+ // walk through row
+ for (IDX_TYPE j(_row_ptr[k]); j<end; j++) {
+ col_idx_new[row_ptr_new_tmp[row_cnt]] = _col_idx[j];
+ data_new[row_ptr_new_tmp[row_cnt]] = _data[j];
+ row_ptr_new_tmp[row_cnt]++;
+ }
+ row_cnt++;
+ } else {
+ erase_row_cnt++;
+ }
+ }
+
+ MatrixBase::_n_rows -= n_rows_cols;
+ BASELIB::swap (row_ptr_new, _row_ptr);
+ BASELIB::swap (col_idx_new, _col_idx);
+ BASELIB::swap (data_new, _data);
+
+ delete [] row_ptr_new;
+ delete [] col_idx_new;
+ delete [] data_new;
+ }
+
+ void transpose (IDX_TYPE n_cols)
+ {
+ // create a helper array row_ptr_nnz
+ IDX_TYPE *row_ptr_nnz(new IDX_TYPE[n_cols+1]);
+ for (IDX_TYPE k(0); k <= n_cols; k++) {
+ row_ptr_nnz[k] = 0;
+ }
+
+ // count entries per row in the transposed matrix
+ IDX_TYPE nnz(_row_ptr[MatrixBase::_n_rows]);
+ for (IDX_TYPE k(0); k < nnz; k++) {
+ row_ptr_nnz[_col_idx[k]]++;
+ }
+
+ // create row_ptr_trans
+ IDX_TYPE *row_ptr_trans(new IDX_TYPE[n_cols + 1]);
+ row_ptr_trans[0] = 0;
+ for (IDX_TYPE k(0); k < n_cols; k++) {
+ row_ptr_trans[k+1] = row_ptr_trans[k] + row_ptr_nnz[k];
+ }
+
+ // make a copy of row_ptr_trans
+ for (IDX_TYPE k(0); k <= n_cols; k++) {
+ row_ptr_nnz[k] = row_ptr_trans[k];
+ }
+
+ // create arrays col_idx_trans and data_trans
+ assert(nnz == row_ptr_trans[n_cols]);
+ IDX_TYPE *col_idx_trans(new IDX_TYPE[nnz]);
+ FP_TYPE *data_trans(new FP_TYPE[nnz]);
+
+ // fill arrays col_idx_trans and data_trans
+ for (IDX_TYPE i(0); i < MatrixBase::_n_rows; i++) {
+ const IDX_TYPE row_end(_row_ptr[i + 1]);
+ for (IDX_TYPE j(_row_ptr[i]); j < row_end; j++) {
+ const IDX_TYPE k(_col_idx[j]);
+ col_idx_trans[row_ptr_nnz[k]] = i;
+ data_trans[row_ptr_nnz[k]] = _data[j];
+ row_ptr_nnz[k]++;
+ }
+ }
+
+ MatrixBase::_n_rows = n_cols;
+ BASELIB::swap(row_ptr_trans, _row_ptr);
+ BASELIB::swap(col_idx_trans, _col_idx);
+ BASELIB::swap(data_trans, _data);
+
+ delete[] row_ptr_nnz;
+ delete[] row_ptr_trans;
+ delete[] col_idx_trans;
+ delete[] data_trans;
}
- void getColumn(IDX_TYPE col, FP_TYPE* column_entries)
+ void correctColumnIndices(IDX_TYPE n_rows_cols, IDX_TYPE const*const cols)
{
- // ToDo
+ // create the mapping
+ const IDX_TYPE size(MatrixBase::_n_rows + n_rows_cols);
+ IDX_TYPE *mapping(new IDX_TYPE[size]);
+ IDX_TYPE k(0);
+ for (IDX_TYPE i(0); i < n_rows_cols; i++) {
+ while (k < cols[i]) {
+ mapping[k] = k - i;
+ k++;
+ }
+ }
+ for (IDX_TYPE k(cols[n_rows_cols - 1]); k < size; k++) {
+ mapping[k] = k - n_rows_cols;
+ }
+
+ // apply mapping to _col_idx
+ const IDX_TYPE nnz(_row_ptr[MatrixBase::_n_rows]);
+ for (IDX_TYPE k(0); k < nnz; k++) {
+ _col_idx[k] = mapping[_col_idx[k]];
+ }
+
+ delete[] mapping;
}
-protected:
IDX_TYPE *_row_ptr;
IDX_TYPE *_col_idx;
FP_TYPE* _data;
diff --git a/MathLib/LinAlg/Sparse/CRSTranspose.h b/MathLib/LinAlg/Sparse/CRSTranspose.h
index 6b3eba98604..010a49e5a31 100644
--- a/MathLib/LinAlg/Sparse/CRSTranspose.h
+++ b/MathLib/LinAlg/Sparse/CRSTranspose.h
@@ -2,17 +2,17 @@
* CRSTranspose.h
*
* Created on: Sep 27, 2011
- * Author: fischeth
+ * Author: TF
*/
#ifndef CRSTRANSPOSE_H_
#define CRSTRANSPOSE_H_
-void CS_transp(unsigned n, double* A, unsigned* jA, unsigned *iA,
- double* B, unsigned *jB, unsigned *iB)
+void CS_transp(unsigned n, unsigned *iA, unsigned* jA, double* A,
+ unsigned *iB, unsigned *jB, double* B)
{
unsigned nnz = iA[n];
- unsigned *inz = new unsigned[n];
+ unsigned *inz(new unsigned[n]);
for (unsigned i=0; i<n; ++i) inz[i] = 0;
diff --git a/MathLib/sparse.h b/MathLib/sparse.h
index fa3f6623659..3278ea93c69 100644
--- a/MathLib/sparse.h
+++ b/MathLib/sparse.h
@@ -4,51 +4,48 @@
#include <iostream>
#include <cassert>
-extern void CS_write(char*, unsigned, unsigned*, unsigned*, double*);
+extern void CS_write(char*, unsigned, unsigned const*, unsigned const*, double const*);
extern void CS_read(char*, unsigned&, unsigned*&, unsigned*&, double*&);
-template<class T> void
-CS_write(std::ostream &os, unsigned n, unsigned* iA, unsigned* jA, T* A)
+template<class T> void CS_write(std::ostream &os, unsigned n, unsigned const* iA, unsigned const* jA, T const* A)
{
- os.write((char*) &n, sizeof(unsigned));
- os.write((char*) iA, (n+1)*sizeof(unsigned));
- os.write((char*) jA, iA[n]*sizeof(unsigned));
- os.write((char*) A, iA[n]*sizeof(T));
+ os.write((char*) &n, sizeof(unsigned));
+ os.write((char*) iA, (n + 1) * sizeof(unsigned));
+ os.write((char*) jA, iA[n] * sizeof(unsigned));
+ os.write((char*) A, iA[n] * sizeof(T));
}
-template<class T> void
-CS_read(std::istream &is, unsigned &n, unsigned* &iA, unsigned* &jA, T* &A)
+template<class T> void CS_read(std::istream &is, unsigned &n, unsigned* &iA, unsigned* &jA, T* &A)
{
- is.read((char*) &n, sizeof(unsigned));
- if (iA != NULL) {
- delete [] iA;
- delete [] jA;
- delete [] A;
- }
- iA = new unsigned[n+1];
- assert(iA!=NULL);
- is.read((char*) iA, (n+1)*sizeof(unsigned));
-
- jA = new unsigned[iA[n]];
- assert(jA!=NULL);
- is.read((char*) jA, iA[n]*sizeof(unsigned));
-
- A = new T[iA[n]];
- assert(A!=NULL);
- is.read((char*) A, iA[n]*sizeof(T));
+ is.read((char*) &n, sizeof(unsigned));
+ if (iA != NULL) {
+ delete[] iA;
+ delete[] jA;
+ delete[] A;
+ }
+ iA = new unsigned[n + 1];
+ assert(iA != NULL);
+ is.read((char*) iA, (n + 1) * sizeof(unsigned));
+
+ jA = new unsigned[iA[n]];
+ assert(jA != NULL);
+ is.read((char*) jA, iA[n] * sizeof(unsigned));
+
+ A = new T[iA[n]];
+ assert(A != NULL);
+ is.read((char*) A, iA[n] * sizeof(T));
#ifndef NDEBUG
- // do simple checks
- if (iA[0]!=0)
- std::cerr << std::endl
- << "CRS matrix: array iA doesn't start with 0" << std::endl;
-
- unsigned i = 0;
- while (i<iA[n] && jA[i]<n) ++i;
- if (i<iA[n])
- std::cerr << std::endl
- << "CRS matrix: the " << i << "th entry of jA has the value "
- << jA[i] << ", which is out of bounds." << std::endl;
+ // do simple checks
+ if (iA[0] != 0) std::cerr << std::endl << "CRS matrix: array iA doesn't start with 0"
+ << std::endl;
+
+ unsigned i = 0;
+ while (i < iA[n] && jA[i] < n)
+ ++i;
+ if (i < iA[n]) std::cerr << std::endl << "CRS matrix: the " << i
+ << "th entry of jA has the value " << jA[i] << ", which is out of bounds."
+ << std::endl;
#endif
}
diff --git a/SimpleTests/MatrixTests/CMakeLists.txt b/SimpleTests/MatrixTests/CMakeLists.txt
index bda3c0ca41a..50238c4a25c 100644
--- a/SimpleTests/MatrixTests/CMakeLists.txt
+++ b/SimpleTests/MatrixTests/CMakeLists.txt
@@ -19,6 +19,13 @@ ADD_EXECUTABLE( MatMult
${HEADERS}
)
+# Create the executable
+ADD_EXECUTABLE( MatTestRemoveRowsCols
+ MatTestRemoveRowsCols.cpp
+ ${SOURCES}
+ ${HEADERS}
+)
+
IF (WIN32)
TARGET_LINK_LIBRARIES(MatMult Winmm.lib)
ENDIF (WIN32)
@@ -28,3 +35,12 @@ TARGET_LINK_LIBRARIES ( MatMult
MathLib
)
+IF (WIN32)
+ TARGET_LINK_LIBRARIES(MatTestRemoveRowsCols Winmm.lib)
+ENDIF (WIN32)
+
+TARGET_LINK_LIBRARIES ( MatTestRemoveRowsCols
+ Base
+ MathLib
+)
+
diff --git a/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp b/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
new file mode 100644
index 00000000000..905f1498bf3
--- /dev/null
+++ b/SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp
@@ -0,0 +1,71 @@
+/*
+ * MatTestRemoveRowsCols.cpp
+ *
+ * Created on: Nov 8, 2011
+ * Author: TF
+ */
+
+#include <fstream>
+#include <iostream>
+
+// Base
+#include "RunTimeTimer.h"
+#include "CPUTimeTimer.h"
+
+// MathLib
+#include "LinAlg/Sparse/CRSMatrix.h"
+
+int main(int argc, char *argv[])
+{
+ if (argc < 3) {
+ std::cout << "Usage: " << argv[0] << " input-matrix output-matrix" << std::endl;
+ return 1;
+ }
+
+ std::string fname_mat (argv[1]);
+ bool verbose (true);
+
+ // *** reading matrix in crs format from file
+ std::ifstream in(fname_mat.c_str(), std::ios::in | std::ios::binary);
+ double *A(NULL);
+ unsigned *iA(NULL), *jA(NULL), n;
+ if (in) {
+ if (verbose) {
+ std::cout << "reading matrix from " << fname_mat << " ... " << std::flush;
+ }
+ RunTimeTimer timer;
+ timer.start();
+ CS_read(in, n, iA, jA, A);
+ in.close();
+ timer.stop();
+ if (verbose) {
+ std::cout << "ok, " << timer.elapsed() << " s)" << std::endl;
+ }
+ } else {
+ std::cout << "error reading matrix from " << fname_mat << std::endl;
+ }
+ unsigned nnz(iA[n]);
+ if (verbose) {
+ std::cout << "Parameters read: n=" << n << ", nnz=" << nnz << std::endl;
+ }
+
+ MathLib::CRSMatrix<double, unsigned> mat (n, iA, jA, A);
+
+ const unsigned n_rows_cols_to_erase(10);
+ unsigned *rows_to_erase(new unsigned[n_rows_cols_to_erase]);
+ unsigned *cols_to_erase(new unsigned[n_rows_cols_to_erase]);
+
+ for (unsigned k(0); k<n_rows_cols_to_erase; k++) {
+ rows_to_erase[k] = k*5;
+ cols_to_erase[k] = k*5;
+ }
+
+ mat.eraseEntries(n_rows_cols_to_erase, rows_to_erase, cols_to_erase);
+ delete[] rows_to_erase;
+ delete[] cols_to_erase;
+
+ fname_mat = argv[2];
+ std::ofstream out (fname_mat.c_str(), std::ios::binary);
+ CS_write (out, mat.getNRows(), mat.getRowPtrArray(), mat.getColIdxArray(), mat.getEntryArray());
+ out.close();
+}
--
GitLab