diff --git a/fastlib/sparse/A.txt b/fastlib/sparse/A.txt new file mode 100755 index 0000000000..86886f5a2f --- /dev/null +++ b/fastlib/sparse/A.txt @@ -0,0 +1,40 @@ +1 1 0.31106 +13 1 0.47969 +2 2 0.734 +3 3 0.4109 +4 3 0.6412 +1 4 0.16853 +3 4 0.39979 +11 4 0.87535 +19 4 0.9667 +20 4 0.31775 +14 5 0.61591 +1 6 0.89665 +15 6 0.61663 +11 7 0.83526 +1 8 0.32272 +6 8 0.69778 +18 8 0.51015 +3 9 0.50552 +6 9 0.46189 +18 9 0.71396 +13 10 0.56082 +14 11 0.6619 +20 11 0.5877 +3 12 0.16931 +6 12 0.082613 +7 12 0.82072 +5 13 0.83685 +15 13 0.68514 +8 14 0.19302 +18 14 0.51521 +12 15 0.88071 +5 16 0.80346 +8 16 0.44535 +10 16 0.30874 +3 17 0.52475 +4 18 0.016197 +11 18 0.3331 +19 18 0.82212 +18 19 0.60587 +9 20 0.012958 diff --git a/fastlib/sparse/AdottimesB.txt b/fastlib/sparse/AdottimesB.txt new file mode 100755 index 0000000000..d7ece2c230 --- /dev/null +++ b/fastlib/sparse/AdottimesB.txt @@ -0,0 +1,6 @@ +3 3 0.37376 +19 4 0.04147 +14 5 0.1847 +1 6 0.20795 +20 11 0.43719 +18 14 0.13314 diff --git a/fastlib/sparse/AminusB.txt b/fastlib/sparse/AminusB.txt new file mode 100755 index 0000000000..2414a239ae --- /dev/null +++ b/fastlib/sparse/AminusB.txt @@ -0,0 +1,73 @@ +1 1 0.31106 +13 1 0.47969 +16 1 -0.11207 +2 2 0.734 +3 3 -0.4987 +4 3 0.6412 +12 3 -0.89951 +16 3 -0.29156 +1 4 0.16853 +3 4 0.39979 +10 4 -0.67427 +11 4 0.87535 +19 4 0.9238 +20 4 0.31775 +10 5 -0.9271 +14 5 0.31602 +19 5 -0.0058848 +1 6 0.66473 +8 6 -0.199 +15 6 0.61663 +17 6 -0.94423 +11 7 0.83526 +12 7 -0.69276 +16 7 -0.097447 +1 8 0.32272 +6 8 0.69778 +18 8 0.51015 +3 9 0.50552 +6 9 0.46189 +10 9 -0.34382 +16 9 -0.39745 +18 9 0.71396 +1 10 -0.47866 +5 10 -0.76755 +13 10 0.56082 +1 11 -0.52652 +14 11 0.6619 +19 11 -0.57442 +20 11 -0.15621 +3 12 0.16931 +6 12 0.082613 +7 12 0.82072 +10 12 -0.59449 +13 12 -0.60971 +16 12 -0.33331 +5 13 0.83685 +6 13 -0.94734 +15 13 0.68514 +3 14 -0.9222 +6 14 -0.81331 +7 14 -0.92383 +8 14 0.19302 +18 14 0.25678 +1 15 -0.79272 +12 15 0.88071 +2 16 -0.19301 +5 16 0.80346 +8 16 0.44535 +10 16 0.30874 +11 16 -0.0033741 +14 16 -0.85604 +3 17 0.52475 +4 18 0.016197 +11 18 0.3331 +12 18 -0.43965 +19 18 0.82212 +4 19 -0.013266 +11 19 -0.98201 +17 19 -0.83856 +18 19 0.60587 +9 20 0.012958 +10 20 -0.61549 +12 20 -0.70102 diff --git a/fastlib/sparse/AplusB.txt b/fastlib/sparse/AplusB.txt new file mode 100755 index 0000000000..6a12d9f139 --- /dev/null +++ b/fastlib/sparse/AplusB.txt @@ -0,0 +1,73 @@ +1 1 0.31106 +13 1 0.47969 +16 1 0.11207 +2 2 0.734 +3 3 1.3205 +4 3 0.6412 +12 3 0.89951 +16 3 0.29156 +1 4 0.16853 +3 4 0.39979 +10 4 0.67427 +11 4 0.87535 +19 4 1.0096 +20 4 0.31775 +10 5 0.9271 +14 5 0.91579 +19 5 0.0058848 +1 6 1.1286 +8 6 0.199 +15 6 0.61663 +17 6 0.94423 +11 7 0.83526 +12 7 0.69276 +16 7 0.097447 +1 8 0.32272 +6 8 0.69778 +18 8 0.51015 +3 9 0.50552 +6 9 0.46189 +10 9 0.34382 +16 9 0.39745 +18 9 0.71396 +1 10 0.47866 +5 10 0.76755 +13 10 0.56082 +1 11 0.52652 +14 11 0.6619 +19 11 0.57442 +20 11 1.3316 +3 12 0.16931 +6 12 0.082613 +7 12 0.82072 +10 12 0.59449 +13 12 0.60971 +16 12 0.33331 +5 13 0.83685 +6 13 0.94734 +15 13 0.68514 +3 14 0.9222 +6 14 0.81331 +7 14 0.92383 +8 14 0.19302 +18 14 0.77364 +1 15 0.79272 +12 15 0.88071 +2 16 0.19301 +5 16 0.80346 +8 16 0.44535 +10 16 0.30874 +11 16 0.0033741 +14 16 0.85604 +3 17 0.52475 +4 18 0.016197 +11 18 0.3331 +12 18 0.43965 +19 18 0.82212 +4 19 0.013266 +11 19 0.98201 +17 19 0.83856 +18 19 0.60587 +9 20 0.012958 +10 20 0.61549 +12 20 0.70102 diff --git a/fastlib/sparse/AtimesB.txt b/fastlib/sparse/AtimesB.txt new file mode 100755 index 0000000000..2afd5e5999 --- /dev/null +++ b/fastlib/sparse/AtimesB.txt @@ -0,0 +1,70 @@ +5 1 0.090046 +8 1 0.049912 +10 1 0.034601 +3 3 0.52605 +4 3 0.58324 +5 3 0.23426 +6 3 0.074311 +7 3 0.73825 +8 3 0.12985 +10 3 0.090017 +13 4 0.37814 +18 4 0.025991 +8 5 0.057884 +13 5 0.51993 +18 5 0.15807 +1 6 0.13636 +3 6 0.49548 +6 6 0.13886 +13 6 0.11125 +18 6 0.10152 +3 7 0.11729 +5 7 0.078295 +6 7 0.05723 +7 7 0.56856 +8 7 0.043398 +10 7 0.030086 +5 9 0.31933 +8 9 0.17701 +10 9 0.12271 +13 9 0.19282 +1 10 0.14889 +13 10 0.22961 +14 10 0.47274 +1 11 0.16378 +9 11 0.0096394 +13 11 0.25257 +18 11 0.34802 +5 12 0.77805 +8 12 0.14844 +10 12 0.10291 +13 12 0.3334 +15 12 0.41774 +1 13 0.84943 +15 13 0.58416 +1 14 0.72925 +3 14 0.37893 +4 14 0.5955 +11 14 0.85772 +15 14 0.50151 +19 14 0.21246 +1 15 0.24658 +13 15 0.38026 +2 16 0.14167 +8 16 0.16523 +14 16 0.0022333 +18 16 0.44104 +20 16 0.0019829 +3 18 0.074436 +6 18 0.036321 +7 18 0.36083 +1 19 0.0022358 +3 19 0.44533 +11 19 0.011613 +14 19 0.64999 +19 19 0.012825 +20 19 0.58134 +3 20 0.11869 +6 20 0.057913 +7 20 0.57534 +13 20 0.34518 diff --git a/fastlib/sparse/B.txt b/fastlib/sparse/B.txt new file mode 100755 index 0000000000..280d8db8f5 --- /dev/null +++ b/fastlib/sparse/B.txt @@ -0,0 +1,39 @@ +16 1 0.11207 +3 3 0.9096 +12 3 0.89951 +16 3 0.29156 +10 4 0.67427 +19 4 0.042898 +10 5 0.9271 +14 5 0.29989 +19 5 0.0058848 +1 6 0.23192 +8 6 0.199 +17 6 0.94423 +12 7 0.69276 +16 7 0.097447 +10 9 0.34382 +16 9 0.39745 +1 10 0.47866 +5 10 0.76755 +1 11 0.52652 +19 11 0.57442 +20 11 0.7439 +10 12 0.59449 +13 12 0.60971 +16 12 0.33331 +6 13 0.94734 +3 14 0.9222 +6 14 0.81331 +7 14 0.92383 +18 14 0.25843 +1 15 0.79272 +2 16 0.19301 +11 16 0.0033741 +14 16 0.85604 +12 18 0.43965 +4 19 0.013266 +11 19 0.98201 +17 19 0.83856 +10 20 0.61549 +12 20 0.70102 diff --git a/fastlib/sparse/build.py b/fastlib/sparse/build.py new file mode 100644 index 0000000000..6c4f11c5c9 --- /dev/null +++ b/fastlib/sparse/build.py @@ -0,0 +1,12 @@ +librule(name="sparse", + headers=["sparse_matrix.h"], +# ] +deplibs=["trilinos:trilinos", "base:base"] + ); + +binrule(name="mtest", + sources=["sparse_matrix_test.cc"], + cflags=" -fexceptions", + deplibs=[":sparse"] + ); + diff --git a/fastlib/sparse/sparse_matrix.h b/fastlib/sparse/sparse_matrix.h new file mode 100644 index 0000000000..dfafce0ac3 --- /dev/null +++ b/fastlib/sparse/sparse_matrix.h @@ -0,0 +1,550 @@ +/* + * ===================================================================================== + * + * Filename: sparse_matrix.h + * + * Description: + * + * Version: 1.0 + * Created: 12/01/2007 04:12:00 PM EST + * Revision: none + * Compiler: gcc + * + * Author: Nikolaos Vasiloglou (NV), nvasil@ieee.org + * Company: Georgia Tech Fastlab-ESP Lab + * + * ===================================================================================== + */ + +#ifndef SPARSE_MATRIX_H_ +#define SPARSE_MATRIX_H_ +#ifndef HAVE_CONFIG_H +#define HAVE_CONFIG_H +#endif +#ifndef USE_TRILINOS +#define USE_TRILINOS +#endif +#include +#include +#include +#include +#include +#include +#include +#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 "trilinos/include/Epetra_CrsMatrix.h" +#include "trilinos/include/Epetra_SerialComm.h" +#include "trilinos/include/Epetra_Map.h" +#include "trilinos/include/Epetra_Vector.h" +#include "trilinos/include/Epetra_MultiVector.h" +#include "trilinos/include/AnasaziBasicEigenproblem.hpp" +#include "trilinos/include/AnasaziEpetraAdapter.hpp" +#include "trilinos/include/AnasaziBlockKrylovSchurSolMgr.hpp" +#include "trilinos/include/AztecOO.h" + +/* class SparseMatrix created by Nick + * This is a sparse matrix wrapper for trilinos Epetra_CrsMatrix + * It is much simpler than Epetra_CrsMatrix. At this time + * it supports eigenvalues (Krylov method) and linear system solution + * I have added matrix addition/subtraction multiplication + * I am also trying to add the submatrices + * Note: There is a restriction on these matrices, the number of rows is + * always greater or equal to the number of columns. The number of rows + * is also called dimension. We pose this restriction because trilinos supports + * square matrices only. In sparse matrices though this is not the problem since + * an mxn matrix where m>n can is equivalent to an mxm matrix where all the + * elements with n MVT; + + SparseMatrix() ; + // Constructor + // num_of_rows: number of rows + // num_of_cols: number of columns + // nnz_per_row: an estimate of the non zero elements per row + // This doesn't need to be accurate. If you need + // more it will automatically resize. Try to be as accurate + // as you can because resizing costs. It is better if your + // estimete if greater than the true non zero elements. So + // it is better to overestimate than underestimate + 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(std::string textfile) { + Init(textfile); + } + ~SparseMatrix() { + Destruct(); + } + void Destruct(); + // Use this initializer like the Constructor + void Init(const index_t num_of_rows, + const index_t num_of_columns, + const index_t nnz_per_row); + // This Initializer is like the previous one with the main difference that + // for every row we give a seperate estimate for the non-zero elements. + void Init(index_t num_of_rows, index_t num_of_columns, index_t *nnz_per_row); + // This Initializer fills the sparse matrix with data. + // row_indices: row indices for non-zero elements + // col_indices: column indices for non-zero elements + // values : values of non-zeros elements + // If the dimension (number of rows)and the expected (nnz elements per row) + // are set to a negative value, the function will automatically detect it + void Init(const std::vector &row_indices, + const std::vector &col_indices, + const Vector &values, + index_t nnz_per_row, + index_t dimension); + // The same as above but we use STL vector for values + void Init(const std::vector &row_indices, + const std::vector &col_indices, + const std::vector &values, + index_t nnz_per_row, + index_t dimension); + + // Initialize from a text file in the following format + // row column value \n + void Init(std::string textfile); + // Initialize the diagonal + void InitDiagonal(const Vector &vec); + // Initialize the diagonal with a constant + void InitDiagonal(const double value); + // It is recomended that you load the matrix row-wise, Before + // you do that call StartLoadingRows() + void StartLoadingRows(); + // All these functions load Rows, with the data in different format + void LoadRow(index_t row, std::vector &columns, Vector &values); + void LoadRow(index_t row, index_t *columns, Vector &values); + void LoadRow(index_t row, index_t num, index_t *columns, double *values); + void LoadRow(index_t row, std::vector &columns, std::vector &values); + // When you are done call this it does some optimization in the storage, no + // further asignment + void EndLoading(); + // It makes the matrix symmetric. It scans the rows of the matrix and for every (i,j) + // element (j,i) equal to (j,i) + void MakeSymmetric(); + // if you know that the matrix is symmetric set the flag + void set_symmetric(bool val) { + issymmetric_ = val; + } + // Not implemented yet + void SetDiagonal(const Vector &vector); + // Copy function, used also by copy constructor + void Copy(const SparseMatrix &other); + // Not implemented yet + void Alias(const Matrix& other); + // Not Implemented yet + void SwapValues(SparseMatrix* other); + // Access values, It will fail if EndLoading() has been called + double get(index_t r, index_t c) const; + // Set Values + void set(index_t r, index_t c, double v); + // For debug purposes you can call it to print the matrix + std::string Print() { + std::ostringstream s1; + matrix_->Print(s1); + return s1.str(); + } + // Get the number of rows + index_t get_num_of_rows() { + return num_of_rows_; + } + // Get the number of columns + index_t get_num_of_columns() { + return num_of_columns_; + } + // Dimension should be equal to the number of rows + index_t get_dimension() { + return dimension_; + } + // The number of non zero elements + index_t get_nnz() { + return matrix_->NumGlobalNonzeros(); + } + // Computes the eignvalues with the Krylov Method + void Eig(index_t num_of_eigvalues, // number of eigenvalues to compute + std::string eigtype, // Choose which eigenvalues to compute + // Choices are: + // LM - target the largest magnitude + // SM - target the smallest magnitude + // LR - target the largest real + // SR - target the smallest real + // LI - target the largest imaginary + // SI - target the smallest imaginary + Matrix *eigvectors, // The eigenvectors computed must not be initialized + std::vector *real_eigvalues, // real part of the eigenvalues + // must be initialized, but should not + // allocate space. The eigenvalues + // returned might actually be less + // than the ones requested + // for example when the matrix has + // rank n< eigenvalues requested + std::vector *imag_eigvalues // imaginary part of the eigenvalues + // must be initialized. If the + // problem is symmetric there is + // no need to initialize. + // The same as real_eigvalues + // hold for the space allocated + // in the non-symmetric case + ); + // Linear System solution, Call Endloading First. + void LinSolve(Vector &b, // must be initialized (space allocated) + Vector *x, // must be initialized (space allocated) + double tolerance, + index_t iterations + ) { + if (matrix_->Filled()==false && matrix_->StorageOptimized()==false){ + FATAL("You should call EndLoading() first\n"); + } + Epetra_Vector tempb(View, *map_, b.ptr()); + Epetra_Vector tempx(View, *map_, x->ptr()); + // create linear problem + Epetra_LinearProblem problem(matrix_.get(), &tempx, &tempb); + // create the AztecOO instance + AztecOO solver(problem); + solver.SetAztecOption( AZ_precond, AZ_Jacobi); + solver.Iterate(iterations, tolerance); + NONFATAL("Solver performed %i iterations, true residual %lg", + solver.NumIters(), solver.TrueResidual()); + } + + // Use this for the general case + void LinSolve(Vector &b, Vector *x) { + LinSolve(b, x, 1E-9, 1000); + } + + // scales the matrix with a scalar; + void Scale(double scalar) { + matrix_->Scale(scalar); + } + // The matrix will be scaled such that A(i,j) = x(j)*A(i,j) + // where i denotes the global row number of A and j denotes the column number + void ColumnScale(const Vector &vec) { + Epetra_Vector temp(View, *map_, (double*)vec.ptr()); + matrix_->RightScale(temp); + } + // The matrix will be scaled such that A(i,j) = x(i)*A(i,j) + // where i denotes the row number of A and j denotes the column number of A. + void RowScale(const Vector &vec) { + Epetra_Vector temp(View, *map_, (double *)vec.ptr()); + matrix_->LeftScale(temp); + } + // computes the L1 norm + double L1Norm() { + return matrix_->NormOne(); + } + // L infinity norm + double LInfNorm() { + return matrix_->NormInf(); + } + // Computes the inverse of the sum of absolute values of the rows + // of the matrix + void InvRowsSums(Vector *result) { + Epetra_Vector temp(View, *map_, result->ptr()); + matrix_->InvRowSums(temp); + } + // Computes the inv of max of absolute values of the rows of the matrix, + void InvRowMaxs(Vector *result) { + Epetra_Vector temp(View, *map_, result->ptr()); + matrix_->InvRowMaxs(temp); + } + // Computes the inverse of the sum of absolute values of the columns of the + // matrix + void InvColSums(Vector *result) { + Epetra_Vector temp(View, *map_, result->ptr()); + matrix_->InvColSums(temp); + } + // Computes the inv of max of absolute values of the columns of the matrix, + void InvColMaxs(Vector *result) { + Epetra_Vector temp(View, *map_, result->ptr()); + matrix_->InvColMaxs(temp); + } + + private: + index_t dimension_; + index_t num_of_rows_; + index_t num_of_columns_; + Epetra_SerialComm comm_; + bool issymmetric_; + Epetra_Map *map_; + Teuchos::RCP matrix_; + index_t *my_global_elements_; + void Load(const std::vector &rows, + const std::vector &columns, + const Vector &values); + void AllRowsLoad(Vector &rows, Vector &columns); + +}; + +class Sparsem { + public: + static inline void Add(const SparseMatrix &a, + const SparseMatrix &b, + SparseMatrix *result) { + DEBUG_ASSERT(a.num_of_rows_==b.num_of_rows_); + DEBUG_ASSERT(a.num_of_columns_==b.num_of_columns_); + DEBUG_ASSERT(a.num_of_rows_==result->num_of_rows_); + DEBUG_ASSERT(a.num_of_columns_==result->num_of_columns_); + result->StartLoadingRows(); + for(index_t r=0; rExtractGlobalRowView(a.my_global_elements_[r], num1, values1, indices1); + b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r], num2, values2, indices2); + std::vector values3; + std::vector indices3; + index_t i=0; + index_t j=0; + while (likely(i=num1)) { + break; + } + } + if ( likely(iLoadRow(r, indices3, values3); + } + } + static inline void Subtract(const SparseMatrix &a, + const SparseMatrix &b, + SparseMatrix *result) { + DEBUG_ASSERT(a.num_of_rows_==b.num_of_rows_); + DEBUG_ASSERT(a.num_of_columns_==b.num_of_columns_); + DEBUG_ASSERT(a.num_of_rows_==result->num_of_rows_); + DEBUG_ASSERT(a.num_of_columns_==result->num_of_columns_); + // If you try assigning the results to an already initialized matrix + // you might get unexpected results. The following assertions + // prevent you partially from that + DEBUG_ASSERT(&a != result); + DEBUG_ASSERT(&b != result); + result->StartLoadingRows(); + for(index_t r=0; rExtractGlobalRowView(a.my_global_elements_[r], num1, values1, indices1); + b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r], num2, values2, indices2); + std::vector values3; + std::vector indices3; + index_t i=0; + index_t j=0; + while (likely(i=num1)) { + break; + } + } + if (likely(iLoadRow(r, indices3, values3); + } + } + /* Multiplication of two matrices A*B in matlab notation + * If B is symmetric then it is much faster, because we can + * multiply rows. Otherwise we have to compute the transpose + * As an advise multiplication of two sparse matrices might + * lead to a dense one, so please be carefull + */ + static inline void Multiply(const SparseMatrix &a, + const SparseMatrix &b, + SparseMatrix *result) { + DEBUG_ASSERT(a.num_of_columns_ == b.num_of_rows_); + DEBUG_ASSERT(a.num_of_rows_ == result->num_of_rows_); + DEBUG_ASSERT(b.num_of_columns_ == result->num_of_columns_); + // If you try assigning the results to an already initialized matrix + // you might get unexpected results. The following assertions + // prevent you partially from that + DEBUG_ASSERT(&a != result); + DEBUG_ASSERT(&b != result); + + if (b.issymmetric_ == true) { + for(index_t r1=0; r1 indices3; + std::vector values3; + index_t num1; + double *values1; + index_t *indices1; + a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r1], + num1, values1, indices1); + for(index_t r2=0; r2ExtractGlobalRowView(b.my_global_elements_[r2], + num2, values2, indices2); + index_t i=0; + index_t j=0; + double dot_product=0; + while (likely(i=num1)) { + break; + } + } + if (likely(iLoadRow(r1, indices3, values3); + indices3.clear(); + values3.clear(); + } + } else { + for(index_t r1=0; r1ExtractGlobalRowView(a.my_global_elements_[r1], + num1, values1, indices1); + double dot_product=0; + std::vector indices3; + std::vector values3; + for(index_t r2=0; r2LoadRow(r1, indices3, values3); + } + indices3.clear(); + values3.clear(); + } + } + } + /* The transpose flag should be set to true if + * we want to use the transpose of mat, otherwise + * set it to false. + * */ + static inline void Multiply(const SparseMatrix &mat, + const Vector &vec, + Vector *result, + bool transpose_flag) { + Epetra_Vector temp_in(View, *(mat.map_), (double *)vec.ptr()); + Epetra_Vector temp_out(View, *(mat.map_), (double *)result->ptr()); + mat.matrix_->Multiply(transpose_flag, temp_in, temp_out); + } + /* Multiply the matrix with a scalar + */ + static inline void Multiply(const SparseMatrix &mat, + const double scalar, + SparseMatrix *result) { + result->Copy(mat); + result->Scale(scalar); + } + /* element wise multiplication of the matrices + * A.*B in matlab notation + */ + static inline void DotMultiply(const SparseMatrix &a, + const SparseMatrix &b, + SparseMatrix *result) { + DEBUG_ASSERT(a.num_of_columns_ == b.num_of_rows_); + DEBUG_ASSERT(a.num_of_rows_ == result->num_of_rows_); + DEBUG_ASSERT(b.num_of_columns_ == result->num_of_columns_); + // If you try assigning the results to an already initialized matrix + // you might get unexpected results. The following assertions + // prevent you partially from that + DEBUG_ASSERT(&a != result); + DEBUG_ASSERT(&b != result); + for(index_t r=0; r indices3; + std::vector values3; + indices3.clear(); + values3.clear(); + index_t num1, num2; + double *values1, *values2; + index_t *indices1, *indices2; + a.matrix_->ExtractGlobalRowView(a.my_global_elements_[r], num1, values1, indices1); + b.matrix_->ExtractGlobalRowView(b.my_global_elements_[r], num2, values2, indices2); + index_t i=0; + index_t j=0; + while (likely(i=num1)) { + break; + } + } + if ( likely(iLoadRow(r, indices3, values3); + } + } +}; +#include "u/nvasil/sparse_matrix/sparse_matrix_impl.h" +#endif diff --git a/fastlib/sparse/sparse_matrix_impl.h b/fastlib/sparse/sparse_matrix_impl.h new file mode 100644 index 0000000000..11ac2e68a2 --- /dev/null +++ b/fastlib/sparse/sparse_matrix_impl.h @@ -0,0 +1,449 @@ +/* + * ===================================================================================== + * + * Filename: sparse_matrix_impl.h + * + * Description: + * + * Version: 1.0 + * Created: 12/02/2007 10:18:02 AM EST + * Revision: none + * Compiler: gcc + * + * Author: Nikolaos Vasiloglou (NV), nvasil@ieee.org + * Company: Georgia Tech Fastlab-ESP Lab + * + * ===================================================================================== + * + */ + +SparseMatrix::SparseMatrix() { + map_ = NULL; + issymmetric_=false; +} + +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) { + issymmetric_ = false; + 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_ = Teuchos::rcp( + new Epetra_CrsMatrix((Epetra_DataAccess)0, *map_, nnz_per_row)); + StartLoadingRows(); +} + +void SparseMatrix::Init(index_t num_of_rows, + index_t num_of_columns, + index_t *nnz_per_row) { + issymmetric_ = false; + 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_ = Teuchos::rcp( + new Epetra_CrsMatrix((Epetra_DataAccess)0 , *map_, nnz_per_row)); + StartLoadingRows(); +} + + +void SparseMatrix::Init(const std::vector &rows, + const std::vector &columns, + const Vector &values, + index_t nnz_per_row, + index_t dimension) { + issymmetric_ = false; + if (nnz_per_row > 0 && dimension > 0) { + Init(dimension, dimension, nnz_per_row); + } else { + num_of_rows_ = 0; + num_of_columns_=0; + std::map frequencies; + for(index_t i=0; i<(index_t)rows.size(); i++) { + frequencies[rows[i]]++; + frequencies[columns[i]]++; + if (rows[i]>num_of_rows_) { + num_of_rows_ = rows[i]; + } + if (columns[i]>num_of_columns_) { + num_of_columns_ = columns[i]; + } + } + num_of_columns_++; + num_of_rows_++; + if (num_of_rows_ < num_of_columns_) { + FATAL("At this point we only support rows (%i) >= columns (%i)\n", + num_of_rows_, num_of_columns_); + } + dimension_ = 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_]; + for(index_t i=0; i &rows, + const std::vector &columns, + const std::vector &values, + index_t nnz_per_row, + index_t dimension) { + issymmetric_ = false; + Vector temp; + temp.Alias((double *)&values[0], values.size()); + Init(rows, columns, temp, nnz_per_row, dimension); + +} + +void SparseMatrix::Init(std::string filename) { + issymmetric_ = false; + FILE *fp = fopen(filename.c_str(), "r"); + if (fp == NULL) { + FATAL("Cannot open %s, error: %s", + filename.c_str(), + strerror(errno)); + } + std::vector rows; + std::vector cols; + std::vector vals; + while (!feof(fp)) { + index_t r, c; + double v; + fscanf(fp,"%i %i %lg\n", &r, &c, &v); + rows.push_back(r); + cols.push_back(c); + 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); +} + +void SparseMatrix::Destruct() { + if (map_!=NULL) { + delete map_; + map_=NULL; + } +} + +void SparseMatrix::Copy(const SparseMatrix &other) { + num_of_rows_ = other.num_of_rows_; + num_of_columns_ = other.num_of_columns_; + dimension_ = other.dimension_; + map_ = new Epetra_Map(num_of_rows_, 0, comm_); + issymmetric_ = other.issymmetric_; + if (other.matrix_->Filled()==false) { + matrix_ = Teuchos::rcp(new Epetra_CrsMatrix(Epetra_DataAccess(0), *map_, 10)); + this->StartLoadingRows(); + for(index_t r=0; rExtractGlobalRowView(global_row, num_of_entries, values, indices); + this->LoadRow(r, num_of_entries, indices, values); + } + } else { + matrix_ = Teuchos::rcp(new Epetra_CrsMatrix(*(other.matrix_.get()))); + } +} +void SparseMatrix::StartLoadingRows() { + my_global_elements_ = map_->MyGlobalElements(); +} + +void SparseMatrix::LoadRow(index_t row, + std::vector &columns, + Vector &values) { + DEBUG_ASSERT(values.length() == (index_t)columns.size()); + matrix_->InsertGlobalValues(my_global_elements_[row], + values.length(), + values.ptr(), + &columns[0]); + +} + +void SparseMatrix::LoadRow(index_t row, + index_t *columns, + Vector &values) { + matrix_->InsertGlobalValues(my_global_elements_[row], + values.length(), + values.ptr(), + &columns[0]); + +} + +void SparseMatrix::LoadRow(index_t row, + std::vector &columns, + std::vector &values) { + matrix_->InsertGlobalValues(my_global_elements_[row], + values.size(), + &values[0], + &columns[0]); + +} + +void SparseMatrix::LoadRow(index_t row, + index_t num, + index_t *columns, + double *values) { + matrix_->InsertGlobalValues(my_global_elements_[row], + num, + values, + columns); + +} + +void SparseMatrix::EndLoading() { + matrix_->FillComplete(); + matrix_->OptimizeStorage(); +} + +void SparseMatrix::Load(const std::vector &rows, + const std::vector &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]; + index_t prev_row= rows[i]; + std::vector indices; + std::vector row_values; + while (true) { + 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_[rows[prev_row]], + row_values.size(), + &row_values[0], + &indices[0]); + prev_row=cur_row; + if (i >= (index_t)rows.size()) { + break; + } + } +} + +double SparseMatrix::get(index_t r, index_t c) const { + DEBUG_BOUNDS(r, num_of_rows_); + DEBUG_BOUNDS(c, num_of_columns_); + index_t global_row = my_global_elements_[r]; + index_t num_of_entries; + double *values; + index_t *indices; + matrix_->ExtractGlobalRowView(global_row, num_of_entries, values, indices); + index_t *pos = std::find(indices, indices+num_of_entries, c); + if (pos==indices+num_of_entries) { + return 0; + } + return values[(ptrdiff_t)(pos-indices)]; +} + +void SparseMatrix::set(index_t r, index_t c, double v) { + DEBUG_BOUNDS(r, num_of_rows_); + DEBUG_BOUNDS(c, num_of_columns_); + if (get(r,c)!=0) { + matrix_->InsertGlobalValues(my_global_elements_[r], 1, &v, &c); + } else { + matrix_->ReplaceGlobalValues(my_global_elements_[r], 1, &v, &c); + } +} +void SparseMatrix::MakeSymmetric() { + index_t num_of_entries; + double *values; + index_t *indices; + for(index_t i=0; iExtractGlobalRowView(global_row, num_of_entries, values, indices); + for(index_t j=0; j *real_eigvalues, + std::vector *imag_eigvalues) { + if (unlikely(!matrix_->Filled())) { + FATAL("You have to call EndLoading before running eigenvalues otherwise " + "it will fail\n"); + } + index_t block_size=2; + Teuchos::RCP ivec = Teuchos::rcp(new + Epetra_MultiVector(*map_, + block_size)); + // Fill it with random numbers + ivec->Random(); + // Setup the eigenproblem, with the matrix A and the initial vectors ivec + Teuchos::RCP >problem = + Teuchos::rcp(new Anasazi::BasicEigenproblem(matrix_, ivec)); + + // The 2-D laplacian is symmetric. Specify this in the eigenproblem. + if (issymmetric_ == true) { + problem->setHermitian(true); + } else { + problem->setHermitian(false); + } + + // Specify the desired number of eigenvalues + problem->setNEV(num_of_eigvalues); + + // Signal that we are done setting up the eigenvalue problem + bool ierr = problem->setProblem(); + + // Check the return from setProblem(). If this is true, there was an + // error. This probably means we did not specify enough information for + // the eigenproblem. + if unlikely(ierr == false) { + FATAL("Trilinos solver error, you probably didn't specify enough information " + "for the eigenvalue problem\n"); + } + + // Specify the verbosity level. Options include: + // Anasazi::Errors + // This option is always set + // Anasazi::Warnings + // Warnings (less severe than errors) + // Anasazi::IterationDetails + // Details at each iteration, such as the current eigenvalues + // Anasazi::OrthoDetails + // Details about orthogonality + // Anasazi::TimingDetails + // A summary of the timing info for the solve() routine + // Anasazi::FinalSummary + // A final summary + // Anasazi::Debug + // Debugging information + int verbosity = Anasazi::Warnings + + Anasazi::Errors + + Anasazi::FinalSummary + + Anasazi::TimingDetails; + + // Choose which eigenvalues to compute + // Choices are: + // LM - target the largest magnitude [default] + // SM - target the smallest magnitude + // LR - target the largest real + // SR - target the smallest real + // LI - target the largest imaginary + // SI - target the smallest imaginary + + // Create the parameter list for the eigensolver + Teuchos::ParameterList my_pl; + my_pl.set( "Verbosity", verbosity); + my_pl.set( "Which", eigtype); + my_pl.set( "Block Size", block_size); + my_pl.set( "Num Blocks", 20); + my_pl.set( "Maximum Restarts", 100); + my_pl.set( "Convergence Tolerance", 1.0e-8); + + // Create the Block Krylov Schur solver + // This takes as inputs the eigenvalue problem and the solver parameters + Anasazi::BlockKrylovSchurSolMgr + my_block_krylov_schur(problem, my_pl); + + // Solve the eigenvalue problem, and save the return code + Anasazi::ReturnType solver_return = my_block_krylov_schur.solve(); + + // Check return code of the solver: Unconverged, Failed, or OK + switch (solver_return) { + // UNCONVERGED + case Anasazi::Unconverged: + NONFATAL("Anasazi::BlockKrylovSchur::solve() did not converge!\n"); + return ; + // CONVERGED + case Anasazi::Converged: + NONFATAL("Anasazi::BlockKrylovSchur::solve() converged!\n"); + } + // Get eigensolution struct + Anasazi::Eigensolution sol = problem->getSolution(); + // Get the number of eigenpairs returned + int num_of_eigvals_returned = sol.numVecs; + if (num_of_eigvals_returned < num_of_eigvalues) { + NONFATAL("The solver returned less eigenvalues (%i) " + "than requested (%i)\n", num_of_eigvalues, + num_of_eigvals_returned); + } + + // Get eigenvectors + eigvectors->Init(dimension_, num_of_eigvals_returned); + Teuchos::RCP evecs = sol.Evecs; + evecs->ExtractCopy(eigvectors->GetColumnPtr(0), dimension_); + + // Get eigenvalues + std::vector > evals = sol.Evals; + real_eigvalues->resize(num_of_eigvals_returned); + for(index_t i=0; iassign(i, evals[i].realpart); + if (issymmetric_ == false) { + imag_eigvalues->resize(num_of_eigvals_returned); + imag_eigvalues->assign(i, evals[i].imagpart); + } + } + + // Test residuals + // Generate a (numev x numev) dense matrix for the eigenvalues... + // This matrix is automatically initialized to zero + Teuchos::SerialDenseMatrix d(num_of_eigvalues, + num_of_eigvalues); + + // Add the eigenvalues on the diagonals (only the real part since problem is Hermitian) + for (int i=0; iApply( *evecs, res); + + // R -= evecs*D + // = A*evecs - evecs*D + MVT::MvTimesMatAddMv( -1.0, *evecs, d, 1.0, res); + + // Compute the 2-norm of each vector in the MultiVector + // and store them to a std::vector + std::vector norm_res(num_of_eigvalues); + MVT::MvNorm(res, &norm_res); +} + diff --git a/fastlib/sparse/sparse_matrix_test.cc b/fastlib/sparse/sparse_matrix_test.cc new file mode 100644 index 0000000000..94c9f21d00 --- /dev/null +++ b/fastlib/sparse/sparse_matrix_test.cc @@ -0,0 +1,270 @@ +/* + * ===================================================================================== + * + * Filename: sparse_matrix_test.cc + * + * Description: + * + * Version: 1.0 + * Created: 12/08/2007 03:50:44 PM EST + * Revision: none + * Compiler: gcc + * + * Author: Nikolaos Vasiloglou (NV), nvasil@ieee.org + * Company: Georgia Tech Fastlab-ESP Lab + * + * ===================================================================================== + */ + +#include +#include +#include +#include +#include +#include +#include +#include "fastlib/fastlib.h" +#include "base/test.h" +#include "sparse/sparse_matrix.h" + +class SparseMatrixTest { + public: + SparseMatrixTest() { + } + ~SparseMatrixTest() { + } + void Init() { + for(index_t i=0; iStartLoadingRows(); + std::vector ind; + std::vector val; + for(index_t i=0; iLoadRow(i, ind, val); + } + + for(index_t i=0; iget(i,j), + mat_[i][j], + std::numeric_limits::epsilon()); + } + } + NONFATAL("TestInit1 sucess!!\n"); + } + void TestInit2() { + std::vector rows; + std::vector cols; + std::vector vals; + std::vector nnz(num_of_rows_); + for(index_t i=0; iInit(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; iget(i,j), mat_[i][j], + std::numeric_limits::epsilon()); + } + } + NONFATAL("TestInit2 success!!\n"); + } + void TestInit3() { + FILE *fp = fopen("temp.txt", "w"); + if (fp==NULL) { + FATAL("Cannot open temp.txt error %s", strerror(errno)); + } + for(index_t i=0; iInit("temp.txt"); + unlink("temp.txt"); + for(index_t i=0; iget(i,j), + mat_[i][j], + std::numeric_limits::epsilon()); + } + } + NONFATAL("TestInit3 success!!"); + } + void TestCopyConstructor() { + smat_ = new SparseMatrix(); + NONFATAL("TestCopyConstructor success!!\n"); + } + void TestMakeSymmetric() { + TestInit1(); + smat_->set(2, 3, 1.44); + smat_->set(3, 2, 0.74); + smat_->set(7, 8, 4.33); + smat_->set(8, 7, 0.22); + smat_->MakeSymmetric(); + for(index_t i=0; iget(i,j), + smat_->get(j, i), + std::numeric_limits::epsilon()); + } + } + NONFATAL("Test MakeSymmetric success!!\n"); + } + void TestEig() { + TestInit1(); + smat_->EndLoading(); + std::vector eigvalues_real; + std::vector eigvalues_imag; + Matrix eigvectors; + smat_->Eig(1, "LM", &eigvectors, &eigvalues_real, &eigvalues_imag); + eigvectors.PrintDebug(); + NONFATAL("Test Eigenvector success!!\n"); + } + void TestLinSolve() { + TestInit1(); + Vector b,x; + b.Init(num_of_cols_); + b.SetZero(); + x.Init(num_of_cols_); + x.SetAll(1); + smat_->MakeSymmetric(); + smat_->EndLoading(); + smat_->LinSolve(b, &x); + x.PrintDebug(); + NONFATAL("Test Linear Solve success!!\n"); + } + void TestBasicOperations() { + SparseMatrix a("A.txt"); + SparseMatrix b("B.txt"); + SparseMatrix a_plus_b("AplusB.txt"); + SparseMatrix a_minus_b("AminusB.txt"); + SparseMatrix a_times_b("AtimesB.txt"); + SparseMatrix a_dot_times_b("AdottimesB.txt"); + SparseMatrix temp; + temp.Init(21,21, 3); + Sparsem::Add(a, b, &temp); + for(index_t i=0; i<20; i++) { + for(index_t j=0; j<20; j++) { + TEST_DOUBLE_APPROX(a_plus_b.get(i,j), temp.get(i,j), 0.01); + } + } + temp.Destruct(); + NONFATAL("Matrix addition sucess!!\n"); + temp.Init(21,21, 3); + Sparsem::Subtract(a, b, &temp); + for(index_t i=0; i<21; i++) { + for(index_t j=0; j<21; j++) { + TEST_DOUBLE_APPROX(a_minus_b.get(i,j), temp.get(i,j), 0.01); + } + } + temp.Destruct(); + NONFATAL("Matrix subtraction success!!\n"); + temp.Init(21,21, 3); + Sparsem::Multiply(a, b, &temp); + for(index_t i=0; i<21; i++) { + for(index_t j=0; j<21; j++) { + TEST_DOUBLE_APPROX(a_times_b.get(i,j), temp.get(i,j), 0.01); + } + } + temp.Destruct(); + NONFATAL("Matrix multiplication success!!\n"); + temp.Init(21, 21, 3); + Sparsem::DotMultiply(a, b, &temp); + for(index_t i=0; i::epsilon()); + } + } + temp.Destruct(); + NONFATAL("Matrix scalar multiplicationn success!!\n"); + + } + void TestAll() { + Init(); + TestInit1(); + Destruct(); + Init(); + TestInit2(); + Destruct(); + Init(); + TestInit3(); + Destruct(); + Init(); + TestCopyConstructor(); + Destruct(); + Init(); + TestMakeSymmetric(); + Destruct(); + Init(); + TestEig(); + Destruct(); + Init(); + TestLinSolve(); + Destruct(); + Init(); + TestBasicOperations(); + } + + private: + SparseMatrix *smat_; + static const index_t num_of_cols_ = 80; + static const index_t num_of_rows_ = 80; + static const index_t num_of_nnz_ = 4; + double mat_[num_of_rows_][num_of_cols_]; + std::vector indices_; + std::vector rows_; +}; + +int main() { + SparseMatrixTest test; + test.TestAll(); +}