[Toon-members] TooN Lapack_Cholesky.h

[Tom Drummond]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Tom Drummond <twd20>    09/04/22 08:05:10

Added files:
        .              : Lapack_Cholesky.h 

Log message:
        added lapack cholesky decomposition

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Lapack_Cholesky.h?cvsroot=toon&rev=1.1

Patches:
Index: Lapack_Cholesky.h
===================================================================
RCS file: Lapack_Cholesky.h
diff -N Lapack_Cholesky.h
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ Lapack_Cholesky.h   22 Apr 2009 08:05:09 -0000      1.1
@@ -0,0 +1,168 @@
+// -*- c++ -*-
+
+//     Copyright (C) 2009 Tom Drummond (address@hidden)
+//
+// This file is part of the TooN Library.  This library is free
+// software; 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, or (at your option)
+// any later version.
+
+// This library 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 General Public License for more details.
+
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING.  If not, write to the Free
+// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
+// USA.
+
+// As a special exception, you may use this file as part of a free software
+// library without restriction.  Specifically, if other files instantiate
+// templates or use macros or inline functions from this file, or you compile
+// this file and link it with other files to produce an executable, this
+// file does not by itself cause the resulting executable to be covered by
+// the GNU General Public License.  This exception does not however
+// invalidate any other reasons why the executable file might be covered by
+// the GNU General Public License.
+
+
+#ifndef TOON_INCLUDE_CHOLESKY_H
+#define TOON_INCLUDE_CHOLESKY_H
+
+#include <TooN/TooN.h>
+
+namespace TooN {
+
+
+/**
+Decomposes a positive-semidefinite symmetric matrix A (such as a covariance) 
into L*L^T, where L is lower-triangular.
+Also can compute A = S*S^T, with S lower triangular.  The LDL^T form is faster 
to compute than the class Cholesky decomposition.
+The decomposition can be used to compute A^-1*x, A^-1*M, M*A^-1*M^T, and A^-1 
itself, though the latter rarely needs to be explicitly represented.
+Also efficiently computes det(A) and rank(A).
+It can be used as follows:
address@hidden
+// Declare some matrices.
+Matrix<3> A = ...; // we'll pretend it is pos-def
+Matrix<2,3> M;
+Matrix<2> B;
+Vector<3> y = make_Vector(2,3,4);
+// create the Cholesky decomposition of A
+Cholesky<3> chol(A);
+// compute x = A^-1 * y
+x = cholA.backsub(y);
+//compute A^-1
+Matrix<3> Ainv = cholA.get_inverse();
address@hidden
address@hidden gDecomps
+
+Cholesky decomposition of a symmetric matrix.
+Only the lower half of the matrix is considered
+This uses the non-sqrt version of the decomposition
+giving symmetric M = L*D*L.T() where the diagonal of L contains ones
address@hidden Size the size of the matrix
address@hidden Precision the precision of the entries in the matrix and its 
decomposition
+**/
+template <int Size, tyepname Precision>
+class Lapack_Cholesky {
+public:
+
+    Lapack_Cholesky(){}
+       
+       template<class P2, class B2>
+       Lapack_Cholesky(const Matrix<Size, Size, P2, B2>& m) {
+               : my_cholesky(m) {
+               SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
+               do_compute();
+       }
+
+       /// Constructor for Size=Dynamic
+       Lapack_Cholesky(int size) : my_cholesky(size,size) {}
+
+       template<class P2, class B2> void compute(const Matrix<Size, Size, P2, 
B2>& m){
+               SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
+               SizeMismatch<Size,Size>::test(m.num_rows(), 
my_cholesky.num_rows());
+               my_cholesky=m;
+               do_compute();
+       }
+
+
+
+       void do_compute(){
+               int N = my_cholesky.num_rows();
+               int info;
+               dpotrf_("L", &N, L.get_data_ptr(), &N, &info);
+               assert(info >= 0);
+               if (info > 0) {
+                       my_rank = info-1;
+               }
+       }
+
+       int rank() const { return my_rank; }
+
+       template <int Size2, typename P2, typename B2>
+               Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>& 
v) {
+               SizeMismatch<Size,Size2>::test(my_cholesky.num_cols(), 
v.size());
+
+               Vector<Size> result(v);
+               int N=L.num_rows();
+               int NRHS=1;
+               int info;
+               dpotrs_("L", &N, &NRHS, my_cholesky.my_data, &N, 
result.my_data, &N, &info);     
+               assert(info==0);
+               return result;
+       }
+
+       template <int Size2, int Cols2, typename P2, typename B2>
+               Matrix<Size, Cols2, Precision, ColMajor> backsub (const 
Matrix<Size2, Cols2, P2, B2>& m) {
+               SizeMismatch<Size,Size2>::test(my_cholesky.num_cols(), 
m.num_rows());
+
+               Matrix<Size, Cols2, Precision, ColMajor> result(m);
+               int N=my_cholesky.num_rows();
+               int NRHS=m.num_cols();
+               int info;
+               dpotrs_("L", &N, &NRHS, my_cholesky.my_data, &N, 
result.my_data, &N, &info);     
+               assert(info==0);
+               return result;
+       }
+
+
+
+
+
+       template <int Size2, typename P2, typename B2>
+               Precision mahalanobis(const Vector<Size2, P2, B2>& v) const {
+               return v * backsub(v);
+       }
+
+       const Matrix<>& get_L() const {
+               return my_cholesky;
+       }
+
+       Precision determinant() const {
+               Precision det = L[0][0];
+               for (int i=1; i<my_cholesky.num_rows(); i++)
+                       det *= L[i][i];
+               return det*det;
+       }
+
+       Matrix<> get_inverse() const {
+               Matrix<Size> M(my_cholesky);
+               int N = my_cholesky.num_rows();
+               int info;
+               dpotri_("L", &N, M.my_data, &N, &info);
+               assert(info == 0);
+               TooN::Symmetrize(M);
+               return M;
+       }
+
+private:
+       Matrix<Size,Size,Precision> my_cholesky;     
+       int rank;
+};
+
+
+}
+
+#endif