diff --git a/Source/C/SparseMatrix_c.h b/Source/C/SparseMatrix_c.h index 4d3ec919..1e397c8e 100644 --- a/Source/C/SparseMatrix_c.h +++ b/Source/C/SparseMatrix_c.h @@ -1,8 +1,6 @@ #ifndef SparseMatrix_ch #define SparseMatrix_ch -// -// void ConstructEmptySparseMatrix_wrp(int *ih_this, const int *columns, -// const int *rows); + void ConstructSparseMatrixFromFile_wrp(int *ih_this, const char *file_name, const int *name_size); void ConstructFromTripletList_wrp(int *ih_this, const int *ih_triplet_list, @@ -12,6 +10,9 @@ void ConstructZeroSparseMatrix_wrp(int *ih_this, const int *rows, void CopySparseMatrix_wrp(const int *ih_matA, int *ih_matB); void GetRows_wrp(const int *ih_this, int *rows); void GetColumns_wrp(const int *ih_this, int *columns); +void ExtractRow_wrp(const int *ih_this, int *row_number, int *ih_row_out); +void ExtractColumn_wrp(const int *ih_this, int *column_number, + int *ih_column_out); void ScaleSparseMatrix_wrp(int *ih_this, const double *constant); void IncrementSparseMatrix_wrp(const int *ih_matA, int *ih_matB, const double *alpha_in, diff --git a/Source/CPlusPlus/SparseMatrix.cc b/Source/CPlusPlus/SparseMatrix.cc index ce5fa01f..9541e030 100644 --- a/Source/CPlusPlus/SparseMatrix.cc +++ b/Source/CPlusPlus/SparseMatrix.cc @@ -49,6 +49,19 @@ int SparseMatrix::GetColumns() const { return temp; } +//////////////////////////////////////////////////////////////////////////////// +void SparseMatrix::ExtractRow(int row_number, SparseMatrix &row_out) const { + int temp = row_number + 1; + ExtractRow_wrp(ih_this, &temp, row_out.ih_this); +} + +//////////////////////////////////////////////////////////////////////////////// +void SparseMatrix::ExtractColumn(int column_number, + SparseMatrix &column_out) const { + int temp = column_number + 1; + ExtractColumn_wrp(ih_this, &temp, column_out.ih_this); +} + //////////////////////////////////////////////////////////////////////////////// void SparseMatrix::Scale(double constant) { ScaleSparseMatrix_wrp(ih_this, &constant); diff --git a/Source/CPlusPlus/SparseMatrix.h b/Source/CPlusPlus/SparseMatrix.h index bfeb3202..177edcb5 100644 --- a/Source/CPlusPlus/SparseMatrix.h +++ b/Source/CPlusPlus/SparseMatrix.h @@ -33,6 +33,16 @@ class SparseMatrix { int GetRows() const; //! Get the number of columns in a matrix. int GetColumns() const; + //! Extract a row from the matrix into the compressed vector representation. + //!\param this the matrix to extract from. + //!\param row_number the row to extract + //!\param row_out the matrix representing that row + void ExtractRow(int row_number, SparseMatrix& row_out) const; + //! Extract a column from the matrix into the compressed vector representation. + //!\param this the matrix to extract from. + //!\param column_number the column to extract + //!\param column_out the matrix representing that column + void ExtractColumn(int column_number, SparseMatrix& column_out) const; public: //! Scale the matrix by a constant. diff --git a/Source/Fortran/SparseMatrixModule.f90 b/Source/Fortran/SparseMatrixModule.f90 index 56ba8e44..0d857de9 100644 --- a/Source/Fortran/SparseMatrixModule.f90 +++ b/Source/Fortran/SparseMatrixModule.f90 @@ -31,6 +31,8 @@ MODULE SparseMatrixModule !! Basic Accessors PUBLIC :: GetRows PUBLIC :: GetColumns + PUBLIC :: ExtractRow + PUBLIC :: ExtractColumn !! Helper routines PUBLIC :: SplitSparseMatrixColumns PUBLIC :: ComposeSparseMatrixColumns @@ -46,8 +48,8 @@ MODULE SparseMatrixModule !! @param[in] rows number of matrix rows. PURE SUBROUTINE ConstructEmptySparseMatrix(this, columns, rows) !! Parameters - TYPE(SparseMatrix_t),INTENT(out) :: this - INTEGER, INTENT(in) :: columns, rows + TYPE(SparseMatrix_t),INTENT(OUT) :: this + INTEGER, INTENT(IN) :: columns, rows CALL DestructSparseMatrix(this) this%rows = rows @@ -61,8 +63,8 @@ END SUBROUTINE ConstructEmptySparseMatrix !! @param[in] file_name name of the file. SUBROUTINE ConstructSparseMatrixFromFile(this, file_name) !! Parameters - TYPE(SparseMatrix_t), INTENT(out) :: this - CHARACTER(len=*), INTENT(in) :: file_name + TYPE(SparseMatrix_t), INTENT(OUT) :: this + CHARACTER(len=*), INTENT(IN) :: file_name !! About the matrix market file. INTEGER :: sparsity_type, data_type, pattern_type !! Local Data @@ -121,8 +123,8 @@ END SUBROUTINE ConstructSparseMatrixFromFile !! @param[in] columns number of matrix columns PURE SUBROUTINE ConstructZeroSparseMatrix(this,rows,columns) !! Parameters - TYPE(SparseMatrix_t), INTENT(out) :: this - INTEGER, INTENT(in) :: rows, columns + TYPE(SparseMatrix_t), INTENT(OUT) :: this + INTEGER, INTENT(IN) :: rows, columns !! Allocate CALL ConstructEmptySparseMatrix(this,columns,rows) @@ -140,9 +142,9 @@ END SUBROUTINE ConstructZeroSparseMatrix !! @param[in] columns number of matrix columns PURE SUBROUTINE ConstructFromTripletList(this,triplet_list,rows,columns) !! Parameters - TYPE(SparseMatrix_t), INTENT(out) :: this - TYPE(TripletList_t), INTENT(in) :: triplet_list - INTEGER, INTENT(in) :: rows, columns + TYPE(SparseMatrix_t), INTENT(OUT) :: this + TYPE(TripletList_t), INTENT(IN) :: triplet_list + INTEGER, INTENT(IN) :: rows, columns INTEGER :: outer_array_ptr INTEGER :: values_counter @@ -187,7 +189,7 @@ END SUBROUTINE ConstructFromTripletList !! @param[inout] this the matrix to free up PURE SUBROUTINE DestructSparseMatrix(this) !! Parameters - TYPE(SparseMatrix_t), INTENT(inout) :: this + TYPE(SparseMatrix_t), INTENT(INOUT) :: this IF (ALLOCATED(this%outer_index)) THEN DEALLOCATE(this%outer_index) @@ -205,8 +207,8 @@ END SUBROUTINE DestructSparseMatrix !! @param[inout] matB = matA PURE SUBROUTINE CopySparseMatrix(matA, matB) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: matA - TYPE(SparseMatrix_t), INTENT(inout) :: matB + TYPE(SparseMatrix_t), INTENT(IN) :: matA + TYPE(SparseMatrix_t), INTENT(INOUT) :: matB CALL DestructSparseMatrix(matB) matB = matA @@ -217,7 +219,7 @@ END SUBROUTINE CopySparseMatrix !! @result number of rows. PURE FUNCTION GetRows(this) RESULT(rows) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: this + TYPE(SparseMatrix_t), INTENT(IN) :: this INTEGER :: rows rows = this%rows END FUNCTION GetRows @@ -227,10 +229,63 @@ END FUNCTION GetRows !! @result number of columns. PURE FUNCTION GetColumns(this) RESULT(columns) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: this + TYPE(SparseMatrix_t), INTENT(IN) :: this INTEGER :: columns columns = this%columns END FUNCTION GetColumns +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Extract a row from the matrix into the compressed vector representation. + !! @param[in] this the matrix to extract from. + !! @param[in] row_number the row to extract + !! @param[out] row_out the matrix representing that row + PURE SUBROUTINE ExtractRow(this, row_number, row_out) + !! Parameters + TYPE(SparseMatrix_t), INTENT(IN) :: this + INTEGER, INTENT(IN) :: row_number + TYPE(SparseMatrix_t), INTENT(INOUT) :: row_out + !! Local variables + INTEGER :: number_of_values + INTEGER :: start_index + INTEGER :: counter + TYPE(SparseMatrix_t) :: temp, temp_c + + CALL TransposeSparseMatrix(this,temp) + CALL ExtractColumn(temp, row_number, temp_c) + CALL TransposeSparseMatrix(temp_c,row_out) + CALL DestructSparseMatrix(temp_c) + CALL DestructSparseMatrix(temp) + END SUBROUTINE ExtractRow +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Extract a column from the matrix into the compressed vector representation. + !! @param[in] this the matrix to extract from. + !! @param[in] column_number the row to extract + !! @param[out] column_out the matrix representing that row + PURE SUBROUTINE ExtractColumn(this, column_number, column_out) + !! Parameters + TYPE(SparseMatrix_t), INTENT(IN) :: this + INTEGER, INTENT(IN) :: column_number + TYPE(SparseMatrix_t), INTENT(INOUT) :: column_out + !! Local variables + INTEGER :: number_of_values + INTEGER :: start_index + INTEGER :: counter + + !! Allocate Memory + CALL ConstructEmptySparseMatrix(column_out, 1, this%rows) + start_index = this%outer_index(column_number) + number_of_values = this%outer_index(column_number+1) - & + & this%outer_index(column_number) + ALLOCATE(column_out%inner_index(number_of_values)) + ALLOCATE(column_out%values(number_of_values)) + + !! Copy Values + column_out%outer_index(1) = 0 + column_out%outer_index(2) = number_of_values + DO counter=1, number_of_values + column_out%inner_index(counter) = this%inner_index(start_index+counter) + column_out%values(counter) = this%values(start_index+counter) + END DO + END SUBROUTINE ExtractColumn !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Transpose a sparse matrix and return it in a separate matrix. !! The current implementation has you go from matrix to triplet list, @@ -240,8 +295,8 @@ END FUNCTION GetColumns !! @param[out] matT the input matrix transposed. PURE SUBROUTINE TransposeSparseMatrix(this, matT) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: this - TYPE(SparseMatrix_t), INTENT(inout) :: matT + TYPE(SparseMatrix_t), INTENT(IN) :: this + TYPE(SparseMatrix_t), INTENT(INOUT) :: matT !! Local Data TYPE(TripletList_t) :: triplet_list TYPE(TripletList_t) :: sorted_triplet_list @@ -285,8 +340,8 @@ END SUBROUTINE TransposeSparseMatrix !! @param[out] out_matrix = [Matrix 1 | Matrix 2, ...] . PURE SUBROUTINE ComposeSparseMatrixColumns(mat_list, out_matrix) !! Parameters - TYPE(SparseMatrix_t), DIMENSION(:), INTENT(in) :: mat_list - TYPE(SparseMatrix_t), INTENT(inout) :: out_matrix + TYPE(SparseMatrix_t), DIMENSION(:), INTENT(IN) :: mat_list + TYPE(SparseMatrix_t), INTENT(INOUT) :: out_matrix !! Local Variables INTEGER :: total_columns, total_values INTEGER :: inner_start, inner_length @@ -342,10 +397,10 @@ END SUBROUTINE ComposeSparseMatrixColumns PURE SUBROUTINE SplitSparseMatrixColumns(this, num_blocks, split_list, & & block_offsets_out) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: this - INTEGER, INTENT(in) :: num_blocks - TYPE(SparseMatrix_t), DIMENSION(num_blocks), INTENT(out) :: split_list - INTEGER, DIMENSION(num_blocks+1), INTENT(out), OPTIONAL :: block_offsets_out + TYPE(SparseMatrix_t), INTENT(IN) :: this + INTEGER, INTENT(IN) :: num_blocks + TYPE(SparseMatrix_t), DIMENSION(num_blocks), INTENT(OUT) :: split_list + INTEGER, DIMENSION(num_blocks+1), INTENT(OUT), OPTIONAL :: block_offsets_out !! Local Data INTEGER, DIMENSION(num_blocks) :: block_sizes INTEGER, DIMENSION(num_blocks+1) :: block_offsets @@ -391,12 +446,12 @@ PURE SUBROUTINE SplitSparseMatrixColumns(this, num_blocks, split_list, & total_values = split_list(split_counter)%outer_index(lcolumns+1) !! Copy Inner Indices and Values IF (total_values .GT. 0) THEN - ALLOCATE(split_list(split_counter)%inner_index(total_values)) - split_list(split_counter)%inner_index = & - & this%inner_index(linner_offset:linner_offset+total_values-1) - ALLOCATE(split_list(split_counter)%values(total_values)) - split_list(split_counter)%values = & - & this%values(linner_offset:linner_offset+total_values-1) + ALLOCATE(split_list(split_counter)%inner_index(total_values)) + split_list(split_counter)%inner_index = & + & this%inner_index(linner_offset:linner_offset+total_values-1) + ALLOCATE(split_list(split_counter)%values(total_values)) + split_list(split_counter)%values = & + & this%values(linner_offset:linner_offset+total_values-1) END IF END DO END IF @@ -410,8 +465,8 @@ END SUBROUTINE SplitSparseMatrixColumns !! @param[out] triplet_list the triplet list we created. PURE SUBROUTINE MatrixToTripletList(this, triplet_list) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: this - TYPE(TripletList_t), INTENT(inout) :: triplet_list + TYPE(SparseMatrix_t), INTENT(IN) :: this + TYPE(TripletList_t), INTENT(INOUT) :: triplet_list !! Helper variables INTEGER :: outer_counter, inner_counter INTEGER :: elements_per_inner @@ -441,8 +496,8 @@ END SUBROUTINE MatrixToTripletList !! @param[in] file_name_in optionally, you can pass a file to print to. SUBROUTINE PrintSparseMatrix(this, file_name_in) !! Parameters - TYPE(SparseMatrix_t), INTENT(in) :: this - CHARACTER(len=*), OPTIONAL, INTENT(in) :: file_name_in + TYPE(SparseMatrix_t), INTENT(IN) :: this + CHARACTER(len=*), OPTIONAL, INTENT(IN) :: file_name_in !! Local Data TYPE(TripletList_t) :: triplet_list INTEGER :: file_handler diff --git a/Source/Wrapper/SparseMatrixModule_wrp.f90 b/Source/Wrapper/SparseMatrixModule_wrp.f90 index 77bb278a..7341db25 100644 --- a/Source/Wrapper/SparseMatrixModule_wrp.f90 +++ b/Source/Wrapper/SparseMatrixModule_wrp.f90 @@ -7,7 +7,8 @@ MODULE SparseMatrixModule_wrp & ConstructEmptySparseMatrix, ConstructSparseMatrixFromFile, & & ConstructFromTripletList, DestructSparseMatrix, CopySparseMatrix, & & TransposeSparseMatrix, PrintSparseMatrix, MatrixToTripletList, & - & GetRows, GetColumns, ComposeSparseMatrixColumns + & GetRows, GetColumns, ComposeSparseMatrixColumns, & + & ExtractRow, ExtractColumn USE TripletListModule_wrp, ONLY : TripletList_wrp USE WrapperModule, ONLY : SIZE_wrp USE ISO_C_BINDING, ONLY : c_int, c_char, c_bool @@ -26,6 +27,8 @@ MODULE SparseMatrixModule_wrp PUBLIC :: DestructSparseMatrix_wrp PUBLIC :: GetRows_wrp PUBLIC :: GetColumns_wrp + PUBLIC :: ExtractRow_wrp + PUBLIC :: ExtractColumn_wrp PUBLIC :: CopySparseMatrix_wrp PUBLIC :: TransposeSparseMatrix_wrp PUBLIC :: PrintSparseMatrix_wrp @@ -87,7 +90,7 @@ END SUBROUTINE ConstructZeroSparseMatrix_wrp !> Explicitly destruct a sparse matrix PURE SUBROUTINE DestructSparseMatrix_wrp(ih_this) & & bind(c,name="DestructSparseMatrix_wrp") - INTEGER(kind=c_int), INTENT(inout) :: ih_this(SIZE_wrp) + INTEGER(kind=c_int), INTENT(INOUT) :: ih_this(SIZE_wrp) TYPE(SparseMatrix_wrp) :: h_this h_this = TRANSFER(ih_this,h_this) @@ -128,6 +131,44 @@ PURE SUBROUTINE GetColumns_wrp(ih_this, columns) bind(c,name="GetColumns_wrp") h_this = TRANSFER(ih_this,h_this) columns = GetColumns(h_this%data) END SUBROUTINE GetColumns_wrp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Extract a row from the matrix into the compressed vector representation. + !! @param[in] ih_this the matrix to extrat from. + !! @param[in] row_number the row to extract + !! @param[out] ih_row_out the matrix representing that row + PURE SUBROUTINE ExtractRow_wrp(ih_this, row_number, ih_row_out) & + & bind(c,name="ExtractRow_wrp") + INTEGER(kind=c_int), INTENT(IN) :: ih_this(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: row_number + INTEGER(kind=c_int), INTENT(INOUT) :: ih_row_out(SIZE_wrp) + TYPE(SparseMatrix_wrp) :: h_this + TYPE(SparseMatrix_wrp) :: h_row_out + + ALLOCATE(h_row_out%data) + h_this = TRANSFER(ih_this,h_this) + CALL ExtractRow(h_this%data, row_number, h_row_out%data) + + ih_row_out= TRANSFER(h_row_out,ih_row_out) + END SUBROUTINE ExtractRow_wrp +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !> Extract a column from the matrix into the compressed vector representation. + !! @param[in] ih_this the matrix to extrat from. + !! @param[in] column_number the row to extract. + !! @param[out] ih_column_out the matrix representing that column. + PURE SUBROUTINE ExtractColumn_wrp(ih_this, column_number, ih_column_out) & + & bind(c,name="ExtractColumn_wrp") + INTEGER(kind=c_int), INTENT(IN) :: ih_this(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: column_number + INTEGER(kind=c_int), INTENT(INOUT) :: ih_column_out(SIZE_wrp) + TYPE(SparseMatrix_wrp) :: h_this + TYPE(SparseMatrix_wrp) :: h_column_out + + ALLOCATE(h_column_out%data) + h_this = TRANSFER(ih_this,h_this) + CALL ExtractColumn(h_this%data, column_number, h_column_out%data) + + ih_column_out= TRANSFER(h_column_out, ih_column_out) + END SUBROUTINE ExtractColumn_wrp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> Wrap the matrix transpose function. PURE SUBROUTINE TransposeSparseMatrix_wrp(ih_matA, ih_matAT) & @@ -164,7 +205,7 @@ END SUBROUTINE PrintSparseMatrixF_wrp !> Warp the routine that prints the sparse matrix to the console. SUBROUTINE PrintSparseMatrix_wrp(ih_this) & & bind(c,name="PrintSparseMatrix_wrp") - INTEGER(kind=c_int), INTENT(in) :: ih_this(SIZE_wrp) + INTEGER(kind=c_int), INTENT(IN) :: ih_this(SIZE_wrp) TYPE(SparseMatrix_wrp) :: h_this h_this = TRANSFER(ih_this,h_this) diff --git a/UnitTests/testSparseMatrix.py b/UnitTests/testSparseMatrix.py index b1f44831..2b797e3d 100644 --- a/UnitTests/testSparseMatrix.py +++ b/UnitTests/testSparseMatrix.py @@ -6,12 +6,11 @@ import scipy from numpy import sum, multiply import scipy.sparse -from random import uniform +from random import uniform, randint from scipy.sparse import random, csr_matrix from scipy.sparse.linalg import norm from scipy.io import mmwrite, mmread -#import random import os from Helpers import THRESHOLD @@ -241,6 +240,46 @@ def test_multiply_zero(self): normval = abs(norm(CheckMat - ResultMat)) self.assertLessEqual(normval, THRESHOLD) + def test_get_row(self): + '''Test function that extracts a row from the matrix''' + for param in self.parameters: + if param.rows == 0: + continue + matrix1 = random(param.rows, param.columns, + param.sparsity, format="csr") + mmwrite(self.scratch_dir + "/matrix1.mtx", csr_matrix(matrix1)) + row_num = randint(0,param.rows-1) + CheckMat = matrix1[row_num, :] + + ntmatrix1 = nt.SparseMatrix(self.scratch_dir + "/matrix1.mtx") + ntmatrix2 = nt.SparseMatrix(ntmatrix1.GetColumns(), 1) + ntmatrix1.ExtractRow(row_num, ntmatrix2) + ntmatrix2.WriteToMatrixMarket(self.scratch_dir + "/matrix2.mtx") + + ResultMat = mmread(self.scratch_dir + "/matrix2.mtx") + normval = abs(norm(CheckMat - ResultMat)) + self.assertLessEqual(normval, THRESHOLD) + + def test_get_column(self): + '''Test function that extracts a column from the matrix''' + for param in self.parameters: + if param.columns == 0: + continue + matrix1 = random(param.rows, param.columns, + param.sparsity, format="csr") + mmwrite(self.scratch_dir + "/matrix1.mtx", csr_matrix(matrix1)) + column_num = randint(0,param.columns-1) + CheckMat = matrix1[:, column_num] + + ntmatrix1 = nt.SparseMatrix(self.scratch_dir + "/matrix1.mtx") + ntmatrix2 = nt.SparseMatrix(1, ntmatrix1.GetRows()) + ntmatrix1.ExtractColumn(column_num, ntmatrix2) + ntmatrix2.WriteToMatrixMarket(self.scratch_dir + "/matrix2.mtx") + + ResultMat = mmread(self.scratch_dir + "/matrix2.mtx") + normval = abs(norm(CheckMat - ResultMat)) + self.assertLessEqual(normval, THRESHOLD) + ############################################################################### if __name__ == '__main__':