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

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.

Ok,

thanks. I will try that when I am back.

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 ?

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()

I found my error: test.size is the number of elements but I need the number of columns /rows.

Sorry and thanks.

Feel free to share your KLU wrapper once it's ready! I'm sure it will be useful to many people.

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_ */