implementation of method eraseEntries() in class template CRSMatrix - backup of work

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;

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;
 

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
 }
 

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
+)
+

SimpleTests/MatrixTests/MatTestRemoveRowsCols.cpp

+/*
+ * 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();
+}