[prev in list] [next in list] [prev in thread] [next in thread] 

List:       kde-commits
Subject:    kdesupport/eigen
From:       Benoît Jacob <jacob () math ! jussieu ! fr>
Date:       2006-08-31 21:26:49
Message-ID: 1157059609.959935.1192.nullmailer () svn ! kde ! org
[Download RAW message or body]

SVN commit 579344 by bjacob:

fixed bugs introduced in recent changes

M    eigen/kvectorn.h
M    eigen/kmatrixn.h
M    eigen/testkmatrixn.cpp
M    eigen/kludecomposition.h


 M  +65 -20    kludecomposition.h  
 M  +1 -1      kmatrixn.h  
 M  +6 -2      kvectorn.h  
 M  +5 -6      testkmatrixn.cpp  


--- trunk/kdesupport/eigen/kludecomposition.h #579343:579344
@@ -60,6 +60,15 @@
 template <class T>
 class KLUDecomposition
 {
+
+protected:
+
+    /**
+      * The method actually computing the LU decomposition.
+      * Called by the constructor.
+      */
+    void perform ( const KMatrixN<T>& A );
+
 public:
 
     /**
@@ -67,7 +76,8 @@
       *
       * @param A the square matrix to decompose
       */
-    KLUDecomposition( const KMatrixN<T>& A );
+    inline KLUDecomposition( const KMatrixN<T>& A )
+    { perform( A ); }
 
     ~KLUDecomposition();
 
@@ -84,7 +94,14 @@
     T determinant() const;
 
     /**
-      * This method returns the dimension of the kernel of the square matrix A
+      * Returns the dimension (size) of square matrix A
+      * of which *this is the LU decomposition -- just in case you forgot it!
+      */
+    inline unsigned int dim() const
+    { return m_dim; }
+
+    /**
+      * Returns the dimension of the kernel of the square matrix A
       * of which *this is the LU decomposition. It is very fast because the
       * computation has already been done during the LU decomposition.
       */
@@ -92,7 +109,7 @@
     { return m_dimKer; }
 
     /**
-      * This method returns the rank of the square matrix A of which *this is
+      * Returns the rank of the square matrix A of which *this is
       * the LU decomposition. It is very fast because the
       * computation has already been done during the LU decomposition.
       */
@@ -100,7 +117,7 @@
     { return( m_dim - m_dimKer ); }
 
     /**
-      * This method returns true if the square matrix A, of which *this is
+      * Returns true if the square matrix A, of which *this is
       * the LU decomposition, is invertible. It returns false if it is singular.
       * It is very fast because the computation has already been done during
       * the LU decomposition.
@@ -108,14 +125,25 @@
     inline bool isInvertible() const
     { return( m_dimKer == 0 ); }
 
-//    KMatrixN solveLinearSystem( const 
+    /**
+      * Use this method to compute a basis of the kernel of the matrix A of
+      * which *this is the LU decomposition.
+      *
+      * @returns null if the kernel has dimension 0, that is, if A is
+      * invertible. Otherwise returns a pointer to a matrix whose column vectors
+      * form a basis of the kernel of A.
+      *
+      * Note: the caller doesn't have to delete the returned pointer. This class
+      * remembers the it and will delete it in the destructor.
+      */
+    KMatrixN<T> * basisKer();
 
     /**
       * Not yet implemented
       */
     KMatrixN<T> inverseMatrix();
 
-protected:
+//protected:
 
     /**
       * The dimension of the matrix A of which *this is the LU decomposition
@@ -168,23 +196,9 @@
       * The eigenvalue of U that has biggest absolute value.
       */
     T m_biggestEigenValueU;
-
-    /**
-      * The method actually computing the LU decomposition.
-      * Called by the constructor.
-      */
-    void perform ( const KMatrixN<T>& A );
 };
 
-
-
 template<class T>
-KLUDecomposition<T>::KLUDecomposition( const KMatrixN<T> & A )
-{
-    perform(A);
-}
-
-template<class T>
 void KLUDecomposition<T>::perform( const KMatrixN<T> & A )
 {
     //assert( A.width() == A.height() );
@@ -361,6 +375,37 @@
     return ret;
 }
 
+#if 0
+template<class T>
+KMatrixN<T> * KLUDecomposition<T>::basisKer()
+{
+    // remember that m_basisKer has been set to 0 by perform()
+    if( isInvertible() // if the matrix A is invertible, we want to return 0
+      || m_basisKer ) // if the computation has already been done
+    {
+        return m_basisKer;
+    }
+
+    // if we're here, we know that m_basisKer == 0 and that m_dimKer > 0
+    // (the latter is because isInvertible() returned false)
+
+    KMatrixN<T> m_basisKer = new KMatrixN<T> ( m_dim, m_dimKer );
+
+    /* now we use the following lemma:
+     *
+     * Lemma: If the matrix A has the LU decomposition PAQ = LU,
+     * then Ker A = Q( Ker U ).
+     *
+     * Proof: trivial: just keep in mind that P, Q, L are invertible.
+     */
+
+    // thus, all we need to do is to compute Ker U, and then multiply by Q.
+
+    
+
+}
+#endif
+
 } // namespace eigen
 
 #endif // KLUDECOMPOSITION_H
--- trunk/kdesupport/eigen/kmatrixn.h #579343:579344
@@ -294,7 +294,7 @@
     if( width * height > m_width * m_height )
     {
         delete[] m_array;
-        m_array  = new T[m_width * m_height];
+        m_array  = new T[width * height];
     }
     m_width  = width;
     m_height = height;
--- trunk/kdesupport/eigen/kvectorn.h #579343:579344
@@ -137,8 +137,12 @@
 void KVectorN<T>::resize( unsigned int size )
 {
     if( size == m_size ) return;
-    delete[] m_array;
-    init( size );
+    if( size > m_size )
+    {
+        delete[] m_array;
+        m_array  = new T[size];
+    }
+    m_size = size;
 }
 
 template<class T>
--- trunk/kdesupport/eigen/testkmatrixn.cpp #579343:579344
@@ -18,7 +18,7 @@
 using namespace eigen;
 
 
-/*
+
 void print ( const std::complex<double> c)
 {
     printf("%g+%gi\n", c.real(), c.imag() );
@@ -48,7 +48,7 @@
         printf("\n");
     }
     printf("---------------------------------\n");
-}*/
+}
 
 /*
 class MyThread : public QThread
@@ -265,10 +265,9 @@
 */
 
 
-  /*  KMatrixN<float> a;
+    KMatrixN<double> a;
     a.loadIdentity(5);
-    a = a* a;
-    printf("determinant of a: %.2f\n", a.determinant());
-*/
+    KLUDecomposition<double> lu(a);
+    print(*(lu.m_LU));
     return 1;
 }
[prev in list] [next in list] [prev in thread] [next in thread]