Reply to topic

Use Eigen together with KLU

mhuck
Registered Member
Posts
5
Karma
0

Use Eigen together with KLU

Mon Jul 08, 2013 12:08 pm
Hello,

in my simulation I need a large (and sparse) electrical equivalent circuit. At the Moment I am using the Eigen library for Managing the matrixes and solving the equations Systems.
However I found that for the case of a unsymmetrical Matrix KLU is much faster than SuperLU. Since Eigen works well I want to keep using it.
Can I use the Eigen::SparseMatrx class with KLU?

My guess is yes, since you can bring the SparseMatrix Class in the compressed-column form. The question is what are the needed Arrays for KLU (I refer to the KLU User Guide chapter 5.4 ):
  1. n: size of the Matrix this one is clear
  2. Ap this may be OuterStarts
  3. Ai this may be InnerIndices
  4. Ax this may be Values
For the Ap,Ai,Ax I am not sure which SparseMatrix methods I need to use. I think it should be outerIndexPtr, innerIndexPtr and valuePtr, but I am really unsure.
I wont answer until 21.07, but be free to discuss
User avatar ggael
Moderator
Posts
2194
Karma
15
OS

Re: Use Eigen together with KLU

Mon Jul 08, 2013 6:29 pm
You guessed right, and you can find similar example in the code of the CholmodSupport and SPQRSupport modules. I'd also be curious to know how does BICGSTAB<MatrixXd,IncompleteLUT<double> > would perform on such problems.
mhuck
Registered Member
Posts
5
Karma
0

Re: Use Eigen together with KLU

Tue Jul 09, 2013 6:50 am
Ok,

thanks. I will try that when I am back.
mhuck
Registered Member
Posts
5
Karma
0

Re: Use Eigen together with KLU

Wed Jul 24, 2013 1:29 pm
I'm running in some trouble with KLU.
I try to do the sample in the KLU Userguide p18 with a Eigen Matrix (the sample itself with C Arrays works fine)
my code:
Code: Select all
int main (void)
{
klu_symbolic *Symbolic ;
klu_numeric *Numeric ;
klu_common Common ;
int i ;
klu_defaults (&Common) ;
test.resize(5,5);
test.reserve(12);
test.insert(0,0)=2;
test.insert(0,1)=3;
test.insert(1,0)=3;
test.insert(1,2)=4;
test.insert(1,4)=6;
test.insert(2,1)=-1;
test.insert(2,2)=-3;
test.insert(2,3)=2;
test.insert(3,2)=1;
test.insert(4,1)=4;
test.insert(4,2)=2;
test.insert(4,4)=1;
test.finalize();
test.makeCompressed();

std::cout<<test;
Symbolic = klu_analyze (test.size(), test.outerIndexPtr(), test.innerIndexPtr(), &Common) ;
Numeric = klu_factor (test.outerIndexPtr(), test.innerIndexPtr(), test.valuePtr(), Symbolic, &Common) ;
klu_solve (Symbolic, Numeric, 5, 1, b, &Common) ;
klu_free_symbolic (&Symbolic, &Common) ;
klu_free_numeric (&Numeric, &Common) ;
for (i = 0 ; i < n ; i++) printf ("x [%d] = %g\n", i, b [i]) ;
return (0) ;
}

where test is a Eigen::SparseMatrix<double>.

The Output ist:
Code: Select all
Nonzero entries:
(2,0) (3,1) (3,0) (-1,2) (4,4) (4,1) (-3,2) (1,3) (2,4) (2,2) (6,1) (1,4)

Outer pointers:
0 2 5 9 10  $

2 3 0 0 0
3 0 4 0 6
0 -1 -3 2 0
0 0 1 0 0
0 4 2 0 1
 


Then it crashes in the klu_analyze function.
What I have noticed ist that in the KLU sample there is an additional outerIndex : 12 , is there a reason why it isnt there in Eigen ?
Dee33
Registered Member
Posts
54
Karma
0
OS

Re: Use Eigen together with KLU

Wed Jul 24, 2013 1:49 pm
Everything seems normal in your output. We don't just output the last term (which is the number of nonzeros). But you can find it in the array test.outerIndexPtr()
mhuck
Registered Member
Posts
5
Karma
0

Re: Use Eigen together with KLU  Topic is solved

Wed Jul 24, 2013 2:16 pm
I found my error: test.size is the number of elements but I need the number of columns /rows.

Sorry and thanks.
User avatar ggael
Moderator
Posts
2194
Karma
15
OS

Re: Use Eigen together with KLU

Wed Jul 24, 2013 2:45 pm
Feel free to share your KLU wrapper once it's ready! I'm sure it will be useful to many people.
mhuck
Registered Member
Posts
5
Karma
0

Re: Use Eigen together with KLU

