Refactoring codes for numeric updates
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 25fad0f..474dfde 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h
@@ -388,6 +388,7 @@ #include "SparseLU_snode_bmod.h" #include "SparseLU_pivotL.h" #include "SparseLU_panel_dfs.h" +#include "SparseLU_kernel_bmod.h" #include "SparseLU_panel_bmod.h" #include "SparseLU_column_dfs.h" #include "SparseLU_column_bmod.h"
diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h index 0078772..1457b6f 100644 --- a/Eigen/src/SparseLU/SparseLU_column_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h
@@ -66,7 +66,7 @@ typedef typename IndexVector::Scalar Index; typedef typename ScalarVector::Scalar Scalar; int jsupno, k, ksub, krep, ksupno; - int lptr, nrow, isub, i, irow, nextlu, new_next, ufirst; + int lptr, nrow, isub, irow, nextlu, new_next, ufirst; int fsupc, nsupc, nsupr, luptr, kfnz, no_zeros; /* krep = representative of current k-th supernode * fsupc = first supernodal column @@ -122,46 +122,7 @@ // Perform a triangular solver and block update, // then scatter the result of sup-col update to dense no_zeros = kfnz - fst_col; - // First, copy U[*,j] segment from dense(*) to tempv(*) - isub = lptr + no_zeros; - for (i = 0; i < segsize; i++) - { - irow = lsub(isub); - tempv(i) = dense(irow); - ++isub; - } - // Dense triangular solve -- start effective triangle - luptr += nsupr * no_zeros + no_zeros; - // Form Eigen matrix and vector - Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); - VectorBlock<ScalarVector> u(tempv, 0, segsize); - - u = A.template triangularView<UnitLower>().solve(u); - - // Dense matrix-vector product y <-- A*x - luptr += segsize; - new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); - VectorBlock<ScalarVector> l(tempv, segsize, nrow); - l= A * u; - - // Scatter tempv[] into SPA dense[] as a temporary storage - isub = lptr + no_zeros; - for (i = 0; i < segsize; i++) - { - irow = lsub(isub); - dense(irow) = tempv(i); - tempv(i) = Scalar(0.0); - ++isub; - } - - // Scatter l into SPA dense[] - for (i = 0; i < nrow; i++) - { - irow = lsub(isub); - dense(irow) -= l(i); - l(i) = Scalar(0.0); - ++isub; - } + LU_kernel_bmod(segsize, dense, tempv, lusup, luptr, nsupr, nrow, lsub, lptr, no_zeros); } // end if jsupno } // end for each segment
diff --git a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h new file mode 100644 index 0000000..d5df70f --- /dev/null +++ b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h
@@ -0,0 +1,92 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef SPARSELU_KERNEL_BMOD_H +#define SPARSELU_KERNEL_BMOD_H + +/** + * \brief Performs numeric block updates from a given supernode to a single column + * + * \param segsize Size of the segment (and blocks ) to use for updates + * \param [in,out]dense Packed values of the original matrix + * \param tempv temporary vector to use for updates + * \param lusup array containing the supernodes + * \param nsupr Number of rows in the supernode + * \param nrow Number of rows in the rectangular part of the supernode + * \param lsub compressed row subscripts of supernodes + * \param lptr pointer to the first column of the current supernode in lsub + * \param no_zeros Number of nonzeros elements before the diagonal part of the supernode + * \return 0 on success + */ +template <typename BlockScalarVector, typename ScalarVector, typename IndexVector> +int LU_kernel_bmod(const int segsize, BlockScalarVector& dense, ScalarVector& tempv, ScalarVector& lusup, int& luptr, const int nsupr, const int nrow, IndexVector& lsub, const int lptr, const int no_zeros) +{ + typedef typename ScalarVector::Scalar Scalar; + // First, copy U[*,j] segment from dense(*) to tempv(*) + // The result of triangular solve is in tempv[*]; + // The result of matric-vector update is in dense[*] + int isub = lptr + no_zeros; + int i, irow; + for (i = 0; i < segsize; i++) + { + irow = lsub(isub); + tempv(i) = dense(irow); + ++isub; + } + // Dense triangular solve -- start effective triangle + luptr += nsupr * no_zeros + no_zeros; + // Form Eigen matrix and vector + Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); + VectorBlock<ScalarVector> u(tempv, 0, segsize); + + u = A.template triangularView<UnitLower>().solve(u); + + // Dense matrix-vector product y <-- A*x + luptr += segsize; + new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); + VectorBlock<ScalarVector> l(tempv, segsize, nrow); + l= A * u; + + // Scatter tempv[] into SPA dense[] as a temporary storage + isub = lptr + no_zeros; + for (i = 0; i < segsize; i++) + { + irow = lsub(isub); + dense(irow) = tempv(i); + tempv(i) = Scalar(0.0); + ++isub; + } + + // Scatter l into SPA dense[] + for (i = 0; i < nrow; i++) + { + irow = lsub(isub); + dense(irow) -= l(i); + l(i) = Scalar(0.0); + ++isub; + } + + return 0; +} +#endif \ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h index 59ec69e..ebff787 100644 --- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h +++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h
@@ -73,12 +73,12 @@ IndexVector& xlusup = glu.xlusup; ScalarVector& lusup = glu.lusup; - int i,ksub,jj,nextl_col,irow; + int ksub,jj,nextl_col; int fsupc, nsupc, nsupr, nrow; int krep, kfnz; int lptr; // points to the row subscripts of a supernode int luptr; // ... - int segsize,no_zeros,isub ; + int segsize,no_zeros ; // For each nonz supernode segment of U[*,j] in topological order int k = nseg - 1; for (ksub = 0; ksub < nseg; ksub++) @@ -118,52 +118,7 @@ // Perform a trianglar solve and block update, // then scatter the result of sup-col update to dense[] no_zeros = kfnz - fsupc; - // First Copy U[*,j] segment from dense[*] to tempv[*] : - // The result of triangular solve is in tempv[*]; - // The result of matric-vector update is in dense_col[*] - isub = lptr + no_zeros; - for (i = 0; i < segsize; ++i) - { - irow = lsub(isub); - tempv(i) = dense_col(irow); // Gather to a compact vector - ++isub; - } - // Start effective triangle - luptr += nsupr * no_zeros + no_zeros; - // triangular solve with Eigen - Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); - VectorBlock<ScalarVector> u(tempv, 0, segsize); - u = A.template triangularView<UnitLower>().solve(u); - - luptr += segsize; - // Dense Matrix vector product y <-- A*x; - new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); - VectorBlock<ScalarVector> l(tempv, segsize, nrow); - l= A * u; - - // Scatter tempv(*) into SPA dense(*) such that tempv(*) - // can be used for the triangular solve of the next - // column of the panel. The y will be copied into ucol(*) - // after the whole panel has been finished... after column_dfs() and column_bmod() - - isub = lptr + no_zeros; - for (i = 0; i < segsize; i++) - { - irow = lsub(isub); - dense_col(irow) = tempv(i); - tempv(i) = Scalar(0.0); - isub++; - } - - // Scatter the update from &tempv[segsize] into SPA dense(*) - // Start dense rectangular L - for (i = 0; i < nrow; i++) - { - irow = lsub(isub); - dense_col(irow) -= l(i); - l(i) = Scalar(0); - ++isub; - } + LU_kernel_bmod(segsize, dense_col, tempv, lusup, luptr, nsupr, nrow, lsub, lptr, no_zeros); } // End for each column in the panel } // End for each updating supernode
diff --git a/bench/spbench/sp_solver.cpp b/bench/spbench/sp_solver.cpp new file mode 100644 index 0000000..e18f2d1 --- /dev/null +++ b/bench/spbench/sp_solver.cpp
@@ -0,0 +1,124 @@ +// Small bench routine for Eigen available in Eigen +// (C) Desire NUENTSA WAKAM, INRIA + +#include <iostream> +#include <fstream> +#include <iomanip> +#include <Eigen/Jacobi> +#include <Eigen/Householder> +#include <Eigen/IterativeLinearSolvers> +#include <Eigen/LU> +#include <unsupported/Eigen/SparseExtra> +//#include <Eigen/SparseLU> +#include <Eigen/SuperLUSupport> +// #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h> +#include <bench/BenchTimer.h> + +using namespace std; +using namespace Eigen; + +int main(int argc, char **args) +{ + SparseMatrix<double, ColMajor> A; + typedef SparseMatrix<double, ColMajor>::Index Index; + typedef Matrix<double, Dynamic, Dynamic> DenseMatrix; + typedef Matrix<double, Dynamic, 1> DenseRhs; + VectorXd b, x, tmp; + BenchTimer timer,totaltime; + //SparseLU<SparseMatrix<double, ColMajor> > solver; + SuperLU<SparseMatrix<double, ColMajor> > solver; + ifstream matrix_file; + string line; + int n; + // Set parameters +// solver.iparm(IPARM_THREAD_NBR) = 4; + /* Fill the matrix with sparse matrix stored in Matrix-Market coordinate column-oriented format */ + if (argc < 2) assert(false && "please, give the matrix market file "); + + timer.start(); + totaltime.start(); + loadMarket(A, args[1]); + cout << "End charging matrix " << endl; + bool iscomplex=false, isvector=false; + int sym; + getMarketHeader(args[1], sym, iscomplex, isvector); + 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<double, ColMajor> temp; + temp = A; + A = temp.selfadjointView<Lower>(); + } + timer.stop(); + + n = A.cols(); + // ====== TESTS FOR SPARSE TUTORIAL ====== +// cout<< "OuterSize " << A.outerSize() << " inner " << A.innerSize() << endl; +// SparseMatrix<double, RowMajor> mat1(A); +// SparseMatrix<double, RowMajor> mat2; +// cout << " norm of A " << mat1.norm() << endl; ; +// PermutationMatrix<Dynamic, Dynamic, int> perm(n); +// perm.resize(n,1); +// perm.indices().setLinSpaced(n, 0, n-1); +// mat2 = perm * mat1; +// mat.subrows(); +// mat2.resize(n,n); +// mat2.reserve(10); +// mat2.setConstant(); +// std::cout<< "NORM " << mat1.squaredNorm()<< endl; + + cout<< "Time to load the matrix " << timer.value() <<endl; + /* Fill the right hand side */ + +// solver.set_restart(374); + if (argc > 2) + loadMarketVector(b, args[2]); + else + { + b.resize(n); + tmp.resize(n); +// tmp.setRandom(); + for (int i = 0; i < n; i++) tmp(i) = i; + b = A * tmp ; + } +// Scaling<SparseMatrix<double> > scal; +// scal.computeRef(A); +// b = scal.LeftScaling().cwiseProduct(b); + + /* Compute the factorization */ + cout<< "Starting the factorization "<< endl; + timer.reset(); + timer.start(); + cout<< "Size of Input Matrix "<< b.size()<<"\n\n"; + cout<< "Rows and columns "<< A.rows() <<" " <<A.cols() <<"\n"; + solver.compute(A); +// solver.analyzePattern(A); +// solver.factorize(A); + if (solver.info() != Success) { + std::cout<< "The solver failed \n"; + return -1; + } + timer.stop(); + float time_comp = timer.value(); + cout <<" Compute Time " << time_comp<< endl; + + timer.reset(); + timer.start(); + x = solver.solve(b); +// x = scal.RightScaling().cwiseProduct(x); + timer.stop(); + float time_solve = timer.value(); + cout<< " Time to solve " << time_solve << endl; + + /* Check the accuracy */ + VectorXd tmp2 = b - A*x; + double tempNorm = tmp2.norm()/b.norm(); + cout << "Relative norm of the computed solution : " << tempNorm <<"\n"; +// cout << "Iterations : " << solver.iterations() << "\n"; + + totaltime.stop(); + cout << "Total time " << totaltime.value() << "\n"; +// std::cout<<x.transpose()<<"\n"; + + return 0; +} \ No newline at end of file
diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp index 08b6c92..ecf254b 100644 --- a/bench/spbench/test_sparseLU.cpp +++ b/bench/spbench/test_sparseLU.cpp
@@ -13,8 +13,8 @@ int main(int argc, char **args) { - typedef complex<double> scalar; -// typedef double scalar; +// typedef complex<double> scalar; + typedef double scalar; SparseMatrix<scalar, ColMajor> A; typedef SparseMatrix<scalar, ColMajor>::Index Index; typedef Matrix<scalar, Dynamic, Dynamic> DenseMatrix;