From f8a0745cb0426eb3095dbea24288a64eddab04f0 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Wed, 13 Jun 2012 18:26:05 +0200 Subject: Build process... --- bench/spbench/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) (limited to 'bench/spbench/CMakeLists.txt') diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index 079912266..4b3c6f8e3 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -63,3 +63,8 @@ endif(RT_LIBRARY) add_executable(spbenchsolver spbenchsolver.cpp) target_link_libraries (spbenchsolver ${SPARSE_LIBS}) +add_executable(spsolver sp_solver.cpp) +target_link_libraries (spsolver ${SPARSE_LIBS}) + +add_executable(test_sparseLU test_sparseLU.cpp) +target_link_libraries (test_sparseLU ${SPARSE_LIBS}) \ No newline at end of file -- cgit v1.2.3 From 59642da88bf83709e918667680e4ed63af4c31e5 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Thu, 19 Jul 2012 18:03:44 +0200 Subject: Add exception handler to memory allocation --- Eigen/src/OrderingMethods/Eigen_Colamd.h | 8 +- Eigen/src/OrderingMethods/Ordering.h | 6 -- Eigen/src/SparseLU/SparseLU.h | 3 - Eigen/src/SparseLU/SparseLU_Coletree.h | 2 - Eigen/src/SparseLU/SparseLU_Matrix.h | 1 - Eigen/src/SparseLU/SparseLU_Memory.h | 147 +++++++++++++++++------------- Eigen/src/SparseLU/SparseLU_column_bmod.h | 1 - Eigen/src/SparseLU/SparseLU_panel_dfs.h | 1 - Eigen/src/SparseLU/SparseLU_snode_bmod.h | 1 - bench/spbench/CMakeLists.txt | 2 +- bench/spbench/test_sparseLU.cpp | 21 +++-- 11 files changed, 99 insertions(+), 94 deletions(-) (limited to 'bench/spbench/CMakeLists.txt') diff --git a/Eigen/src/OrderingMethods/Eigen_Colamd.h b/Eigen/src/OrderingMethods/Eigen_Colamd.h index 39701d0af..0af137d54 100644 --- a/Eigen/src/OrderingMethods/Eigen_Colamd.h +++ b/Eigen/src/OrderingMethods/Eigen_Colamd.h @@ -155,7 +155,6 @@ #endif /* MATLAB_MEX_FILE */ // == Row and Column structures == - typedef struct EIGEN_Colamd_Col_struct { int start ; /* index for A of first row in this column, or EIGEN_DEAD */ @@ -248,11 +247,9 @@ void eigen_colamd_report (int stats [EIGEN_COLAMD_STATS]); int eigen_init_rows_cols (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col col [], int A [], int p [], int stats[EIGEN_COLAMD_STATS] ); -void eigen_init_scoring (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], - double knobs[EIGEN_COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg); +void eigen_init_scoring (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], double knobs[EIGEN_COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg); -int eigen_find_ordering (int n_row, int n_col, int Alen, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], - int n_col2, int max_deg, int pfree); +int eigen_find_ordering (int n_row, int n_col, int Alen, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], int n_col2, int max_deg, int pfree); void eigen_order_children (int n_col, EIGEN_Colamd_Col Col [], int p []); @@ -2514,5 +2511,4 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E } #endif /* NDEBUG */ - #endif diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index cbd2e5d34..47cd6f169 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h @@ -124,9 +124,6 @@ class COLAMDOrdering typedef PermutationMatrix PermutationType; typedef Matrix IndexVector; /** Compute the permutation vector form a sparse matrix */ - - - template void operator() (const MatrixType& mat, PermutationType& perm) { @@ -152,9 +149,6 @@ class COLAMDOrdering for (int i = 0; i < n; i++) perm.indices()(p(i)) = i; } - - private: - }; diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index bb1decc4c..25fad0f29 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -339,9 +339,6 @@ void SparseLU::analyzePattern(const MatrixType& mat) //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. - // Compute the fill-reducing ordering - // TODO Currently, the only available ordering method is AMD. - OrderingType ord; ord(mat,m_perm_c); //FIXME Check the right semantic behind m_perm_c diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h index 1329d383f..142f4995e 100644 --- a/Eigen/src/SparseLU/SparseLU_Coletree.h +++ b/Eigen/src/SparseLU/SparseLU_Coletree.h @@ -94,7 +94,6 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent) int rset, cset, rroot; for (col = 0; col < nc; col++) { -// cset = pp(col) = col; // Initially, each element is in its own set //FIXME pp(col) = col; cset = col; root(cset) = col; @@ -108,7 +107,6 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent) if (rroot != col) { parent(rroot) = col; -// cset = pp(cset) = rset; // Get the union of cset and rset //FIXME pp(cset) = rset; cset = rset; root(cset) = col; diff --git a/Eigen/src/SparseLU/SparseLU_Matrix.h b/Eigen/src/SparseLU/SparseLU_Matrix.h index 90a0f2740..9f2dcaa56 100644 --- a/Eigen/src/SparseLU/SparseLU_Matrix.h +++ b/Eigen/src/SparseLU/SparseLU_Matrix.h @@ -192,7 +192,6 @@ class SuperNodalMatrix protected: Index m_row; // Number of rows Index m_col; // Number of columns -// Index m_nnz; // Number of nonzero values Index m_nsuper; // Number of supernodes Scalar* m_nzval; //array of nonzero values packed by column Index* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h index a17079199..7a2ab93df 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h @@ -78,41 +78,82 @@ int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex VectorType old_vec; // Temporary vector to hold the previous values if (nbElts > 0 ) - old_vec = vec.segment(0,nbElts); // old_vec should be of size nbElts... to be checked + old_vec = vec.segment(0,nbElts); - //expand the current vector //FIXME Should be in a try ... catch region - vec.resize(new_len); - /* - * Test if the memory has been well allocated - * otherwise reduce the size and try to reallocate - * copy data from previous vector (if exists) to the newly allocated vector - */ - if ( num_expansions != 0 ) // The memory has been expanded before + //Allocate or expand the current vector + try + { + vec.resize(new_len); + } + catch(std::bad_alloc& ) { - int tries = 0; - if (keep_prev) + if ( !num_expansions ) { - if (!vec.size()) return new_len ; + // First time to allocate from LUMemInit() + throw; // Pass the exception to LUMemInit() which has a try... catch block + } + if (keep_prev) + { + // In this case, the memory length should not not be reduced + return new_len; } else { - while (!vec.size()) + // Reduce the size and increase again + int tries = 0; // Number of attempts + do { - // Reduce the size and allocate again - if ( ++tries > 10) return new_len; - alpha = LU_Reduce(alpha); + alpha = LU_Reduce(alpha); new_len = alpha * length ; - vec.resize(new_len); //FIXME Should be in a try catch section - } - } // end allocation - - //Copy the previous values to the newly allocated space - if (nbElts > 0) - vec.segment(0, nbElts) = old_vec; - } // end expansion + try + { + vec.resize(new_len); + } + catch(std::bad_alloc& ) + { + tries += 1; + if ( tries > 10) return new_len; + } + } while (!vec.size()); + } + } + //Copy the previous values to the newly allocated space + if (nbElts > 0) + vec.segment(0, nbElts) = old_vec; + + length = new_len; if(num_expansions) ++num_expansions; return 0; + + /* + * Test if the memory has been well allocated + * otherwise reduce the size and try to reallocate + * copy data from previous vector (if exists) to the newly allocated vector + */ +// if ( num_expansions != 0 ) // The memory has been expanded before +// { +// int tries = 0; +// if (keep_prev) +// { +// if (!vec.size()) return new_len ; +// } +// else +// { +// while (!vec.size()) +// { +// // Reduce the size and allocate again +// if ( ++tries > 10) return new_len; +// alpha = LU_Reduce(alpha); +// new_len = alpha * length ; +// vec.resize(new_len); //FIXME Should be in a try catch section +// } +// } // end allocation +// +// //Copy the previous values to the newly allocated space +// if (nbElts > 0) +// vec.segment(0, nbElts) = old_vec; +// } // end expansion } /** @@ -122,8 +163,8 @@ int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex * \param annz number of initial nonzeros in the matrix * \param lwork if lwork=-1, this routine returns an estimated size of the required memory * \param glu persistent data to facilitate multiple factors : will be deleted later ?? - * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated when memory allocation failed - * NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the row permutation + * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated memory when allocation failed, and 0 on success + * NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation */ template int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, LU_GlobalLU_t& glu) @@ -159,27 +200,26 @@ int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size, glu.xusub.resize(n+1); // Reserve memory for L/U factors - expand(glu.lusup, nzlumax, 0, 0, num_expansions); - expand(glu.ucol,nzumax, 0, 0, num_expansions); - expand(glu.lsub,nzlmax, 0, 0, num_expansions); - expand(glu.usub,nzumax, 0, 1, num_expansions); - - // Check if the memory is correctly allocated, - // FIXME Should be a try... catch section here - while ( !glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size()) + do { - //Reduce the estimated size and retry - nzlumax /= 2; - nzumax /= 2; - nzlmax /= 2; - - if (nzlumax < annz ) return nzlumax; + try + { + expand(glu.lusup, nzlumax, 0, 0, num_expansions); + expand(glu.ucol,nzumax, 0, 0, num_expansions); + expand(glu.lsub,nzlmax, 0, 0, num_expansions); + expand(glu.usub,nzumax, 0, 1, num_expansions); + } + catch(std::bad_alloc& ) + { + //Reduce the estimated size and retry + nzlumax /= 2; + nzumax /= 2; + nzlmax /= 2; + if (nzlumax < annz ) return nzlumax; + } - expand(glu.lusup, nzlumax, 0, 0, num_expansions); - expand(glu.ucol, nzumax, 0, 0, num_expansions); - expand(glu.lsub, nzlmax, 0, 0, num_expansions); - expand(glu.usub, nzumax, 0, 1, num_expansions); - } + } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size()); + ++num_expansions; @@ -207,23 +247,6 @@ int LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int if (failed_size) return failed_size; - // The following code is not really needed since maxlen is passed by reference - // and correspond to the appropriate field in glu -// switch ( mem_type ) { -// case LUSUP: -// glu.nzlumax = maxlen; -// break; -// case UCOL: -// glu.nzumax = maxlen; -// break; -// case LSUB: -// glu.nzlmax = maxlen; -// break; -// case USUB: -// glu.nzumax = maxlen; -// break; -// } - return 0 ; } diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h index 3042eb5f8..00787721b 100644 --- a/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h @@ -133,7 +133,6 @@ int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, Sca // Dense triangular solve -- start effective triangle luptr += nsupr * no_zeros + no_zeros; // Form Eigen matrix and vector -// std::cout<< "jcol " << jcol << " rows " << segsize << std::endl; Map, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); VectorBlock u(tempv, 0, segsize); diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h index 908ee67ac..70ea0f51f 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h +++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h @@ -92,7 +92,6 @@ void LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, Index int xdfs, maxdfs, kpar; // Initialize pointers -// IndexVector& marker1 = marker.block(m, m); VectorBlock marker1(marker, m, m); nseg = 0; IndexVector& xsup = glu.xsup; diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h index 44438d037..fc8042f52 100644 --- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h @@ -80,7 +80,6 @@ int LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, LU_Glob // Update the trailing part of the column jcol U(jcol:jcol+nrow, jcol) using L(jcol:jcol+nrow, fsupc:jcol) and U(fsupc:jcol) new (&A) Map,0,OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); -// Map > l(&(lusup.data()[ufirst+nsupc], nrow); VectorBlock l(lusup, ufirst+nsupc, nrow); l = l - A * u; } diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index 4b3c6f8e3..a093cc5d9 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -67,4 +67,4 @@ add_executable(spsolver sp_solver.cpp) target_link_libraries (spsolver ${SPARSE_LIBS}) add_executable(test_sparseLU test_sparseLU.cpp) -target_link_libraries (test_sparseLU ${SPARSE_LIBS}) \ No newline at end of file +target_link_libraries (test_sparseLU ${SPARSE_LIBS}) diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp index 6fbf03454..31273add5 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp @@ -13,13 +13,14 @@ using namespace Eigen; int main(int argc, char **args) { - SparseMatrix A; - typedef SparseMatrix::Index Index; - typedef Matrix DenseMatrix; - typedef Matrix DenseRhs; - VectorXd b, x, tmp; -// SparseLU, AMDOrdering > solver; - SparseLU, COLAMDOrdering > solver; + typedef complex scalar; + SparseMatrix A; + typedef SparseMatrix::Index Index; + typedef Matrix DenseMatrix; + typedef Matrix DenseRhs; + Matrix b, x, tmp; +// SparseLU, AMDOrdering > solver; + SparseLU, COLAMDOrdering > solver; ifstream matrix_file; string line; int n; @@ -36,7 +37,7 @@ int main(int argc, char **args) if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; } if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;} if (sym != 0) { // symmetric matrices, only the lower part is stored - SparseMatrix temp; + SparseMatrix temp; temp = A; A = temp.selfadjointView(); } @@ -72,8 +73,8 @@ int main(int argc, char **args) timer.stop(); cout << "solve time " << timer.value() << std::endl; /* Check the accuracy */ - VectorXd tmp2 = b - A*x; - double tempNorm = tmp2.norm()/b.norm(); + Matrix tmp2 = b - A*x; + scalar tempNorm = tmp2.norm()/b.norm(); cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; cout << "Number of nonzeros in the factor : " << solver.nnzL() + solver.nnzU() << std::endl; -- cgit v1.2.3 From 4d3b7e2a1351d60b9ee26d0fe3442cd5b3a1f8a9 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Mon, 6 Aug 2012 14:55:02 +0200 Subject: Add support for Metis fill-reducing ordering ; it is generally more efficient than COLAMD ordering --- Eigen/MetisSupport | 26 +++++++ Eigen/src/MetisSupport/CMakeLists.txt | 6 ++ Eigen/src/MetisSupport/MetisSupport.h | 138 ++++++++++++++++++++++++++++++++++ bench/spbench/CMakeLists.txt | 6 ++ bench/spbench/test_sparseLU.cpp | 8 ++ cmake/FindMetis.cmake | 3 +- 6 files changed, 186 insertions(+), 1 deletion(-) create mode 100644 Eigen/MetisSupport create mode 100644 Eigen/src/MetisSupport/CMakeLists.txt create mode 100644 Eigen/src/MetisSupport/MetisSupport.h (limited to 'bench/spbench/CMakeLists.txt') diff --git a/Eigen/MetisSupport b/Eigen/MetisSupport new file mode 100644 index 000000000..a44086ad9 --- /dev/null +++ b/Eigen/MetisSupport @@ -0,0 +1,26 @@ +#ifndef EIGEN_METISSUPPORT_MODULE_H +#define EIGEN_METISSUPPORT_MODULE_H + +#include "SparseCore" + +#include "src/Core/util/DisableStupidWarnings.h" + +extern "C" { +#include +} + + +/** \ingroup Support_modules + * \defgroup MetisSupport_Module MetisSupport module + * + * \code + * #include + * \endcode + */ + + +#include "src/MetisSupport/MetisSupport.h" + +#include "src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_METISSUPPORT_MODULE_H diff --git a/Eigen/src/MetisSupport/CMakeLists.txt b/Eigen/src/MetisSupport/CMakeLists.txt new file mode 100644 index 000000000..2bad31416 --- /dev/null +++ b/Eigen/src/MetisSupport/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_MetisSupport_SRCS "*.h") + +INSTALL(FILES + ${Eigen_MetisSupport_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/MetisSupport COMPONENT Devel + ) diff --git a/Eigen/src/MetisSupport/MetisSupport.h b/Eigen/src/MetisSupport/MetisSupport.h new file mode 100644 index 000000000..a762d96f6 --- /dev/null +++ b/Eigen/src/MetisSupport/MetisSupport.h @@ -0,0 +1,138 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#ifndef METIS_SUPPORT_H +#define METIS_SUPPORT_H + +namespace Eigen { +/** + * Get the fill-reducing ordering from the METIS package + * + * If A is the original matrix and Ap is the permuted matrix, + * the fill-reducing permutation is defined as follows : + * Row (column) i of A is the matperm(i) row (column) of Ap. + * WARNING: As computed by METIS, this corresponds to the vector iperm (instead of perm) + */ +template +class MetisOrdering +{ +public: + typedef PermutationMatrix PermutationType; + typedef Matrix IndexVector; + + template + void get_symmetrized_graph(const MatrixType& A) + { + Index m = A.cols(); + + // Get the transpose of the input matrix + MatrixType At = A.transpose(); + // Get the number of nonzeros elements in each row/col of At+A + Index TotNz = 0; + IndexVector visited(m); + visited.setConstant(-1); + for (int j = 0; j < m; j++) + { + // Compute the union structure of of A(j,:) and At(j,:) + visited(j) = j; // Do not include the diagonal element + // Get the nonzeros in row/column j of A + for (typename MatrixType::InnerIterator it(A, j); it; ++it) + { + Index idx = it.index(); // Get the row index (for column major) or column index (for row major) + if (visited(idx) != j ) + { + visited(idx) = j; + ++TotNz; + } + } + //Get the nonzeros in row/column j of At + for (typename MatrixType::InnerIterator it(At, j); it; ++it) + { + Index idx = it.index(); + if(visited(idx) != j) + { + visited(idx) = j; + ++TotNz; + } + } + } + // Reserve place for A + At + m_indexPtr.resize(m+1); + m_innerIndices.resize(TotNz); + + // Now compute the real adjacency list of each column/row + visited.setConstant(-1); + Index CurNz = 0; + for (int j = 0; j < m; j++) + { + m_indexPtr(j) = CurNz; + + visited(j) = j; // Do not include the diagonal element + // Add the pattern of row/column j of A to A+At + for (typename MatrixType::InnerIterator it(A,j); it; ++it) + { + Index idx = it.index(); // Get the row index (for column major) or column index (for row major) + if (visited(idx) != j ) + { + visited(idx) = j; + m_innerIndices(CurNz) = idx; + CurNz++; + } + } + //Add the pattern of row/column j of At to A+At + for (typename MatrixType::InnerIterator it(At, j); it; ++it) + { + Index idx = it.index(); + if(visited(idx) != j) + { + visited(idx) = j; + m_innerIndices(CurNz) = idx; + ++CurNz; + } + } + } + m_indexPtr(m) = CurNz; + } + + template + void operator() (const MatrixType& A, PermutationType& matperm) + { + Index m = A.cols(); + IndexVector perm(m),iperm(m); + // First, symmetrize the matrix graph. + get_symmetrized_graph(A); + int output_error; + + // Call the fill-reducing routine from METIS + output_error = METIS_NodeND(&m, m_indexPtr.data(), m_innerIndices.data(), NULL, NULL, perm.data(), iperm.data()); + + if(output_error != METIS_OK) + { + //FIXME The ordering interface should define a class of possible errors + std::cerr << "ERROR WHILE CALLING THE METIS PACKAGE \n"; + return; + } + + // Get the fill-reducing permutation + //NOTE: If Ap is the permuted matrix then perm and iperm vectors are defined as follows + // Row (column) i of Ap is the perm(i) row(column) of A, and row (column) i of A is the iperm(i) row(column) of Ap + + // To be consistent with the use of the permutation in SparseLU module, we thus keep the iperm vector + matperm.resize(m); + for (int j = 0; j < m; j++) + matperm.indices()(j) = iperm(j); + + } + + protected: + IndexVector m_indexPtr; // Pointer to the adjacenccy list of each row/column + IndexVector m_innerIndices; // Adjacency list +}; + +}// end namespace eigen +#endif \ No newline at end of file diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index a093cc5d9..2eb0befa9 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -66,5 +66,11 @@ target_link_libraries (spbenchsolver ${SPARSE_LIBS}) add_executable(spsolver sp_solver.cpp) target_link_libraries (spsolver ${SPARSE_LIBS}) +if(METIS_FOUND) + include_directories(${METIS_INCLUDES}) + set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) + add_definitions("-DEIGEN_METIS_SUPPORT") +endif(METIS_FOUND) + add_executable(test_sparseLU test_sparseLU.cpp) target_link_libraries (test_sparseLU ${SPARSE_LIBS}) diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp index 59f8252d0..8c78b0c9b 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp @@ -7,6 +7,9 @@ #include #include #include +#ifdef EIGEN_METIS_SUPPORT +#include +#endif using namespace std; using namespace Eigen; @@ -21,7 +24,12 @@ int main(int argc, char **args) typedef Matrix DenseRhs; Matrix b, x, tmp; // SparseLU, AMDOrdering > solver; +#ifdef EIGEN_METIS_SUPPORT + SparseLU, MetisOrdering > solver; +#else SparseLU, COLAMDOrdering > solver; +#endif + ifstream matrix_file; string line; int n; diff --git a/cmake/FindMetis.cmake b/cmake/FindMetis.cmake index e4d6ef258..627c3e9ae 100644 --- a/cmake/FindMetis.cmake +++ b/cmake/FindMetis.cmake @@ -12,10 +12,11 @@ find_path(METIS_INCLUDES ${INCLUDE_INSTALL_DIR} PATH_SUFFIXES metis + include ) -find_library(METIS_LIBRARIES metis PATHS $ENV{METISDIR} ${LIB_INSTALL_DIR}) +find_library(METIS_LIBRARIES metis PATHS $ENV{METISDIR} ${LIB_INSTALL_DIR} PATH_SUFFIXES lib) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(METIS DEFAULT_MSG -- cgit v1.2.3 From 063705b5be5a41e324773887d3d5ae065321a719 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Fri, 7 Sep 2012 13:14:57 +0200 Subject: Add tutorial for sparse solvers --- Eigen/src/SparseCore/SparseMatrix.h | 2 +- Eigen/src/SparseLU/SparseLU.h | 4 +- Eigen/src/SuperLUSupport/SuperLUSupport.h | 4 +- bench/spbench/CMakeLists.txt | 1 + doc/I17_SparseLinearSystems.dox | 110 ++++++++++++++++++++++++++++++ 5 files changed, 115 insertions(+), 6 deletions(-) create mode 100644 doc/I17_SparseLinearSystems.dox (limited to 'bench/spbench/CMakeLists.txt') diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 52a9dab70..87f3fb873 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -478,7 +478,7 @@ class SparseMatrix } /** Turns the matrix into the uncompressed mode */ - void Uncompress() + void uncompress() { if(m_innerNonZeros != 0) return; diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 6a6579493..e2076138a 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -271,7 +271,7 @@ void SparseLU::analyzePattern(const MatrixType& mat) //First copy the whole input matrix. m_mat = mat; - m_mat.Uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. + m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. //Then, permute only the column pointers for (int i = 0; i < mat.cols(); i++) { @@ -356,7 +356,7 @@ void SparseLU::factorize(const MatrixType& matrix) // Apply the column permutation computed in analyzepattern() // m_mat = matrix * m_perm_c.inverse(); m_mat = matrix; - m_mat.Uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. + m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. //Then, permute only the column pointers for (int i = 0; i < matrix.cols(); i++) { diff --git a/Eigen/src/SuperLUSupport/SuperLUSupport.h b/Eigen/src/SuperLUSupport/SuperLUSupport.h index e3fae4a36..faefd8169 100644 --- a/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -627,9 +627,7 @@ void SuperLU::factorize(const MatrixType& a) this->initFactorization(a); - //DEBUG -// m_sluOptions.ColPerm = COLAMD; - m_sluOptions.Equil = NO; + m_sluOptions.ColPerm = COLAMD; int info = 0; RealScalar recip_pivot_growth, rcond; RealScalar ferr, berr; diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index 2eb0befa9..5451843b9 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -74,3 +74,4 @@ endif(METIS_FOUND) add_executable(test_sparseLU test_sparseLU.cpp) target_link_libraries (test_sparseLU ${SPARSE_LIBS}) + diff --git a/doc/I17_SparseLinearSystems.dox b/doc/I17_SparseLinearSystems.dox new file mode 100644 index 000000000..740bee18e --- /dev/null +++ b/doc/I17_SparseLinearSystems.dox @@ -0,0 +1,110 @@ +namespace Eigen { +/** \page TopicSparseSystems Solving Sparse Linear Systems +In Eigen, there are several methods available to solve linear systems when the coefficient matrix is sparse. Because of the special representation of this class of matrices, special care should be taken in order to get a good performance. See \ref TutorialSparse for a detailed introduction about sparse matrices in Eigen. In this page, we briefly present the main steps that are common to all the linear solvers in Eigen together with the main concepts behind them. Depending on the properties of the matrix, the desired accuracy, the end-user is able to tune these steps in order to improve the performance of its code. However, an impatient user does not need to know deeply what's hiding behind these steps: the last section presents a benchmark routine that can be easily used to get an insight on the performance of all the available solvers. + +\b Table \b of \b contents \n + - \ref TheSparseCompute + - \ref TheSparseSolve + - \ref BenchmarkRoutine + + As summarized in \ref TutorialSparseDirectSolvers, there are many built-in solvers in Eigen as well as interface to external solvers libraries. All these solvers follow the same calling sequence. The basic steps are as follows : +\code +#include +// ... +SparseMatrix A; +// fill A +VectorXd b, x; +// fill b +// solve Ax = b +SolverClassName > solver; +solver.compute(A); +if(solver.info()!=Succeeded) { + // decomposition failed + return; +} +x = solver.solve(b); +if(solver.info()!=Succeeded) { + // solving failed + return; +} +\endcode + +\section TheSparseCompute The Compute Step +In the compute() function, the matrix is generally factorized: LLT for self-adjoint matrices, LDLT for general hermitian matrices and LU for non hermitian matrices. These are the results of using direct solvers. For this class of solvers precisely, the compute step is further subdivided into analyzePattern() and factorize(). + +The goal of analyzePattern() is to reorder the nonzero elements of the matrix, such that the factorization step creates less fill-in. This step exploits only the structure of the matrix. Hence, the results of this step can be used for other linear systems where the matrix has the same structure. Note however that sometimes, some external solvers (like SuperLU) require that the values of the matrix are set in this step, for instance to equilibrate the rows and columns of the matrix. In this situation, the results of this step can note be used with other matrices. + +Eigen provides a limited set of methods to reorder the matrix in this step, either built-in (COLAMD, AMD) or external (METIS). These methods are set in template parameter list of the solver : +\code +DirectSolverClassName, OrderingMethod > solver; +\endcode + +See \link Ordering_Modules the Ordering module \endlink for the list of available methods and the associated options. + +In factorize(), the factors of the coefficient matrix are computed. This step should be called each time the values of the matrix change. However, the structural pattern of the matrix should not change between multiple calls. + +For iterative solvers, the compute step is used to eventually setup a preconditioner. Remember that, basically, the goal of the preconditioner is to speedup the convergence of an iterative method by solving a modified linear system where the coefficient matrix has more clustered eigenvalues. For real problems, an iterative solver should always be used with a preconditioner. In Eigen, a preconditioner is selected by simply adding it as a template parameter to the iterative solver object. +\code +IterativeSolverClassName, PreconditionerName > solver; +\endcode + +FIXME How to get a reference to the preconditioner, in order to set the parameters + +For instance, with the ILUT preconditioner, the incomplete factors L and U are computed in this step. +See \link Sparse_modules the Sparse module \endlink for the list of available preconditioners in Eigen. +\section TheSparseSolve The Solve step +The solve() function computes the solution of the linear systems with one or many right hand sides. +\code +X = solver.solve(B); +\endcode +Here, B can be a vector or a matrix where the columns form the different right hand sides. The solve() function can be called several times as well, for instance When all the right hand sides are not available at once. +\code +x1 = solver.solve(b1); +// Get the second right hand side b2 +x2 = solver.solve(b2); +// ... +\endcode +For direct methods, the solution are computed at the machine precision. Sometimes, the solution need not be too accurate. In this case, the iterative methods are more suitable and the desired accuracy can be set before the solve step using setTolerance(). For all the available functions, please, refer to the documentation of the \link IterativeLinearSolvers_module Iterative solvers module \endlink. + +\section BenchmarkRoutine +Most of the time, all you need is to know how much time it will take to qolve your system, and hopefully, what is the most suitable solver. In Eigen, we provide a benchmark routine that can be used for this purpose. It is very easy to use. First, it should be activated at the configuration step with the flag TEST_REAL_CASES. Then, in bench/spbench, you can compile the routine by typing \b make \e spbenchsolver. You can then run it with --help option to get the list of all available options. Basically, the matrices to test should be in \link http://math.nist.gov/MatrixMarket/formats.html MatrixMarket Coordinate format \endlink, and the routine returns the statistics from all available solvers in Eigen. + +The following table gives an example of XHTML statistics from several Eigen built-in and external solvers. + + +
Matrix N NNZ UMFPACK SUPERLU PASTIX LU BiCGSTAB BiCGSTAB+ILUT GMRES+ILUT LDLT CHOLMOD LDLT PASTIX LDLT LLT CHOLMOD SP LLT CHOLMOD LLT PASTIX LLT CG
vector_graphics 12855 72069 Compute Time 0.02545490.02156770.07018270.0001533880.01401070.01537090.01016010.009305020.0649689 +
Solve Time 0.003378350.0009518260.004843730.03748860.00464450.008477540.0005418130.0002936960.00485376 +
Total Time 0.02883330.02251950.07502650.0376420.01865520.02384840.01070190.009598710.0698227 +
Error(Iter) 1.299e-16 2.04207e-16 4.83393e-15 3.94856e-11 (80) 1.03861e-12 (3) 5.81088e-14 (6) 1.97578e-16 1.83927e-16 4.24115e-15 +
poisson_SPD 19788 308232 Compute Time 0.4250261.823780.6173670.0004789211.340011.334710.7964190.8575730.4730070.8148260.1847190.8615550.4705590.000458188 +
Solve Time 0.02800530.01944020.02687470.2494370.05484440.09269910.008502040.00531710.02589320.008746030.005781550.005303610.02489420.239093 +
Total Time 0.4530311.843220.6442410.2499161.394861.427410.8049210.8628910.49890.8235720.1905010.8668590.4954530.239551 +
Error(Iter) 4.67146e-16 1.068e-15 1.3397e-15 6.29233e-11 (201) 3.68527e-11 (6) 3.3168e-15 (16) 1.86376e-15 1.31518e-16 1.42593e-15 3.45361e-15 3.14575e-16 2.21723e-15 7.21058e-16 9.06435e-12 (261) +
sherman2 1080 23094 Compute Time 0.006317540.0150520.0247514 -0.02144250.0217988 +
Solve Time 0.0004784240.0003379980.0010291 -0.002431520.00246152 +
Total Time 0.006795970.015390.0257805 -0.0238740.0242603 +
Error(Iter) 1.83099e-15 8.19351e-15 2.625e-14 1.3678e+69 (1080) 4.1911e-12 (7) 5.0299e-13 (12) +
bcsstk01_SPD 48 400 Compute Time 0.0001690790.000107890.0005725381.425e-069.1612e-058.3985e-055.6489e-057.0913e-050.0004682515.7389e-058.0212e-055.8394e-050.0004630171.333e-06 +
Solve Time 1.2288e-051.1124e-050.0002863878.5896e-051.6381e-051.6984e-053.095e-064.115e-060.0003254383.504e-067.369e-063.454e-060.0002940956.0516e-05 +
Total Time 0.0001813670.0001190140.0008589258.7321e-050.0001079930.0001009695.9584e-057.5028e-050.0007936896.0893e-058.7581e-056.1848e-050.0007571126.1849e-05 +
Error(Iter) 1.03474e-16 2.23046e-16 2.01273e-16 4.87455e-07 (48) 1.03553e-16 (2) 3.55965e-16 (2) 2.48189e-16 1.88808e-16 1.97976e-16 2.37248e-16 1.82701e-16 2.71474e-16 2.11322e-16 3.547e-09 (48) +
sherman1 1000 3750 Compute Time 0.002288050.002092310.005282689.846e-060.001635220.001621550.0007892590.0008044950.00438269 +
Solve Time 0.0002137889.7983e-050.0009388310.006298350.0003617640.000787944.3989e-052.5331e-050.000917166 +
Total Time 0.002501840.002190290.006221510.00630820.001996980.002409490.0008332480.0008298260.00529986 +
Error(Iter) 1.16839e-16 2.25968e-16 2.59116e-16 3.76779e-11 (248) 4.13343e-11 (4) 2.22347e-14 (10) 2.05861e-16 1.83555e-16 1.02917e-15 +
young1c 841 4089 Compute Time 0.002358430.002172280.005680751.2735e-050.002648660.00258236 +
Solve Time 0.0003295990.0001686340.000801180.05347380.001871930.00450211 +
Total Time 0.002688030.002340910.006481930.05348650.004520590.00708447 +
Error(Iter) 1.27029e-16 2.81321e-16 5.0492e-15 8.0507e-11 (706) 3.00447e-12 (8) 1.46532e-12 (16) +
mhd1280b 1280 22778 Compute Time 0.002348980.002070790.005709182.5976e-050.003025630.002980360.001445250.0009199220.00426444 +
Solve Time 0.001033920.0002119110.001050.01104320.0006282870.003920890.0001383036.2446e-050.00097564 +
Total Time 0.00338290.00228270.006759180.01106920.003653920.006901240.001583550.0009823680.00524008 +
Error(Iter) 1.32953e-16 3.08646e-16 6.734e-16 8.83132e-11 (40) 1.51153e-16 (1) 6.08556e-16 (8) 1.89264e-16 1.97477e-16 6.68126e-09 +
crashbasis 160000 1750416 Compute Time 3.20195.789215.75730.003835153.10063.09921 +
Solve Time 0.2619150.1062250.4021411.490890.248880.443673 +
Total Time 3.463815.8954216.15941.494733.349483.54288 +
Error(Iter) 1.76348e-16 4.58395e-16 1.67982e-14 8.64144e-11 (61) 8.5996e-12 (2) 6.04042e-14 (5) + +
+*/ +} \ No newline at end of file -- cgit v1.2.3 From 2c99d8413316c97e771a37c7ff04ab38d7cd158a Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Mon, 10 Sep 2012 12:41:26 +0200 Subject: add SparseLU in sparse bench --- Eigen/src/OrderingMethods/Eigen_Colamd.h | 1091 ++++++++---------------------- Eigen/src/OrderingMethods/Ordering.h | 14 +- Eigen/src/SparseLU/SparseLU.h | 2 +- bench/spbench/CMakeLists.txt | 11 +- bench/spbench/spbenchsolver.h | 72 +- 5 files changed, 363 insertions(+), 827 deletions(-) (limited to 'bench/spbench/CMakeLists.txt') diff --git a/Eigen/src/OrderingMethods/Eigen_Colamd.h b/Eigen/src/OrderingMethods/Eigen_Colamd.h index 0af137d54..686c0f9f9 100644 --- a/Eigen/src/OrderingMethods/Eigen_Colamd.h +++ b/Eigen/src/OrderingMethods/Eigen_Colamd.h @@ -7,7 +7,7 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. -// This file is modified from the eigen_colamd/symamd library. The copyright is below +// This file is modified from the colamd/symamd library. The copyright is below // The authors of the code itself are Stefan I. Larimore and Timothy A. // Davis (davis@cise.ufl.edu), University of Florida. The algorithm was @@ -39,18 +39,19 @@ // // Availability: // -// The eigen_colamd/symamd library is available at +// The colamd/symamd library is available at // -// http://www.cise.ufl.edu/research/sparse/eigen_colamd/ +// http://www.cise.ufl.edu/research/sparse/colamd/ -// This is the http://www.cise.ufl.edu/research/sparse/eigen_colamd/eigen_colamd.h -// file. It is required by the eigen_colamd.c, colamdmex.c, and symamdmex.c +// This is the http://www.cise.ufl.edu/research/sparse/colamd/colamd.h +// file. It is required by the colamd.c, colamdmex.c, and symamdmex.c // files, and by any C code that calls the routines whose prototypes are -// listed below, or that uses the eigen_colamd/symamd definitions listed below. +// listed below, or that uses the colamd/symamd definitions listed below. #ifndef EIGEN_COLAMD_H #define EIGEN_COLAMD_H +namespace internal { /* Ensure that debugging is turned off: */ #ifndef COLAMD_NDEBUG #define COLAMD_NDEBUG @@ -61,42 +62,42 @@ /* ========================================================================== */ /* size of the knobs [ ] array. Only knobs [0..1] are currently used. */ -#define EIGEN_COLAMD_KNOBS 20 +#define COLAMD_KNOBS 20 /* number of output statistics. Only stats [0..6] are currently used. */ -#define EIGEN_COLAMD_STATS 20 +#define COLAMD_STATS 20 /* knobs [0] and stats [0]: dense row knob and output statistic. */ -#define EIGEN_COLAMD_DENSE_ROW 0 +#define COLAMD_DENSE_ROW 0 /* knobs [1] and stats [1]: dense column knob and output statistic. */ -#define EIGEN_COLAMD_DENSE_COL 1 +#define COLAMD_DENSE_COL 1 /* stats [2]: memory defragmentation count output statistic */ -#define EIGEN_COLAMD_DEFRAG_COUNT 2 +#define COLAMD_DEFRAG_COUNT 2 -/* stats [3]: eigen_colamd status: zero OK, > 0 warning or notice, < 0 error */ -#define EIGEN_COLAMD_STATUS 3 +/* stats [3]: colamd status: zero OK, > 0 warning or notice, < 0 error */ +#define COLAMD_STATUS 3 /* stats [4..6]: error info, or info on jumbled columns */ -#define EIGEN_COLAMD_INFO1 4 -#define EIGEN_COLAMD_INFO2 5 -#define EIGEN_COLAMD_INFO3 6 +#define COLAMD_INFO1 4 +#define COLAMD_INFO2 5 +#define COLAMD_INFO3 6 /* error codes returned in stats [3]: */ -#define EIGEN_COLAMD_OK (0) -#define EIGEN_COLAMD_OK_BUT_JUMBLED (1) -#define EIGEN_COLAMD_ERROR_A_not_present (-1) -#define EIGEN_COLAMD_ERROR_p_not_present (-2) -#define EIGEN_COLAMD_ERROR_nrow_negative (-3) -#define EIGEN_COLAMD_ERROR_ncol_negative (-4) -#define EIGEN_COLAMD_ERROR_nnz_negative (-5) -#define EIGEN_COLAMD_ERROR_p0_nonzero (-6) -#define EIGEN_COLAMD_ERROR_A_too_small (-7) -#define EIGEN_COLAMD_ERROR_col_length_negative (-8) -#define EIGEN_COLAMD_ERROR_row_index_out_of_bounds (-9) -#define EIGEN_COLAMD_ERROR_out_of_memory (-10) -#define EIGEN_COLAMD_ERROR_internal_error (-999) +#define COLAMD_OK (0) +#define COLAMD_OK_BUT_JUMBLED (1) +#define COLAMD_ERROR_A_not_present (-1) +#define COLAMD_ERROR_p_not_present (-2) +#define COLAMD_ERROR_nrow_negative (-3) +#define COLAMD_ERROR_ncol_negative (-4) +#define COLAMD_ERROR_nnz_negative (-5) +#define COLAMD_ERROR_p0_nonzero (-6) +#define COLAMD_ERROR_A_too_small (-7) +#define COLAMD_ERROR_col_length_negative (-8) +#define COLAMD_ERROR_row_index_out_of_bounds (-9) +#define COLAMD_ERROR_out_of_memory (-10) +#define COLAMD_ERROR_internal_error (-999) /* ========================================================================== */ /* === Definitions ========================================================== */ @@ -105,30 +106,30 @@ #define COLAMD_MAX(a,b) (((a) > (b)) ? (a) : (b)) #define COLAMD_MIN(a,b) (((a) < (b)) ? (a) : (b)) -#define EIGEN_ONES_COMPLEMENT(r) (-(r)-1) +#define ONES_COMPLEMENT(r) (-(r)-1) /* -------------------------------------------------------------------------- */ -#define EIGEN_COLAMD_EMPTY (-1) +#define COLAMD_EMPTY (-1) /* Row and column status */ -#define EIGEN_ALIVE (0) -#define EIGEN_DEAD (-1) +#define ALIVE (0) +#define DEAD (-1) /* Column status */ -#define EIGEN_DEAD_PRINCIPAL (-1) -#define EIGEN_DEAD_NON_PRINCIPAL (-2) +#define DEAD_PRINCIPAL (-1) +#define DEAD_NON_PRINCIPAL (-2) /* Macros for row and column status update and checking. */ -#define EIGEN_ROW_IS_DEAD(r) EIGEN_ROW_IS_MARKED_DEAD (Row[r].shared2.mark) -#define EIGEN_ROW_IS_MARKED_DEAD(row_mark) (row_mark < EIGEN_ALIVE) -#define EIGEN_ROW_IS_ALIVE(r) (Row [r].shared2.mark >= EIGEN_ALIVE) -#define EIGEN_COL_IS_DEAD(c) (Col [c].start < EIGEN_ALIVE) -#define EIGEN_COL_IS_ALIVE(c) (Col [c].start >= EIGEN_ALIVE) -#define EIGEN_EIGEN_COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == EIGEN_DEAD_PRINCIPAL) -#define EIGEN_KILL_ROW(r) { Row [r].shared2.mark = EIGEN_DEAD ; } -#define EIGEN_KILL_PRINCIPAL_COL(c) { Col [c].start = EIGEN_DEAD_PRINCIPAL ; } -#define EIGEN_KILL_NON_PRINCIPAL_COL(c) { Col [c].start = EIGEN_DEAD_NON_PRINCIPAL ; } +#define ROW_IS_DEAD(r) ROW_IS_MARKED_DEAD (Row[r].shared2.mark) +#define ROW_IS_MARKED_DEAD(row_mark) (row_mark < ALIVE) +#define ROW_IS_ALIVE(r) (Row [r].shared2.mark >= ALIVE) +#define COL_IS_DEAD(c) (Col [c].start < ALIVE) +#define COL_IS_ALIVE(c) (Col [c].start >= ALIVE) +#define COL_IS_DEAD_PRINCIPAL(c) (Col [c].start == DEAD_PRINCIPAL) +#define KILL_ROW(r) { Row [r].shared2.mark = DEAD ; } +#define KILL_PRINCIPAL_COL(c) { Col [c].start = DEAD_PRINCIPAL ; } +#define KILL_NON_PRINCIPAL_COL(c) { Col [c].start = DEAD_NON_PRINCIPAL ; } /* ========================================================================== */ /* === Colamd reporting mechanism =========================================== */ @@ -146,7 +147,7 @@ /* Use printf in standard C environment, for debugging and statistics output. */ /* Output is generated only if debugging is enabled at compile time, or if */ -/* the caller explicitly calls eigen_colamd_report or symamd_report. */ +/* the caller explicitly calls colamd_report or symamd_report. */ #define PRINTF printf /* In C, matrices are 0-based and indices are reported as such in *_report */ @@ -155,9 +156,9 @@ #endif /* MATLAB_MEX_FILE */ // == Row and Column structures == -typedef struct EIGEN_Colamd_Col_struct +typedef struct colamd_col_struct { - int start ; /* index for A of first row in this column, or EIGEN_DEAD */ + int start ; /* index for A of first row in this column, or DEAD */ /* if column is dead */ int length ; /* number of rows in this column */ union @@ -186,16 +187,16 @@ typedef struct EIGEN_Colamd_Col_struct int hash_next ; /* next column, if col is in a hash list */ } shared4 ; -} EIGEN_Colamd_Col ; +} colamd_col ; -typedef struct EIGEN_Colamd_Row_struct +typedef struct Colamd_Row_struct { int start ; /* index for A of first col in this row */ int length ; /* number of principal columns in this row */ union { int degree ; /* number of principal & non-principal columns in row */ - int p ; /* used as a row pointer in eigen_init_rows_cols () */ + int p ; /* used as a row pointer in init_rows_cols () */ } shared1 ; union { @@ -203,19 +204,19 @@ typedef struct EIGEN_Colamd_Row_struct int first_column ;/* first column in row (used in garbage collection) */ } shared2 ; -} EIGEN_Colamd_Row ; +} Colamd_Row ; /* ========================================================================== */ /* === Colamd recommended memory size ======================================= */ /* ========================================================================== */ /* - The recommended length Alen of the array A passed to eigen_colamd is given by - the EIGEN_COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any + The recommended length Alen of the array A passed to colamd is given by + the COLAMD_RECOMMENDED (nnz, n_row, n_col) macro. It returns -1 if any argument is negative. 2*nnz space is required for the row and column - indices of the matrix. EIGEN_COLAMD_C (n_col) + EIGEN_COLAMD_R (n_row) space is + indices of the matrix. colamd_c (n_col) + colamd_r (n_row) space is required for the Col and Row arrays, respectively, which are internal to - eigen_colamd. An additional n_col space is the minimal amount of "elbow room", + colamd. An additional n_col space is the minimal amount of "elbow room", and nnz/5 more space is recommended for run time efficiency. This macro is not needed when using symamd. @@ -224,120 +225,41 @@ typedef struct EIGEN_Colamd_Row_struct gcc -pedantic warning messages. */ -#define EIGEN_COLAMD_C(n_col) ((int) (((n_col) + 1) * sizeof (EIGEN_Colamd_Col) / sizeof (int))) -#define EIGEN_COLAMD_R(n_row) ((int) (((n_row) + 1) * sizeof (EIGEN_Colamd_Row) / sizeof (int))) +inline int colamd_c(int n_col) +{ return int( ((n_col) + 1) * sizeof (colamd_col) / sizeof (int) ) ; } -#define EIGEN_COLAMD_RECOMMENDED(nnz, n_row, n_col) \ -( \ -((nnz) < 0 || (n_row) < 0 || (n_col) < 0) \ -? \ - (-1) \ -: \ - (2 * (nnz) + EIGEN_COLAMD_C (n_col) + EIGEN_COLAMD_R (n_row) + (n_col) + ((nnz) / 5)) \ -) +inline int colamd_r(int n_row) +{ return int(((n_row) + 1) * sizeof (Colamd_Row) / sizeof (int)); } // Various routines -int eigen_colamd_recommended (int nnz, int n_row, int n_col) ; +inline int colamd_recommended (int nnz, int n_row, int n_col) ; -void eigen_colamd_set_defaults (double knobs [EIGEN_COLAMD_KNOBS]) ; +static inline void colamd_set_defaults (double knobs [COLAMD_KNOBS]) ; -bool eigen_colamd (int n_row, int n_col, int Alen, int A [], int p [], double knobs[EIGEN_COLAMD_KNOBS], int stats [EIGEN_COLAMD_STATS]) ; +static bool colamd (int n_row, int n_col, int Alen, int A [], int p [], double knobs[COLAMD_KNOBS], int stats [COLAMD_STATS]) ; -void eigen_colamd_report (int stats [EIGEN_COLAMD_STATS]); +static inline void colamd_report (int stats [COLAMD_STATS]); -int eigen_init_rows_cols (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col col [], int A [], int p [], int stats[EIGEN_COLAMD_STATS] ); +static int init_rows_cols (int n_row, int n_col, Colamd_Row Row [], colamd_col col [], int A [], int p [], int stats[COLAMD_STATS] ); -void eigen_init_scoring (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], double knobs[EIGEN_COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg); +static void init_scoring (int n_row, int n_col, Colamd_Row Row [], colamd_col Col [], int A [], int head [], double knobs[COLAMD_KNOBS], int *p_n_row2, int *p_n_col2, int *p_max_deg); -int eigen_find_ordering (int n_row, int n_col, int Alen, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int head [], int n_col2, int max_deg, int pfree); +static int find_ordering (int n_row, int n_col, int Alen, Colamd_Row Row [], colamd_col Col [], int A [], int head [], int n_col2, int max_deg, int pfree); -void eigen_order_children (int n_col, EIGEN_Colamd_Col Col [], int p []); +static void order_children (int n_col, colamd_col Col [], int p []); -void eigen_detect_super_cols ( -#ifndef COLAMD_NDEBUG - int n_col, - EIGEN_Colamd_Row Row [], -#endif /* COLAMD_NDEBUG */ - EIGEN_Colamd_Col Col [], +static void detect_super_cols ( + colamd_col Col [], int A [], int head [], int row_start, int row_length ) ; - int eigen_garbage_collection (int n_row, int n_col, EIGEN_Colamd_Row Row [], EIGEN_Colamd_Col Col [], int A [], int *pfree) ; - - int eigen_clear_mark (int n_row, EIGEN_Colamd_Row Row [] ) ; - - void eigen_print_report (char *method, int stats [EIGEN_COLAMD_STATS]) ; - -/* ========================================================================== */ -/* === Debugging prototypes and definitions ================================= */ -/* ========================================================================== */ - -#ifndef COLAMD_NDEBUG - -/* colamd_debug is the *ONLY* global variable, and is only */ -/* present when debugging */ - - int colamd_debug ; /* debug print level */ - -#define COLAMD_DEBUG0(params) { (void) PRINTF params ; } -#define COLAMD_DEBUG1(params) { if (colamd_debug >= 1) (void) PRINTF params ; } -#define COLAMD_DEBUG2(params) { if (colamd_debug >= 2) (void) PRINTF params ; } -#define COLAMD_DEBUG3(params) { if (colamd_debug >= 3) (void) PRINTF params ; } -#define COLAMD_DEBUG4(params) { if (colamd_debug >= 4) (void) PRINTF params ; } - -#ifdef MATLAB_MEX_FILE -#define COLAMD_ASSERT(expression) (mxAssert ((expression), "")) -#else -#define COLAMD_ASSERT(expression) (assert (expression)) -#endif /* MATLAB_MEX_FILE */ - - void eigen_colamd_get_debug /* gets the debug print level from getenv */ -( - char *method -) ; - - void eigen_debug_deg_lists -( - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int head [], - int min_score, - int should, - int max_deg -) ; - - void eigen_debug_mark -( - int n_row, - EIGEN_Colamd_Row Row [], - int tag_mark, - int max_mark -) ; +static int garbage_collection (int n_row, int n_col, Colamd_Row Row [], colamd_col Col [], int A [], int *pfree) ; - void eigen_debug_matrix -( - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [] -) ; - - void eigen_debug_structures -( - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [], - int n_col2 -) ; +static inline int clear_mark (int n_row, Colamd_Row Row [] ) ; -#else /* COLAMD_NDEBUG */ +static void print_report (const char *method, int stats [COLAMD_STATS]) ; /* === No debugging ========================================================= */ @@ -349,51 +271,50 @@ void eigen_detect_super_cols ( #define COLAMD_ASSERT(expression) ((void) 0) -#endif /* COLAMD_NDEBUG */ - - /** * \brief Returns the recommended value of Alen * - * Returns recommended value of Alen for use by eigen_colamd. + * Returns recommended value of Alen for use by colamd. * Returns -1 if any input argument is negative. * The use of this routine or macro is optional. * Note that the macro uses its arguments more than once, - * so be careful for side effects, if you pass expressions as arguments to EIGEN_COLAMD_RECOMMENDED. + * so be careful for side effects, if you pass expressions as arguments to COLAMD_RECOMMENDED. * * \param nnz nonzeros in A * \param n_row number of rows in A * \param n_col number of columns in A - * \return recommended value of Alen for use by eigen_colamd + * \return recommended value of Alen for use by colamd */ -int eigen_colamd_recommended ( int nnz, int n_row, int n_col) +inline int colamd_recommended ( int nnz, int n_row, int n_col) { - - return (EIGEN_COLAMD_RECOMMENDED (nnz, n_row, n_col)) ; + if ((nnz) < 0 || (n_row) < 0 || (n_col) < 0) + return (-1); + else + return (2 * (nnz) + colamd_c (n_col) + colamd_r (n_row) + (n_col) + ((nnz) / 5)); } /** * \brief set default parameters The use of this routine is optional. * - * Colamd: rows with more than (knobs [EIGEN_COLAMD_DENSE_ROW] * n_col) + * Colamd: rows with more than (knobs [COLAMD_DENSE_ROW] * n_col) * entries are removed prior to ordering. Columns with more than - * (knobs [EIGEN_COLAMD_DENSE_COL] * n_row) entries are removed prior to + * (knobs [COLAMD_DENSE_COL] * n_row) entries are removed prior to * ordering, and placed last in the output column ordering. * - * EIGEN_COLAMD_DENSE_ROW and EIGEN_COLAMD_DENSE_COL are defined as 0 and 1, - * respectively, in eigen_colamd.h. Default values of these two knobs + * COLAMD_DENSE_ROW and COLAMD_DENSE_COL are defined as 0 and 1, + * respectively, in colamd.h. Default values of these two knobs * are both 0.5. Currently, only knobs [0] and knobs [1] are * used, but future versions may use more knobs. If so, they will * be properly set to their defaults by the future version of - * eigen_colamd_set_defaults, so that the code that calls eigen_colamd will + * colamd_set_defaults, so that the code that calls colamd will * not need to change, assuming that you either use - * eigen_colamd_set_defaults, or pass a (double *) NULL pointer as the - * knobs array to eigen_colamd or symamd. + * colamd_set_defaults, or pass a (double *) NULL pointer as the + * knobs array to colamd or symamd. * - * \param knobs parameter settings for eigen_colamd + * \param knobs parameter settings for colamd */ -void eigen_colamd_set_defaults(double knobs[EIGEN_COLAMD_KNOBS]) +static inline void colamd_set_defaults(double knobs[COLAMD_KNOBS]) { /* === Local variables ================================================== */ @@ -403,12 +324,12 @@ void eigen_colamd_set_defaults(double knobs[EIGEN_COLAMD_KNOBS]) { return ; /* no knobs to initialize */ } - for (i = 0 ; i < EIGEN_COLAMD_KNOBS ; i++) + for (i = 0 ; i < COLAMD_KNOBS ; i++) { knobs [i] = 0 ; } - knobs [EIGEN_COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ - knobs [EIGEN_COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ + knobs [COLAMD_DENSE_ROW] = 0.5 ; /* ignore rows over 50% dense */ + knobs [COLAMD_DENSE_COL] = 0.5 ; /* ignore columns over 50% dense */ } /** @@ -425,10 +346,10 @@ void eigen_colamd_set_defaults(double knobs[EIGEN_COLAMD_KNOBS]) * \param Alen, size of the array A * \param A row indices of the matrix, of size ALen * \param p column pointers of A, of size n_col+1 - * \param knobs parameter settings for eigen_colamd - * \param stats eigen_colamd output statistics and error codes + * \param knobs parameter settings for colamd + * \param stats colamd output statistics and error codes */ -bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[EIGEN_COLAMD_KNOBS], int stats[EIGEN_COLAMD_STATS]) +static bool colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[COLAMD_KNOBS], int stats[COLAMD_STATS]) { /* === Local variables ================================================== */ @@ -437,77 +358,74 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E int Row_size ; /* size of Row [], in integers */ int Col_size ; /* size of Col [], in integers */ int need ; /* minimum required length of A */ - EIGEN_Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ - EIGEN_Colamd_Col *Col ; /* pointer into A of Col [0..n_col] array */ + Colamd_Row *Row ; /* pointer into A of Row [0..n_row] array */ + colamd_col *Col ; /* pointer into A of Col [0..n_col] array */ int n_col2 ; /* number of non-dense, non-empty columns */ int n_row2 ; /* number of non-dense, non-empty rows */ int ngarbage ; /* number of garbage collections performed */ int max_deg ; /* maximum row degree */ - double default_knobs [EIGEN_COLAMD_KNOBS] ; /* default knobs array */ + double default_knobs [COLAMD_KNOBS] ; /* default knobs array */ -#ifndef COLAMD_NDEBUG - eigen_colamd_get_debug ("eigen_colamd") ; -#endif /* COLAMD_NDEBUG */ /* === Check the input arguments ======================================== */ if (!stats) { - COLAMD_DEBUG0 (("eigen_colamd: stats not present\n")) ; + COLAMD_DEBUG0 (("colamd: stats not present\n")) ; return (false) ; } - for (i = 0 ; i < EIGEN_COLAMD_STATS ; i++) + for (i = 0 ; i < COLAMD_STATS ; i++) { stats [i] = 0 ; } - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_OK ; - stats [EIGEN_COLAMD_INFO1] = -1 ; - stats [EIGEN_COLAMD_INFO2] = -1 ; + stats [COLAMD_STATUS] = COLAMD_OK ; + stats [COLAMD_INFO1] = -1 ; + stats [COLAMD_INFO2] = -1 ; if (!A) /* A is not present */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_A_not_present ; - COLAMD_DEBUG0 (("eigen_colamd: A not present\n")) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_A_not_present ; + COLAMD_DEBUG0 (("colamd: A not present\n")) ; return (false) ; } if (!p) /* p is not present */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_p_not_present ; - COLAMD_DEBUG0 (("eigen_colamd: p not present\n")) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_p_not_present ; + COLAMD_DEBUG0 (("colamd: p not present\n")) ; return (false) ; } if (n_row < 0) /* n_row must be >= 0 */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_nrow_negative ; - stats [EIGEN_COLAMD_INFO1] = n_row ; - COLAMD_DEBUG0 (("eigen_colamd: nrow negative %d\n", n_row)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_nrow_negative ; + stats [COLAMD_INFO1] = n_row ; + COLAMD_DEBUG0 (("colamd: nrow negative %d\n", n_row)) ; return (false) ; } if (n_col < 0) /* n_col must be >= 0 */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_ncol_negative ; - stats [EIGEN_COLAMD_INFO1] = n_col ; - COLAMD_DEBUG0 (("eigen_colamd: ncol negative %d\n", n_col)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_ncol_negative ; + stats [COLAMD_INFO1] = n_col ; + COLAMD_DEBUG0 (("colamd: ncol negative %d\n", n_col)) ; return (false) ; } nnz = p [n_col] ; if (nnz < 0) /* nnz must be >= 0 */ { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_nnz_negative ; - stats [EIGEN_COLAMD_INFO1] = nnz ; - COLAMD_DEBUG0 (("eigen_colamd: number of entries negative %d\n", nnz)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_nnz_negative ; + stats [COLAMD_INFO1] = nnz ; + COLAMD_DEBUG0 (("colamd: number of entries negative %d\n", nnz)) ; return (false) ; } if (p [0] != 0) { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_p0_nonzero ; - stats [EIGEN_COLAMD_INFO1] = p [0] ; - COLAMD_DEBUG0 (("eigen_colamd: p[0] not zero %d\n", p [0])) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_p0_nonzero ; + stats [COLAMD_INFO1] = p [0] ; + COLAMD_DEBUG0 (("colamd: p[0] not zero %d\n", p [0])) ; return (false) ; } @@ -515,72 +433,73 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E if (!knobs) { - eigen_colamd_set_defaults (default_knobs) ; + colamd_set_defaults (default_knobs) ; knobs = default_knobs ; } /* === Allocate the Row and Col arrays from array A ===================== */ - Col_size = EIGEN_COLAMD_C (n_col) ; - Row_size = EIGEN_COLAMD_R (n_row) ; + Col_size = colamd_c (n_col) ; + Row_size = colamd_r (n_row) ; need = 2*nnz + n_col + Col_size + Row_size ; if (need > Alen) { /* not enough space in array A to perform the ordering */ - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_A_too_small ; - stats [EIGEN_COLAMD_INFO1] = need ; - stats [EIGEN_COLAMD_INFO2] = Alen ; - COLAMD_DEBUG0 (("eigen_colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); + stats [COLAMD_STATUS] = COLAMD_ERROR_A_too_small ; + stats [COLAMD_INFO1] = need ; + stats [COLAMD_INFO2] = Alen ; + COLAMD_DEBUG0 (("colamd: Need Alen >= %d, given only Alen = %d\n", need,Alen)); return (false) ; } Alen -= Col_size + Row_size ; - Col = (EIGEN_Colamd_Col *) &A [Alen] ; - Row = (EIGEN_Colamd_Row *) &A [Alen + Col_size] ; + Col = (colamd_col *) &A [Alen] ; + Row = (Colamd_Row *) &A [Alen + Col_size] ; /* === Construct the row and column data structures ===================== */ - if (!eigen_init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) + if (!init_rows_cols (n_row, n_col, Row, Col, A, p, stats)) { /* input matrix is invalid */ - COLAMD_DEBUG0 (("eigen_colamd: Matrix invalid\n")) ; + COLAMD_DEBUG0 (("colamd: Matrix invalid\n")) ; return (false) ; } /* === Initialize scores, kill dense rows/columns ======================= */ - eigen_init_scoring (n_row, n_col, Row, Col, A, p, knobs, + init_scoring (n_row, n_col, Row, Col, A, p, knobs, &n_row2, &n_col2, &max_deg) ; /* === Order the supercolumns =========================================== */ - ngarbage = eigen_find_ordering (n_row, n_col, Alen, Row, Col, A, p, + ngarbage = find_ordering (n_row, n_col, Alen, Row, Col, A, p, n_col2, max_deg, 2*nnz) ; /* === Order the non-principal columns ================================== */ - eigen_order_children (n_col, Col, p) ; + order_children (n_col, Col, p) ; /* === Return statistics in stats ======================================= */ - stats [EIGEN_COLAMD_DENSE_ROW] = n_row - n_row2 ; - stats [EIGEN_COLAMD_DENSE_COL] = n_col - n_col2 ; - stats [EIGEN_COLAMD_DEFRAG_COUNT] = ngarbage ; - COLAMD_DEBUG0 (("eigen_colamd: done.\n")) ; + stats [COLAMD_DENSE_ROW] = n_row - n_row2 ; + stats [COLAMD_DENSE_COL] = n_col - n_col2 ; + stats [COLAMD_DEFRAG_COUNT] = ngarbage ; + COLAMD_DEBUG0 (("colamd: done.\n")) ; return (true) ; } /* ========================================================================== */ -/* === eigen_colamd_report ======================================================== */ +/* === colamd_report ======================================================== */ /* ========================================================================== */ - void eigen_colamd_report + static inline void colamd_report ( - int stats [EIGEN_COLAMD_STATS] + int stats [COLAMD_STATS] ) { - eigen_print_report ("eigen_colamd", stats) ; + const char *method = "colamd"; + print_report (method, stats) ; } @@ -592,7 +511,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_init_rows_cols ======================================================= */ +/* === init_rows_cols ======================================================= */ /* ========================================================================== */ /* @@ -604,17 +523,17 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E true otherwise. Not user-callable. */ - int eigen_init_rows_cols /* returns true if OK, or false otherwise */ + static int init_rows_cols /* returns true if OK, or false otherwise */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + Colamd_Row Row [], /* of size n_row+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* row indices of A, of size Alen */ int p [], /* pointers to columns in A, of size n_col+1 */ - int stats [EIGEN_COLAMD_STATS] /* eigen_colamd statistics */ + int stats [COLAMD_STATS] /* colamd statistics */ ) { /* === Local variables ================================================== */ @@ -637,24 +556,24 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E if (Col [col].length < 0) { /* column pointers must be non-decreasing */ - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_col_length_negative ; - stats [EIGEN_COLAMD_INFO1] = col ; - stats [EIGEN_COLAMD_INFO2] = Col [col].length ; - COLAMD_DEBUG0 (("eigen_colamd: col %d length %d < 0\n", col, Col [col].length)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_col_length_negative ; + stats [COLAMD_INFO1] = col ; + stats [COLAMD_INFO2] = Col [col].length ; + COLAMD_DEBUG0 (("colamd: col %d length %d < 0\n", col, Col [col].length)) ; return (false) ; } Col [col].shared1.thickness = 1 ; Col [col].shared2.score = 0 ; - Col [col].shared3.prev = EIGEN_COLAMD_EMPTY ; - Col [col].shared4.degree_next = EIGEN_COLAMD_EMPTY ; + Col [col].shared3.prev = COLAMD_EMPTY ; + Col [col].shared4.degree_next = COLAMD_EMPTY ; } /* p [0..n_col] no longer needed, used as "head" in subsequent routines */ /* === Scan columns, compute row degrees, and check row indices ========= */ - stats [EIGEN_COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ + stats [COLAMD_INFO3] = 0 ; /* number of duplicate or unsorted row indices*/ for (row = 0 ; row < n_row ; row++) { @@ -676,11 +595,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* make sure row indices within range */ if (row < 0 || row >= n_row) { - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_ERROR_row_index_out_of_bounds ; - stats [EIGEN_COLAMD_INFO1] = col ; - stats [EIGEN_COLAMD_INFO2] = row ; - stats [EIGEN_COLAMD_INFO3] = n_row ; - COLAMD_DEBUG0 (("eigen_colamd: row %d col %d out of bounds\n", row, col)) ; + stats [COLAMD_STATUS] = COLAMD_ERROR_row_index_out_of_bounds ; + stats [COLAMD_INFO1] = col ; + stats [COLAMD_INFO2] = row ; + stats [COLAMD_INFO3] = n_row ; + COLAMD_DEBUG0 (("colamd: row %d col %d out of bounds\n", row, col)) ; return (false) ; } @@ -688,11 +607,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { /* row index are unsorted or repeated (or both), thus col */ /* is jumbled. This is a notice, not an error condition. */ - stats [EIGEN_COLAMD_STATUS] = EIGEN_COLAMD_OK_BUT_JUMBLED ; - stats [EIGEN_COLAMD_INFO1] = col ; - stats [EIGEN_COLAMD_INFO2] = row ; - (stats [EIGEN_COLAMD_INFO3]) ++ ; - COLAMD_DEBUG1 (("eigen_colamd: row %d col %d unsorted/duplicate\n",row,col)); + stats [COLAMD_STATUS] = COLAMD_OK_BUT_JUMBLED ; + stats [COLAMD_INFO1] = col ; + stats [COLAMD_INFO2] = row ; + (stats [COLAMD_INFO3]) ++ ; + COLAMD_DEBUG1 (("colamd: row %d col %d unsorted/duplicate\n",row,col)); } if (Row [row].shared2.mark != col) @@ -729,7 +648,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* === Create row form ================================================== */ - if (stats [EIGEN_COLAMD_STATUS] == EIGEN_COLAMD_OK_BUT_JUMBLED) + if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) { /* if cols jumbled, watch for repeated row indices */ for (col = 0 ; col < n_col ; col++) @@ -771,31 +690,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* === See if we need to re-create columns ============================== */ - if (stats [EIGEN_COLAMD_STATUS] == EIGEN_COLAMD_OK_BUT_JUMBLED) + if (stats [COLAMD_STATUS] == COLAMD_OK_BUT_JUMBLED) { - COLAMD_DEBUG0 (("eigen_colamd: reconstructing column form, matrix jumbled\n")) ; + COLAMD_DEBUG0 (("colamd: reconstructing column form, matrix jumbled\n")) ; -#ifndef COLAMD_NDEBUG - /* make sure column lengths are correct */ - for (col = 0 ; col < n_col ; col++) - { - p [col] = Col [col].length ; - } - for (row = 0 ; row < n_row ; row++) - { - rp = &A [Row [row].start] ; - rp_end = rp + Row [row].length ; - while (rp < rp_end) - { - p [*rp++]-- ; - } - } - for (col = 0 ; col < n_col ; col++) - { - COLAMD_ASSERT (p [col] == 0) ; - } - /* now p is all zero (different than when debugging is turned off) */ -#endif /* COLAMD_NDEBUG */ /* === Compute col pointers ========================================= */ @@ -833,7 +731,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_init_scoring ========================================================= */ +/* === init_scoring ========================================================= */ /* ========================================================================== */ /* @@ -841,17 +739,17 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E each column, and places all columns in the degree lists. Not user-callable. */ - void eigen_init_scoring +static void init_scoring ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + Colamd_Row Row [], /* of size n_row+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* column form and row form of A */ int head [], /* of size n_col+1 */ - double knobs [EIGEN_COLAMD_KNOBS],/* parameters */ + double knobs [COLAMD_KNOBS],/* parameters */ int *p_n_row2, /* number of non-dense, non-empty rows */ int *p_n_col2, /* number of non-dense, non-empty columns */ int *p_max_deg /* maximum row degree */ @@ -875,15 +773,12 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E int max_deg ; /* maximum row degree */ int next_col ; /* Used to add to degree list.*/ -#ifndef COLAMD_NDEBUG - int debug_count ; /* debug only. */ -#endif /* COLAMD_NDEBUG */ /* === Extract knobs ==================================================== */ - dense_row_count = COLAMD_MAX (0, COLAMD_MIN (knobs [EIGEN_COLAMD_DENSE_ROW] * n_col, n_col)) ; - dense_col_count = COLAMD_MAX (0, COLAMD_MIN (knobs [EIGEN_COLAMD_DENSE_COL] * n_row, n_row)) ; - COLAMD_DEBUG1 (("eigen_colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; + dense_row_count = COLAMD_MAX (0, COLAMD_MIN (knobs [COLAMD_DENSE_ROW] * n_col, n_col)) ; + dense_col_count = COLAMD_MAX (0, COLAMD_MIN (knobs [COLAMD_DENSE_COL] * n_row, n_row)) ; + COLAMD_DEBUG1 (("colamd: densecount: %d %d\n", dense_row_count, dense_col_count)) ; max_deg = 0 ; n_col2 = n_col ; n_row2 = n_row ; @@ -899,10 +794,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { /* this is a empty column, kill and order it last */ Col [c].shared2.order = --n_col2 ; - EIGEN_KILL_PRINCIPAL_COL (c) ; + KILL_PRINCIPAL_COL (c) ; } } - COLAMD_DEBUG1 (("eigen_colamd: null columns killed: %d\n", n_col - n_col2)) ; + COLAMD_DEBUG1 (("colamd: null columns killed: %d\n", n_col - n_col2)) ; /* === Kill dense columns =============================================== */ @@ -910,7 +805,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (c = n_col-1 ; c >= 0 ; c--) { /* skip any dead columns */ - if (EIGEN_COL_IS_DEAD (c)) + if (COL_IS_DEAD (c)) { continue ; } @@ -926,10 +821,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { Row [*cp++].shared1.degree-- ; } - EIGEN_KILL_PRINCIPAL_COL (c) ; + KILL_PRINCIPAL_COL (c) ; } } - COLAMD_DEBUG1 (("eigen_colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; + COLAMD_DEBUG1 (("colamd: Dense and null columns killed: %d\n", n_col - n_col2)) ; /* === Kill dense and empty rows ======================================== */ @@ -940,7 +835,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E if (deg > dense_row_count || deg == 0) { /* kill a dense or empty row */ - EIGEN_KILL_ROW (r) ; + KILL_ROW (r) ; --n_row2 ; } else @@ -949,7 +844,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E max_deg = COLAMD_MAX (max_deg, deg) ; } } - COLAMD_DEBUG1 (("eigen_colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; + COLAMD_DEBUG1 (("colamd: Dense and null rows killed: %d\n", n_row - n_row2)) ; /* === Compute initial column scores ==================================== */ @@ -962,7 +857,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (c = n_col-1 ; c >= 0 ; c--) { /* skip dead column */ - if (EIGEN_COL_IS_DEAD (c)) + if (COL_IS_DEAD (c)) { continue ; } @@ -975,7 +870,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* get a row */ row = *cp++ ; /* skip if dead */ - if (EIGEN_ROW_IS_DEAD (row)) + if (ROW_IS_DEAD (row)) { continue ; } @@ -994,7 +889,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* and have already been killed) */ COLAMD_DEBUG2 (("Newly null killed: %d\n", c)) ; Col [c].shared2.order = --n_col2 ; - EIGEN_KILL_PRINCIPAL_COL (c) ; + KILL_PRINCIPAL_COL (c) ; } else { @@ -1005,7 +900,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E Col [c].shared2.score = score ; } } - COLAMD_DEBUG1 (("eigen_colamd: Dense, null, and newly-null columns killed: %d\n", + COLAMD_DEBUG1 (("colamd: Dense, null, and newly-null columns killed: %d\n", n_col-n_col2)) ; /* At this point, all empty rows and columns are dead. All live columns */ @@ -1013,20 +908,13 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* yet). Rows may contain dead columns, but all live rows contain at */ /* least one live column. */ -#ifndef COLAMD_NDEBUG - eigen_debug_structures (n_row, n_col, Row, Col, A, n_col2) ; -#endif /* COLAMD_NDEBUG */ - /* === Initialize degree lists ========================================== */ -#ifndef COLAMD_NDEBUG - debug_count = 0 ; -#endif /* COLAMD_NDEBUG */ /* clear the hash buckets */ for (c = 0 ; c <= n_col ; c++) { - head [c] = EIGEN_COLAMD_EMPTY ; + head [c] = COLAMD_EMPTY ; } min_score = n_col ; /* place in reverse order, so low column indices are at the front */ @@ -1034,7 +922,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (c = n_col-1 ; c >= 0 ; c--) { /* only add principal columns to degree lists */ - if (EIGEN_COL_IS_ALIVE (c)) + if (COL_IS_ALIVE (c)) { COLAMD_DEBUG4 (("place %d score %d minscore %d ncol %d\n", c, Col [c].shared2.score, min_score, n_col)) ; @@ -1047,16 +935,16 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E COLAMD_ASSERT (min_score <= n_col) ; COLAMD_ASSERT (score >= 0) ; COLAMD_ASSERT (score <= n_col) ; - COLAMD_ASSERT (head [score] >= EIGEN_COLAMD_EMPTY) ; + COLAMD_ASSERT (head [score] >= COLAMD_EMPTY) ; /* now add this column to dList at proper score location */ next_col = head [score] ; - Col [c].shared3.prev = EIGEN_COLAMD_EMPTY ; + Col [c].shared3.prev = COLAMD_EMPTY ; Col [c].shared4.degree_next = next_col ; /* if there already was a column with the same score, set its */ /* previous pointer to this new column */ - if (next_col != EIGEN_COLAMD_EMPTY) + if (next_col != COLAMD_EMPTY) { Col [next_col].shared3.prev = c ; } @@ -1065,19 +953,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* see if this score is less than current min */ min_score = COLAMD_MIN (min_score, score) ; -#ifndef COLAMD_NDEBUG - debug_count++ ; -#endif /* COLAMD_NDEBUG */ } } -#ifndef COLAMD_NDEBUG - COLAMD_DEBUG1 (("eigen_colamd: Live cols %d out of %d, non-princ: %d\n", - debug_count, n_col, n_col-debug_count)) ; - COLAMD_ASSERT (debug_count == n_col2) ; - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, min_score, n_col2, max_deg) ; -#endif /* COLAMD_NDEBUG */ /* === Return number of remaining columns, and max row degree =========== */ @@ -1088,7 +967,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_find_ordering ======================================================== */ +/* === find_ordering ======================================================== */ /* ========================================================================== */ /* @@ -1097,15 +976,15 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E degree ordering method. Not user-callable. */ - int eigen_find_ordering /* return the number of garbage collections */ +static int find_ordering /* return the number of garbage collections */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows of A */ int n_col, /* number of columns of A */ int Alen, /* size of A, 2*nnz + n_col or larger */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + Colamd_Row Row [], /* of size n_row+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* column form and row form of A */ int head [], /* of size n_col+1 */ int n_col2, /* Remaining columns to order */ @@ -1147,55 +1026,29 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E int next_col ; /* Used by Dlist operations. */ int ngarbage ; /* number of garbage collections performed */ -#ifndef COLAMD_NDEBUG - int debug_d ; /* debug loop counter */ - int debug_step = 0 ; /* debug loop counter */ -#endif /* COLAMD_NDEBUG */ /* === Initialization and clear mark ==================================== */ max_mark = INT_MAX - n_col ; /* INT_MAX defined in */ - tag_mark = eigen_clear_mark (n_row, Row) ; + tag_mark = clear_mark (n_row, Row) ; min_score = 0 ; ngarbage = 0 ; - COLAMD_DEBUG1 (("eigen_colamd: Ordering, n_col2=%d\n", n_col2)) ; + COLAMD_DEBUG1 (("colamd: Ordering, n_col2=%d\n", n_col2)) ; /* === Order the columns ================================================ */ for (k = 0 ; k < n_col2 ; /* 'k' is incremented below */) { -#ifndef COLAMD_NDEBUG - if (debug_step % 100 == 0) - { - COLAMD_DEBUG2 (("\n... Step k: %d out of n_col2: %d\n", k, n_col2)) ; - } - else - { - COLAMD_DEBUG3 (("\n----------Step k: %d out of n_col2: %d\n", k, n_col2)) ; - } - debug_step++ ; - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, - min_score, n_col2-k, max_deg) ; - eigen_debug_matrix (n_row, n_col, Row, Col, A) ; -#endif /* COLAMD_NDEBUG */ - /* === Select pivot column, and order it ============================ */ /* make sure degree list isn't empty */ COLAMD_ASSERT (min_score >= 0) ; COLAMD_ASSERT (min_score <= n_col) ; - COLAMD_ASSERT (head [min_score] >= EIGEN_COLAMD_EMPTY) ; - -#ifndef COLAMD_NDEBUG - for (debug_d = 0 ; debug_d < min_score ; debug_d++) - { - COLAMD_ASSERT (head [debug_d] == EIGEN_COLAMD_EMPTY) ; - } -#endif /* COLAMD_NDEBUG */ + COLAMD_ASSERT (head [min_score] >= COLAMD_EMPTY) ; /* get pivot column from head of minimum degree list */ - while (head [min_score] == EIGEN_COLAMD_EMPTY && min_score < n_col) + while (head [min_score] == COLAMD_EMPTY && min_score < n_col) { min_score++ ; } @@ -1203,12 +1056,12 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E COLAMD_ASSERT (pivot_col >= 0 && pivot_col <= n_col) ; next_col = Col [pivot_col].shared4.degree_next ; head [min_score] = next_col ; - if (next_col != EIGEN_COLAMD_EMPTY) + if (next_col != COLAMD_EMPTY) { - Col [next_col].shared3.prev = EIGEN_COLAMD_EMPTY ; + Col [next_col].shared3.prev = COLAMD_EMPTY ; } - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (pivot_col)) ; + COLAMD_ASSERT (COL_IS_ALIVE (pivot_col)) ; COLAMD_DEBUG3 (("Pivot col: %d\n", pivot_col)) ; /* remember score for defrag check */ @@ -1227,16 +1080,13 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E needed_memory = COLAMD_MIN (pivot_col_score, n_col - k) ; if (pfree + needed_memory >= Alen) { - pfree = eigen_garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; + pfree = garbage_collection (n_row, n_col, Row, Col, A, &A [pfree]) ; ngarbage++ ; /* after garbage collection we will have enough */ COLAMD_ASSERT (pfree + needed_memory < Alen) ; /* garbage collection has wiped out the Row[].shared2.mark array */ - tag_mark = eigen_clear_mark (n_row, Row) ; + tag_mark = clear_mark (n_row, Row) ; -#ifndef COLAMD_NDEBUG - eigen_debug_matrix (n_row, n_col, Row, Col, A) ; -#endif /* COLAMD_NDEBUG */ } /* === Compute pivot row pattern ==================================== */ @@ -1258,9 +1108,9 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { /* get a row */ row = *cp++ ; - COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", EIGEN_ROW_IS_ALIVE (row), row)) ; + COLAMD_DEBUG4 (("Pivot col pattern %d %d\n", ROW_IS_ALIVE (row), row)) ; /* skip if row is dead */ - if (EIGEN_ROW_IS_DEAD (row)) + if (ROW_IS_DEAD (row)) { continue ; } @@ -1272,7 +1122,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E col = *rp++ ; /* add the column, if alive and untagged */ col_thickness = Col [col].shared1.thickness ; - if (col_thickness > 0 && EIGEN_COL_IS_ALIVE (col)) + if (col_thickness > 0 && COL_IS_ALIVE (col)) { /* tag column in pivot row */ Col [col].shared1.thickness = -col_thickness ; @@ -1288,10 +1138,6 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E Col [pivot_col].shared1.thickness = pivot_col_thickness ; max_deg = COLAMD_MAX (max_deg, pivot_row_degree) ; -#ifndef COLAMD_NDEBUG - COLAMD_DEBUG3 (("check2\n")) ; - eigen_debug_mark (n_row, Row, tag_mark, max_mark) ; -#endif /* COLAMD_NDEBUG */ /* === Kill all rows used to construct pivot row ==================== */ @@ -1303,7 +1149,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* may be killing an already dead row */ row = *cp++ ; COLAMD_DEBUG3 (("Kill row in pivot col: %d\n", row)) ; - EIGEN_KILL_ROW (row) ; + KILL_ROW (row) ; } /* === Select a row index to use as the new pivot row =============== */ @@ -1318,7 +1164,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E else { /* there is no pivot row, since it is of zero length */ - pivot_row = EIGEN_COLAMD_EMPTY ; + pivot_row = COLAMD_EMPTY ; COLAMD_ASSERT (pivot_row_length == 0) ; } COLAMD_ASSERT (Col [pivot_col].length > 0 || pivot_row_length == 0) ; @@ -1340,7 +1186,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* context, is the column "length", or the number of row indices */ /* in that column). The number of row indices in a column is */ /* monotonically non-decreasing, from the length of the original */ - /* column on input to eigen_colamd. */ + /* column on input to colamd. */ /* === Compute set differences ====================================== */ @@ -1355,7 +1201,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E while (rp < rp_end) { col = *rp++ ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col) && col != pivot_col) ; + COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; COLAMD_DEBUG3 (("Col: %d\n", col)) ; /* clear tags used to construct pivot row pattern */ @@ -1370,8 +1216,8 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E next_col = Col [col].shared4.degree_next ; COLAMD_ASSERT (cur_score >= 0) ; COLAMD_ASSERT (cur_score <= n_col) ; - COLAMD_ASSERT (cur_score >= EIGEN_COLAMD_EMPTY) ; - if (prev_col == EIGEN_COLAMD_EMPTY) + COLAMD_ASSERT (cur_score >= COLAMD_EMPTY) ; + if (prev_col == COLAMD_EMPTY) { head [cur_score] = next_col ; } @@ -1379,7 +1225,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { Col [prev_col].shared4.degree_next = next_col ; } - if (next_col != EIGEN_COLAMD_EMPTY) + if (next_col != COLAMD_EMPTY) { Col [next_col].shared3.prev = prev_col ; } @@ -1394,7 +1240,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E row = *cp++ ; row_mark = Row [row].shared2.mark ; /* skip if dead */ - if (EIGEN_ROW_IS_MARKED_DEAD (row_mark)) + if (ROW_IS_MARKED_DEAD (row_mark)) { continue ; } @@ -1413,7 +1259,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E if (set_difference == 0) { COLAMD_DEBUG3 (("aggressive absorption. Row: %d\n", row)) ; - EIGEN_KILL_ROW (row) ; + KILL_ROW (row) ; } else { @@ -1423,10 +1269,6 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E } } -#ifndef COLAMD_NDEBUG - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, - min_score, n_col2-k-pivot_row_degree, max_deg) ; -#endif /* COLAMD_NDEBUG */ /* === Add up set differences for each column ======================= */ @@ -1439,7 +1281,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { /* get a column */ col = *rp++ ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col) && col != pivot_col) ; + COLAMD_ASSERT (COL_IS_ALIVE (col) && col != pivot_col) ; hash = 0 ; cur_score = 0 ; cp = &A [Col [col].start] ; @@ -1456,7 +1298,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E COLAMD_ASSERT(row >= 0 && row < n_row) ; row_mark = Row [row].shared2.mark ; /* skip if dead */ - if (EIGEN_ROW_IS_MARKED_DEAD (row_mark)) + if (ROW_IS_MARKED_DEAD (row_mark)) { continue ; } @@ -1480,7 +1322,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { COLAMD_DEBUG4 (("further mass elimination. Col: %d\n", col)) ; /* nothing left but the pivot row in this column */ - EIGEN_KILL_PRINCIPAL_COL (col) ; + KILL_PRINCIPAL_COL (col) ; pivot_row_degree -= Col [col].shared1.thickness ; COLAMD_ASSERT (pivot_row_degree >= 0) ; /* order it */ @@ -1504,7 +1346,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E COLAMD_ASSERT (hash <= n_col) ; head_column = head [hash] ; - if (head_column > EIGEN_COLAMD_EMPTY) + if (head_column > COLAMD_EMPTY) { /* degree list "hash" is non-empty, use prev (shared3) of */ /* first column in degree list as head of hash bucket */ @@ -1521,7 +1363,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* save hash function in Col [col].shared3.hash */ Col [col].shared3.hash = (int) hash ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col)) ; + COLAMD_ASSERT (COL_IS_ALIVE (col)) ; } } @@ -1531,17 +1373,13 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E COLAMD_DEBUG3 (("** Supercolumn detection phase. **\n")) ; - eigen_detect_super_cols ( - -#ifndef COLAMD_NDEBUG - n_col, Row, -#endif /* COLAMD_NDEBUG */ + detect_super_cols ( Col, A, head, pivot_row_start, pivot_row_length) ; /* === Kill the pivotal column ====================================== */ - EIGEN_KILL_PRINCIPAL_COL (pivot_col) ; + KILL_PRINCIPAL_COL (pivot_col) ; /* === Clear mark =================================================== */ @@ -1549,14 +1387,9 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E if (tag_mark >= max_mark) { COLAMD_DEBUG2 (("clearing tag_mark\n")) ; - tag_mark = eigen_clear_mark (n_row, Row) ; + tag_mark = clear_mark (n_row, Row) ; } -#ifndef COLAMD_NDEBUG - COLAMD_DEBUG3 (("check3\n")) ; - eigen_debug_mark (n_row, Row, tag_mark, max_mark) ; -#endif /* COLAMD_NDEBUG */ - /* === Finalize the new pivot row, and column scores ================ */ COLAMD_DEBUG3 (("** Finalize scores phase. **\n")) ; @@ -1570,7 +1403,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { col = *rp++ ; /* skip dead columns */ - if (EIGEN_COL_IS_DEAD (col)) + if (COL_IS_DEAD (col)) { continue ; } @@ -1604,11 +1437,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E COLAMD_ASSERT (min_score <= n_col) ; COLAMD_ASSERT (cur_score >= 0) ; COLAMD_ASSERT (cur_score <= n_col) ; - COLAMD_ASSERT (head [cur_score] >= EIGEN_COLAMD_EMPTY) ; + COLAMD_ASSERT (head [cur_score] >= COLAMD_EMPTY) ; next_col = head [cur_score] ; Col [col].shared4.degree_next = next_col ; - Col [col].shared3.prev = EIGEN_COLAMD_EMPTY ; - if (next_col != EIGEN_COLAMD_EMPTY) + Col [col].shared3.prev = COLAMD_EMPTY ; + if (next_col != COLAMD_EMPTY) { Col [next_col].shared3.prev = col ; } @@ -1619,11 +1452,6 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E } -#ifndef COLAMD_NDEBUG - eigen_debug_deg_lists (n_row, n_col, Row, Col, head, - min_score, n_col2-k, max_deg) ; -#endif /* COLAMD_NDEBUG */ - /* === Resurrect the new pivot row ================================== */ if (pivot_row_degree > 0) @@ -1645,11 +1473,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_order_children ======================================================= */ +/* === order_children ======================================================= */ /* ========================================================================== */ /* - The eigen_find_ordering routine has ordered all of the principal columns (the + The find_ordering routine has ordered all of the principal columns (the representatives of the supercolumns). The non-principal columns have not yet been ordered. This routine orders those columns by walking up the parent tree (a column is a child of the column which absorbed it). The @@ -1661,12 +1489,12 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E columns. Not user-callable. */ - void eigen_order_children +static inline void order_children ( /* === Parameters ======================================================= */ int n_col, /* number of columns of A */ - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + colamd_col Col [], /* of size n_col+1 */ int p [] /* p [0 ... n_col-1] is the column permutation*/ ) { @@ -1682,15 +1510,15 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (i = 0 ; i < n_col ; i++) { /* find an un-ordered non-principal column */ - COLAMD_ASSERT (EIGEN_COL_IS_DEAD (i)) ; - if (!EIGEN_EIGEN_COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == EIGEN_COLAMD_EMPTY) + COLAMD_ASSERT (COL_IS_DEAD (i)) ; + if (!COL_IS_DEAD_PRINCIPAL (i) && Col [i].shared2.order == COLAMD_EMPTY) { parent = i ; /* once found, find its principal parent */ do { parent = Col [parent].shared1.parent ; - } while (!EIGEN_EIGEN_COL_IS_DEAD_PRINCIPAL (parent)) ; + } while (!COL_IS_DEAD_PRINCIPAL (parent)) ; /* now, order all un-ordered non-principal columns along path */ /* to this parent. collapse tree at the same time */ @@ -1700,7 +1528,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E do { - COLAMD_ASSERT (Col [c].shared2.order == EIGEN_COLAMD_EMPTY) ; + COLAMD_ASSERT (Col [c].shared2.order == COLAMD_EMPTY) ; /* order this column */ Col [c].shared2.order = order++ ; @@ -1713,7 +1541,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* continue until we hit an ordered column. There are */ /* guarranteed not to be anymore unordered columns */ /* above an ordered column */ - } while (Col [c].shared2.order == EIGEN_COLAMD_EMPTY) ; + } while (Col [c].shared2.order == COLAMD_EMPTY) ; /* re-order the super_col parent to largest order for this group */ Col [parent].shared2.order = order ; @@ -1730,7 +1558,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_detect_super_cols ==================================================== */ +/* === detect_super_cols ==================================================== */ /* ========================================================================== */ /* @@ -1762,17 +1590,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E Not user-callable. */ - void eigen_detect_super_cols +static void detect_super_cols ( /* === Parameters ======================================================= */ -#ifndef COLAMD_NDEBUG - /* these two parameters are only needed when debugging is enabled: */ - int n_col, /* number of columns of A */ - EIGEN_Colamd_Row Row [], /* of size n_row+1 */ -#endif /* COLAMD_NDEBUG */ - - EIGEN_Colamd_Col Col [], /* of size n_col+1 */ + colamd_col Col [], /* of size n_col+1 */ int A [], /* row indices of A */ int head [], /* head of degree lists and hash buckets */ int row_start, /* pointer to set of columns to check */ @@ -1802,7 +1624,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E while (rp < rp_end) { col = *rp++ ; - if (EIGEN_COL_IS_DEAD (col)) + if (COL_IS_DEAD (col)) { continue ; } @@ -1814,7 +1636,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* === Get the first column in this hash bucket ===================== */ head_column = head [hash] ; - if (head_column > EIGEN_COLAMD_EMPTY) + if (head_column > COLAMD_EMPTY) { first_col = Col [head_column].shared3.headhash ; } @@ -1825,10 +1647,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* === Consider each column in the hash bucket ====================== */ - for (super_c = first_col ; super_c != EIGEN_COLAMD_EMPTY ; + for (super_c = first_col ; super_c != COLAMD_EMPTY ; super_c = Col [super_c].shared4.hash_next) { - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (super_c)) ; + COLAMD_ASSERT (COL_IS_ALIVE (super_c)) ; COLAMD_ASSERT (Col [super_c].shared3.hash == hash) ; length = Col [super_c].length ; @@ -1838,10 +1660,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* === Compare super_c with all columns after it ================ */ for (c = Col [super_c].shared4.hash_next ; - c != EIGEN_COLAMD_EMPTY ; c = Col [c].shared4.hash_next) + c != COLAMD_EMPTY ; c = Col [c].shared4.hash_next) { COLAMD_ASSERT (c != super_c) ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (c)) ; + COLAMD_ASSERT (COL_IS_ALIVE (c)) ; COLAMD_ASSERT (Col [c].shared3.hash == hash) ; /* not identical if lengths or scores are different */ @@ -1859,8 +1681,8 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (i = 0 ; i < length ; i++) { /* the columns are "clean" (no dead rows) */ - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (*cp1)) ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (*cp2)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (*cp1)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (*cp2)) ; /* row indices will same order for both supercols, */ /* no gather scatter nessasary */ if (*cp1++ != *cp2++) @@ -1882,9 +1704,9 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E Col [super_c].shared1.thickness += Col [c].shared1.thickness ; Col [c].shared1.parent = super_c ; - EIGEN_KILL_NON_PRINCIPAL_COL (c) ; - /* order c later, in eigen_order_children() */ - Col [c].shared2.order = EIGEN_COLAMD_EMPTY ; + KILL_NON_PRINCIPAL_COL (c) ; + /* order c later, in order_children() */ + Col [c].shared2.order = COLAMD_EMPTY ; /* remove c from hash bucket */ Col [prev_c].shared4.hash_next = Col [c].shared4.hash_next ; } @@ -1892,22 +1714,22 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* === Empty this hash bucket ======================================= */ - if (head_column > EIGEN_COLAMD_EMPTY) + if (head_column > COLAMD_EMPTY) { /* corresponding degree list "hash" is not empty */ - Col [head_column].shared3.headhash = EIGEN_COLAMD_EMPTY ; + Col [head_column].shared3.headhash = COLAMD_EMPTY ; } else { /* corresponding degree list "hash" is empty */ - head [hash] = EIGEN_COLAMD_EMPTY ; + head [hash] = COLAMD_EMPTY ; } } } /* ========================================================================== */ -/* === eigen_garbage_collection =================================================== */ +/* === garbage_collection =================================================== */ /* ========================================================================== */ /* @@ -1919,14 +1741,14 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E Not user-callable. */ - int eigen_garbage_collection /* returns the new value of pfree */ +static int garbage_collection /* returns the new value of pfree */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows */ int n_col, /* number of columns */ - EIGEN_Colamd_Row Row [], /* row info */ - EIGEN_Colamd_Col Col [], /* column info */ + Colamd_Row Row [], /* row info */ + colamd_col Col [], /* column info */ int A [], /* A [0 ... Alen-1] holds the matrix */ int *pfree /* &A [0] ... pfree is in use */ ) @@ -1940,19 +1762,12 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E int c ; /* a column index */ int length ; /* length of a row or column */ -#ifndef COLAMD_NDEBUG - int debug_rows ; - COLAMD_DEBUG2 (("Defrag..\n")) ; - for (psrc = &A[0] ; psrc < pfree ; psrc++) COLAMD_ASSERT (*psrc >= 0) ; - debug_rows = 0 ; -#endif /* COLAMD_NDEBUG */ - /* === Defragment the columns =========================================== */ pdest = &A[0] ; for (c = 0 ; c < n_col ; c++) { - if (EIGEN_COL_IS_ALIVE (c)) + if (COL_IS_ALIVE (c)) { psrc = &A [Col [c].start] ; @@ -1963,7 +1778,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (j = 0 ; j < length ; j++) { r = *psrc++ ; - if (EIGEN_ROW_IS_ALIVE (r)) + if (ROW_IS_ALIVE (r)) { *pdest++ = r ; } @@ -1976,26 +1791,22 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (r = 0 ; r < n_row ; r++) { - if (EIGEN_ROW_IS_ALIVE (r)) + if (ROW_IS_ALIVE (r)) { if (Row [r].length == 0) { /* this row is of zero length. cannot compact it, so kill it */ COLAMD_DEBUG3 (("Defrag row kill\n")) ; - EIGEN_KILL_ROW (r) ; + KILL_ROW (r) ; } else { /* save first column index in Row [r].shared2.first_column */ psrc = &A [Row [r].start] ; Row [r].shared2.first_column = *psrc ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (r)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; /* flag the start of the row with the one's complement of row */ - *psrc = EIGEN_ONES_COMPLEMENT (r) ; - -#ifndef COLAMD_NDEBUG - debug_rows++ ; -#endif /* COLAMD_NDEBUG */ + *psrc = ONES_COMPLEMENT (r) ; } } @@ -2011,11 +1822,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E { psrc-- ; /* get the row index */ - r = EIGEN_ONES_COMPLEMENT (*psrc) ; + r = ONES_COMPLEMENT (*psrc) ; COLAMD_ASSERT (r >= 0 && r < n_row) ; /* restore first column index */ *psrc = Row [r].shared2.first_column ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (r)) ; + COLAMD_ASSERT (ROW_IS_ALIVE (r)) ; /* move and compact the row */ COLAMD_ASSERT (pdest <= psrc) ; @@ -2024,17 +1835,13 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (j = 0 ; j < length ; j++) { c = *psrc++ ; - if (EIGEN_COL_IS_ALIVE (c)) + if (COL_IS_ALIVE (c)) { *pdest++ = c ; } } Row [r].length = (int) (pdest - &A [Row [r].start]) ; -#ifndef COLAMD_NDEBUG - debug_rows-- ; -#endif /* COLAMD_NDEBUG */ - } } /* ensure we found all the rows */ @@ -2047,7 +1854,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_clear_mark =========================================================== */ +/* === clear_mark =========================================================== */ /* ========================================================================== */ /* @@ -2055,12 +1862,12 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E Return value is the new tag_mark. Not user-callable. */ - int eigen_clear_mark /* return the new value for tag_mark */ +static inline int clear_mark /* return the new value for tag_mark */ ( /* === Parameters ======================================================= */ int n_row, /* number of rows in A */ - EIGEN_Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ + Colamd_Row Row [] /* Row [0 ... n_row-1].shared2.mark is set to zero */ ) { /* === Local variables ================================================== */ @@ -2069,7 +1876,7 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E for (r = 0 ; r < n_row ; r++) { - if (EIGEN_ROW_IS_ALIVE (r)) + if (ROW_IS_ALIVE (r)) { Row [r].shared2.mark = 0 ; } @@ -2080,13 +1887,13 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* ========================================================================== */ -/* === eigen_print_report ========================================================= */ +/* === print_report ========================================================= */ /* ========================================================================== */ - void eigen_print_report +static void print_report ( - char *method, - int stats [EIGEN_COLAMD_STATS] + const char *method, + int stats [COLAMD_STATS] ) { @@ -2098,11 +1905,11 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E return ; } - i1 = stats [EIGEN_COLAMD_INFO1] ; - i2 = stats [EIGEN_COLAMD_INFO2] ; - i3 = stats [EIGEN_COLAMD_INFO3] ; + i1 = stats [COLAMD_INFO1] ; + i2 = stats [COLAMD_INFO2] ; + i3 = stats [COLAMD_INFO3] ; - if (stats [EIGEN_COLAMD_STATUS] >= 0) + if (stats [COLAMD_STATUS] >= 0) { PRINTF ("%s: OK. ", method) ; } @@ -2111,10 +1918,10 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E PRINTF ("%s: ERROR. ", method) ; } - switch (stats [EIGEN_COLAMD_STATUS]) + switch (stats [COLAMD_STATUS]) { - case EIGEN_COLAMD_OK_BUT_JUMBLED: + case COLAMD_OK_BUT_JUMBLED: PRINTF ("Matrix has unsorted or duplicate row indices.\n") ; @@ -2129,77 +1936,77 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E /* no break - fall through to next case instead */ - case EIGEN_COLAMD_OK: + case COLAMD_OK: PRINTF ("\n") ; PRINTF ("%s: number of dense or empty rows ignored: %d\n", - method, stats [EIGEN_COLAMD_DENSE_ROW]) ; + method, stats [COLAMD_DENSE_ROW]) ; PRINTF ("%s: number of dense or empty columns ignored: %d\n", - method, stats [EIGEN_COLAMD_DENSE_COL]) ; + method, stats [COLAMD_DENSE_COL]) ; PRINTF ("%s: number of garbage collections performed: %d\n", - method, stats [EIGEN_COLAMD_DEFRAG_COUNT]) ; + method, stats [COLAMD_DEFRAG_COUNT]) ; break ; - case EIGEN_COLAMD_ERROR_A_not_present: + case COLAMD_ERROR_A_not_present: PRINTF ("Array A (row indices of matrix) not present.\n") ; break ; - case EIGEN_COLAMD_ERROR_p_not_present: + case COLAMD_ERROR_p_not_present: PRINTF ("Array p (column pointers for matrix) not present.\n") ; break ; - case EIGEN_COLAMD_ERROR_nrow_negative: + case COLAMD_ERROR_nrow_negative: PRINTF ("Invalid number of rows (%d).\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_ncol_negative: + case COLAMD_ERROR_ncol_negative: PRINTF ("Invalid number of columns (%d).\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_nnz_negative: + case COLAMD_ERROR_nnz_negative: PRINTF ("Invalid number of nonzero entries (%d).\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_p0_nonzero: + case COLAMD_ERROR_p0_nonzero: PRINTF ("Invalid column pointer, p [0] = %d, must be zero.\n", i1) ; break ; - case EIGEN_COLAMD_ERROR_A_too_small: + case COLAMD_ERROR_A_too_small: PRINTF ("Array A too small.\n") ; PRINTF (" Need Alen >= %d, but given only Alen = %d.\n", i1, i2) ; break ; - case EIGEN_COLAMD_ERROR_col_length_negative: + case COLAMD_ERROR_col_length_negative: PRINTF ("Column %d has a negative number of nonzero entries (%d).\n", INDEX (i1), i2) ; break ; - case EIGEN_COLAMD_ERROR_row_index_out_of_bounds: + case COLAMD_ERROR_row_index_out_of_bounds: PRINTF ("Row index (row %d) out of bounds (%d to %d) in column %d.\n", INDEX (i2), INDEX (0), INDEX (i3-1), INDEX (i1)) ; break ; - case EIGEN_COLAMD_ERROR_out_of_memory: + case COLAMD_ERROR_out_of_memory: PRINTF ("Out of memory.\n") ; break ; - case EIGEN_COLAMD_ERROR_internal_error: + case COLAMD_ERROR_internal_error: /* if this happens, there is a bug in the code */ PRINTF @@ -2208,307 +2015,5 @@ bool eigen_colamd(int n_row, int n_col, int Alen, int *A, int *p, double knobs[E } } - - - -/* ========================================================================== */ -/* === eigen_colamd debugging routines ============================================ */ -/* ========================================================================== */ - -/* When debugging is disabled, the remainder of this file is ignored. */ - -#ifndef COLAMD_NDEBUG - - -/* ========================================================================== */ -/* === eigen_debug_structures ===================================================== */ -/* ========================================================================== */ - -/* - At this point, all empty rows and columns are dead. All live columns - are "clean" (containing no dead rows) and simplicial (no supercolumns - yet). Rows may contain dead columns, but all live rows contain at - least one live column. -*/ - - void eigen_debug_structures -( - /* === Parameters ======================================================= */ - - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [], - int n_col2 -) -{ - /* === Local variables ================================================== */ - - int i ; - int c ; - int *cp ; - int *cp_end ; - int len ; - int score ; - int r ; - int *rp ; - int *rp_end ; - int deg ; - - /* === Check A, Row, and Col ============================================ */ - - for (c = 0 ; c < n_col ; c++) - { - if (EIGEN_COL_IS_ALIVE (c)) - { - len = Col [c].length ; - score = Col [c].shared2.score ; - COLAMD_DEBUG4 (("initial live col %5d %5d %5d\n", c, len, score)) ; - COLAMD_ASSERT (len > 0) ; - COLAMD_ASSERT (score >= 0) ; - COLAMD_ASSERT (Col [c].shared1.thickness == 1) ; - cp = &A [Col [c].start] ; - cp_end = cp + len ; - while (cp < cp_end) - { - r = *cp++ ; - COLAMD_ASSERT (EIGEN_ROW_IS_ALIVE (r)) ; - } - } - else - { - i = Col [c].shared2.order ; - COLAMD_ASSERT (i >= n_col2 && i < n_col) ; - } - } - - for (r = 0 ; r < n_row ; r++) - { - if (EIGEN_ROW_IS_ALIVE (r)) - { - i = 0 ; - len = Row [r].length ; - deg = Row [r].shared1.degree ; - COLAMD_ASSERT (len > 0) ; - COLAMD_ASSERT (deg > 0) ; - rp = &A [Row [r].start] ; - rp_end = rp + len ; - while (rp < rp_end) - { - c = *rp++ ; - if (EIGEN_COL_IS_ALIVE (c)) - { - i++ ; - } - } - COLAMD_ASSERT (i > 0) ; - } - } -} - - -/* ========================================================================== */ -/* === eigen_debug_deg_lists ====================================================== */ -/* ========================================================================== */ - -/* - Prints the contents of the degree lists. Counts the number of columns - in the degree list and compares it to the total it should have. Also - checks the row degrees. -*/ - - void eigen_debug_deg_lists -( - /* === Parameters ======================================================= */ - - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int head [], - int min_score, - int should, - int max_deg -) -{ - /* === Local variables ================================================== */ - - int deg ; - int col ; - int have ; - int row ; - - /* === Check the degree lists =========================================== */ - - if (n_col > 10000 && colamd_debug <= 0) - { - return ; - } - have = 0 ; - COLAMD_DEBUG4 (("Degree lists: %d\n", min_score)) ; - for (deg = 0 ; deg <= n_col ; deg++) - { - col = head [deg] ; - if (col == EIGEN_COLAMD_EMPTY) - { - continue ; - } - COLAMD_DEBUG4 (("%d:", deg)) ; - while (col != EIGEN_COLAMD_EMPTY) - { - COLAMD_DEBUG4 ((" %d", col)) ; - have += Col [col].shared1.thickness ; - COLAMD_ASSERT (EIGEN_COL_IS_ALIVE (col)) ; - col = Col [col].shared4.degree_next ; - } - COLAMD_DEBUG4 (("\n")) ; - } - COLAMD_DEBUG4 (("should %d have %d\n", should, have)) ; - COLAMD_ASSERT (should == have) ; - - /* === Check the row degrees ============================================ */ - - if (n_row > 10000 && colamd_debug <= 0) - { - return ; - } - for (row = 0 ; row < n_row ; row++) - { - if (EIGEN_ROW_IS_ALIVE (row)) - { - COLAMD_ASSERT (Row [row].shared1.degree <= max_deg) ; - } - } -} - - -/* ========================================================================== */ -/* === eigen_debug_mark =========================================================== */ -/* ========================================================================== */ - -/* - Ensures that the tag_mark is less that the maximum and also ensures that - each entry in the mark array is less than the tag mark. -*/ - - void eigen_debug_mark -( - /* === Parameters ======================================================= */ - - int n_row, - EIGEN_Colamd_Row Row [], - int tag_mark, - int max_mark -) -{ - /* === Local variables ================================================== */ - - int r ; - - /* === Check the Row marks ============================================== */ - - COLAMD_ASSERT (tag_mark > 0 && tag_mark <= max_mark) ; - if (n_row > 10000 && colamd_debug <= 0) - { - return ; - } - for (r = 0 ; r < n_row ; r++) - { - COLAMD_ASSERT (Row [r].shared2.mark < tag_mark) ; - } -} - - -/* ========================================================================== */ -/* === eigen_debug_matrix ========================================================= */ -/* ========================================================================== */ - -/* - Prints out the contents of the columns and the rows. -*/ - - void eigen_debug_matrix -( - /* === Parameters ======================================================= */ - - int n_row, - int n_col, - EIGEN_Colamd_Row Row [], - EIGEN_Colamd_Col Col [], - int A [] -) -{ - /* === Local variables ================================================== */ - - int r ; - int c ; - int *rp ; - int *rp_end ; - int *cp ; - int *cp_end ; - - /* === Dump the rows and columns of the matrix ========================== */ - - if (colamd_debug < 3) - { - return ; - } - COLAMD_DEBUG3 (("DUMP MATRIX:\n")) ; - for (r = 0 ; r < n_row ; r++) - { - COLAMD_DEBUG3 (("Row %d alive? %d\n", r, EIGEN_ROW_IS_ALIVE (r))) ; - if (EIGEN_ROW_IS_DEAD (r)) - { - continue ; - } - COLAMD_DEBUG3 (("start %d length %d degree %d\n", - Row [r].start, Row [r].length, Row [r].shared1.degree)) ; - rp = &A [Row [r].start] ; - rp_end = rp + Row [r].length ; - while (rp < rp_end) - { - c = *rp++ ; - COLAMD_DEBUG4 ((" %d col %d\n", EIGEN_COL_IS_ALIVE (c), c)) ; - } - } - - for (c = 0 ; c < n_col ; c++) - { - COLAMD_DEBUG3 (("Col %d alive? %d\n", c, EIGEN_COL_IS_ALIVE (c))) ; - if (EIGEN_COL_IS_DEAD (c)) - { - continue ; - } - COLAMD_DEBUG3 (("start %d length %d shared1 %d shared2 %d\n", - Col [c].start, Col [c].length, - Col [c].shared1.thickness, Col [c].shared2.score)) ; - cp = &A [Col [c].start] ; - cp_end = cp + Col [c].length ; - while (cp < cp_end) - { - r = *cp++ ; - COLAMD_DEBUG4 ((" %d row %d\n", EIGEN_ROW_IS_ALIVE (r), r)) ; - } - } -} - - void eigen_colamd_get_debug -( - char *method -) -{ - colamd_debug = 0 ; /* no debug printing */ - - /* get "D" environment variable, which gives the debug printing level */ - if (getenv ("D")) - { - colamd_debug = atoi (getenv ("D")) ; - } - - COLAMD_DEBUG0 (("%s: debug version, D = %d (THIS WILL BE SLOW!)\n", - method, colamd_debug)) ; -} - -#endif /* NDEBUG */ +} // namespace internal #endif diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h index 47cd6f169..f5757b319 100644 --- a/Eigen/src/OrderingMethods/Ordering.h +++ b/Eigen/src/OrderingMethods/Ordering.h @@ -27,8 +27,10 @@ #define EIGEN_ORDERING_H #include "Amd.h" -#include "Eigen_Colamd.h" namespace Eigen { + +#include "Eigen_Colamd.h" + namespace internal { /** @@ -131,18 +133,18 @@ class COLAMDOrdering int n = mat.cols(); int nnz = mat.nonZeros(); // Get the recommended value of Alen to be used by colamd - int Alen = eigen_colamd_recommended(nnz, m, n); + int Alen = internal::colamd_recommended(nnz, m, n); // Set the default parameters - double knobs [EIGEN_COLAMD_KNOBS]; - int stats [EIGEN_COLAMD_STATS]; - eigen_colamd_set_defaults(knobs); + double knobs [COLAMD_KNOBS]; + int stats [COLAMD_STATS]; + internal::colamd_set_defaults(knobs); int info; IndexVector p(n+1), A(Alen); for(int i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; // Call Colamd routine to compute the ordering - info = eigen_colamd(m, n, Alen, A.data(), p.data(), knobs, stats); + info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); eigen_assert( info && "COLAMD failed " ); perm.resize(n); diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index e2076138a..77df091c3 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -205,7 +205,7 @@ class SparseLU void initperfvalues() { m_perfv.panel_size = 12; - m_perfv.relax = 6; + m_perfv.relax = 1; m_perfv.maxsuper = 100; m_perfv.rowblk = 200; m_perfv.colblk = 60; diff --git a/bench/spbench/CMakeLists.txt b/bench/spbench/CMakeLists.txt index 5451843b9..6e0e1b103 100644 --- a/bench/spbench/CMakeLists.txt +++ b/bench/spbench/CMakeLists.txt @@ -55,6 +55,12 @@ if(PASTIX_FOUND AND BLAS_FOUND) set(PASTIX_ALL_LIBS ${PASTIX_LIBRARIES} ${BLAS_LIBRARIES}) endif(PASTIX_FOUND AND BLAS_FOUND) +if(METIS_FOUND) + include_directories(${METIS_INCLUDES}) + set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) + add_definitions("-DEIGEN_METIS_SUPPORT") +endif(METIS_FOUND) + find_library(RT_LIBRARY rt) if(RT_LIBRARY) set(SPARSE_LIBS ${SPARSE_LIBS} ${RT_LIBRARY}) @@ -66,11 +72,6 @@ target_link_libraries (spbenchsolver ${SPARSE_LIBS}) add_executable(spsolver sp_solver.cpp) target_link_libraries (spsolver ${SPARSE_LIBS}) -if(METIS_FOUND) - include_directories(${METIS_INCLUDES}) - set (SPARSE_LIBS ${SPARSE_LIBS} ${METIS_LIBRARIES}) - add_definitions("-DEIGEN_METIS_SUPPORT") -endif(METIS_FOUND) add_executable(test_sparseLU test_sparseLU.cpp) target_link_libraries (test_sparseLU ${SPARSE_LIBS}) diff --git a/bench/spbench/spbenchsolver.h b/bench/spbench/spbenchsolver.h index c48ed7aa7..19c719c04 100644 --- a/bench/spbench/spbenchsolver.h +++ b/bench/spbench/spbenchsolver.h @@ -21,9 +21,14 @@ #include #include #include +#include #include "spbenchstyle.h" +#ifdef EIGEN_METIS_SUPPORT +#include +#endif + #ifdef EIGEN_CHOLMOD_SUPPORT #include #endif @@ -45,26 +50,27 @@ #endif // CONSTANTS -#define EIGEN_UMFPACK 0 -#define EIGEN_SUPERLU 1 -#define EIGEN_PASTIX 2 -#define EIGEN_PARDISO 3 -#define EIGEN_BICGSTAB 4 -#define EIGEN_BICGSTAB_ILUT 5 -#define EIGEN_GMRES 6 -#define EIGEN_GMRES_ILUT 7 -#define EIGEN_SIMPLICIAL_LDLT 8 -#define EIGEN_CHOLMOD_LDLT 9 -#define EIGEN_PASTIX_LDLT 10 -#define EIGEN_PARDISO_LDLT 11 -#define EIGEN_SIMPLICIAL_LLT 12 -#define EIGEN_CHOLMOD_SUPERNODAL_LLT 13 -#define EIGEN_CHOLMOD_SIMPLICIAL_LLT 14 -#define EIGEN_PASTIX_LLT 15 -#define EIGEN_PARDISO_LLT 16 -#define EIGEN_CG 17 -#define EIGEN_CG_PRECOND 18 -#define EIGEN_ALL_SOLVERS 19 +#define EIGEN_UMFPACK 10 +#define EIGEN_SUPERLU 20 +#define EIGEN_PASTIX 30 +#define EIGEN_PARDISO 40 +#define EIGEN_SPARSELU_COLAMD 50 +#define EIGEN_SPARSELU_METIS 51 +#define EIGEN_BICGSTAB 60 +#define EIGEN_BICGSTAB_ILUT 61 +#define EIGEN_GMRES 70 +#define EIGEN_GMRES_ILUT 71 +#define EIGEN_SIMPLICIAL_LDLT 80 +#define EIGEN_CHOLMOD_LDLT 90 +#define EIGEN_PASTIX_LDLT 100 +#define EIGEN_PARDISO_LDLT 110 +#define EIGEN_SIMPLICIAL_LLT 120 +#define EIGEN_CHOLMOD_SUPERNODAL_LLT 130 +#define EIGEN_CHOLMOD_SIMPLICIAL_LLT 140 +#define EIGEN_PASTIX_LLT 150 +#define EIGEN_PARDISO_LLT 160 +#define EIGEN_CG 170 +#define EIGEN_CG_PRECOND 180 using namespace Eigen; using namespace std; @@ -188,6 +194,17 @@ void printStatheader(std::ofstream& out) out << " EIGEN \n"; out << " \n"; + out <<" \n"; + out << " LU_COLAMD \n"; + out << " EIGEN \n"; + out << " \n"; + +#ifdef EIGEN_METIS_SUPPORT + out <<" \n"; + out << " LU_METIS \n"; + out << " EIGEN \n"; + out << " \n"; +#endif out << " \n"; } @@ -325,8 +342,19 @@ void SelectSolvers(const SparseMatrix&A, unsigned int sym, Matrix > solver; + call_directsolver(solver, EIGEN_SPARSELU_COLAMD, A, b, refX, statFile); + // Eigen SparseLU METIS + #ifdef EIGEN_METIS_SUPPORT + { + cout << "\n Solving with Sparse LU AND METIS ... \n"; + SparseLU > solver; + call_directsolver(solver, EIGEN_SPARSELU_METIS, A, b, refX, statFile); + } + #endif //BiCGSTAB { -- cgit v1.2.3