Fri Aug 09, 2013 1:20 pm
Ok here is my wrapper. It is not sophisticated, because I have no time to test it with various templates.
At the current state it works only for x,b as vector (in Ax=b) but I think you can expand it to matrices if you Change the "1" in :
Code: Select all
m_ok = klu_solve(m_symbolic_klu, m_nummeric_klu, m_MatSize, 1, _RhS.data(), &m_common_klu);


I hope it helps some one.
Code: Select all
/*
 * KLUWrapper.hpp
 *
 *  Edited on: Aug 9, 2013
 *      Author: Moritz Huck
 */

#ifndef KLUWRAPPER_HPP_
#define KLUWRAPPER_HPP_

#include  "Eigen/src/Core/util/Meta.h"
#include "klu.h"

namespace Eigen {
template<typename _valuetype>
class KLU_Wrapper: public Eigen::internal::noncopyable {
public:
   typedef Eigen::SparseMatrix<_valuetype> _SparseMatrix;




   void analyzePattern( _SparseMatrix Matrix){
      eigen_assert(Matrix.isCompressed() && "call makeCompressed() on the Matrix first!");
      m_MatSize=Matrix.innerSize();
      printf("NodesSize Solver = %i",m_MatSize);
      analyze_symbolic(m_MatSize,Matrix);

      m_analised_Pattern=true;
   }

   void factorize(_SparseMatrix Matrix){
      eigen_assert(Matrix.isCompressed() && "call makeCompressed() on the Matrix first!");
      eigen_assert(m_analised_Pattern && "You must first call analyzePattern()");
      if(m_is_factored){refactor_nummeric(Matrix);}
      else{factor_nummeric(Matrix);}
      m_factorized=true;
   }


   Eigen::VectorXd solve(Eigen::VectorXd _RhS){
      eigen_assert(m_factorized && "You must first call analyzePattern() and factorize()" );
      m_ok = klu_solve(m_symbolic_klu, m_nummeric_klu, m_MatSize, 1, _RhS.data(), &m_common_klu);
            if (m_ok == 0) throw std::runtime_error("klu_solvefehlgeschlagen");
            return _RhS;
   }



   KLU_Wrapper():m_symbolic_is_done(false),m_is_factored(false),m_analised_Pattern(false),m_factorized(false),m_ok(1)
   {
      m_ok = klu_defaults(&m_common_klu);
      if (m_ok == 0) throw std::runtime_error("klu_defaults failed");
   }
   virtual ~KLU_Wrapper(){
      klu_free_symbolic(&m_symbolic_klu,&m_common_klu);
      klu_free_numeric(&m_nummeric_klu,&m_common_klu);
   }


protected:

            klu_numeric* m_nummeric_klu;
            klu_symbolic* m_symbolic_klu;
            klu_common m_common_klu;
            int m_ok;
            int m_MatSize;
            bool m_symbolic_is_done;
            bool m_is_factored;
            bool m_analised_Pattern;
            bool m_factorized;


            void factor_nummeric( _SparseMatrix Matrix){
                  m_nummeric_klu = klu_factor(Matrix.outerIndexPtr(), Matrix.innerIndexPtr(), Matrix.valuePtr(),
                           m_symbolic_klu, &m_common_klu);
                     if (m_nummeric_klu == NULL) {throw std::runtime_error("klu_factor failed");}
                     else {m_is_factored=true;} //::TODO:: write real exception
               }


            void analyze_symbolic(int MatSize, _SparseMatrix Matrix){
                  m_symbolic_klu = klu_analyze(MatSize, Matrix.outerIndexPtr(), Matrix.innerIndexPtr(), &m_common_klu);
                  if (m_symbolic_klu == NULL) {
                        std::cout << m_symbolic_klu << std::endl;
                        std::cout << m_common_klu.status << std::endl;
                        throw std::runtime_error("klu_analyze failed");
                     } //::TODO:: write real exception
                  else{m_symbolic_is_done=true;}
               }

            void refactor_nummeric(_SparseMatrix Matrix){
               m_ok = klu_refactor( Matrix.outerIndexPtr(),  Matrix.innerIndexPtr(),  Matrix.valuePtr(), m_symbolic_klu,
                           m_nummeric_klu, &m_common_klu);
                     if (m_ok == 0) throw std::runtime_error("klu_refactor failed");
                     //::TODO:: write real exception
            }


};

} /* namespace Eigen */
#endif /* KLUWRAPPER_HPP_ */

 
Reply to topic

Bookmarks



Who is online

Registered users: Abdelrahman Magdy, Baidu [Spider], Bing [Bot], colomar, cylverbak, edmael, Exabot [Bot], Google [Bot], google01103, jeanliu, ken300, koriun, Majestic-12 [Bot], MSNbot Media, navalgupte, paulus3005, psbot [Picsearch], robinh, scummos, Sentynel, TheraHedwig, tshearer, Uri_Herrera, vascobasque, wolfi323, Yahoo [Bot]