matrixMat.h Source File

00001 /*=============================================================================
00002         File: matrixMat.h
00003      Purpose:       
00004     Revision: $Id: matrixMat.h,v 1.2 2002/05/13 21:07:45 philosophil Exp $
00005   Created by: Philippe Lavoie          (22 Oct, 1997)
00006  Modified by: 
00007 
00008  Copyright notice:
00009           Copyright (C) 1996-1997 Philippe Lavoie
00010  
00011           This library is free software; you can redistribute it and/or
00012           modify it under the terms of the GNU Library General Public
00013           License as published by the Free Software Foundation; either
00014           version 2 of the License, or (at your option) any later version.
00015  
00016           This library is distributed in the hope that it will be useful,
00017           but WITHOUT ANY WARRANTY; without even the implied warranty of
00018           MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019           Library General Public License for more details.
00020  
00021           You should have received a copy of the GNU Library General Public
00022           License along with this library; if not, write to the Free
00023           Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00024 =============================================================================*/
00025 #ifndef _matrix_matrixMat_h_
00026 #define _matrix_matrixMat_h_
00027 
00028 #include "matrix.h"
00029 
00032 namespace PLib {
00033 
00061   template <class T>
00062   class LUMatrix : public Matrix<T> 
00063   {
00064   public:
00065     LUMatrix(int r, int c) : Matrix<T>(r,c),pivot(pivot_) { pivot_.resize(r) ; }
00066     LUMatrix()  : Matrix<T>(),pivot(pivot_)    {; }
00067     LUMatrix(const LUMatrix<T>& lu): Matrix<T>(lu),pivot(pivot_) { pivot_ = lu.pivot_;}
00068     LUMatrix(const Matrix<T>& a): Matrix<T>(a.rows(),a.cols()),pivot(pivot_) { decompose(a) ; }
00069     
00070     void resize(const int r, const int c) { Matrix<T>::resize(r,c) ; pivot_.resize(r) ; }
00071     LUMatrix& operator=(const LUMatrix<T>&);
00072     LUMatrix& decompose(const Matrix<T> &a);
00073     T determinant() ;
00074     
00075     Matrix<T> inverse();
00076     void inverseIn(Matrix<T>&);
00077     const Vector<int>& pivot ;
00078     
00079   private:
00080     Vector<int> pivot_ ;
00081     
00082   protected:
00083     int errval ;
00084     int sign ;
00085   };
00086   
00087   
00178   template <class T>
00179   class SVDMatrix 
00180   {
00181   public:
00182     SVDMatrix():U(U_), V(V_), sig(sig_) { ; }
00183     SVDMatrix(const Matrix<T>& A);
00184     
00185     const Matrix<T>& U ;
00186     const Matrix<T>& V ;
00187     const Vector<T>& sig ;
00188     
00189     int decompose(const Matrix<T>& A) ;
00190     void minMax(T& min_sig, T& max_sig) const;
00191     double q_cond_number(void) const;   
00192     
00193     void cut_off(const double min_sig);
00194     void inverseIn(Matrix<T>& inv, double tau=0) ;
00195     Matrix<T> inverse(double tau=0) ;
00196     int solve(const Matrix<T>& B, Matrix<T>& X, double tau=0) ;
00197     
00198   protected:
00199     int M,N;    
00200     Matrix<T> U_;         
00201     Matrix<T> V_;         
00202     Vector<T> sig_; 
00203     
00204     // Internal procedures used in SVD
00205     double left_householder(Matrix<T>& A, const int i);
00206     double right_householder(Matrix<T>& A, const int i);
00207     double bidiagonalize(Vector<T>& super_diag, const Matrix<T>& _A);
00208     
00209     void rotate(Matrix<T>& U, const int i, const int j,
00210                 const double cos_ph, const double sin_ph);
00211     void rip_through(Vector<T>& super_diag, const int k, const int l, const double eps);
00212     int get_submatrix_to_work_on(Vector<T>& super_diag, const int k, const double eps);
00213     void diagonalize(Vector<T>& super_diag, const double eps);
00214     
00215   };
00216 
00217   template <class T> int solve(const Matrix<T>& A, const Matrix<T>& B, Matrix<T>& C) ;
00218   template <class T> Matrix<T> inverse(const Matrix<T>& A) ;
00219 
00220 } // end namespace
00221 
00222 
00223 #ifdef INCLUDE_TEMPLATE_SOURCE
00224 #include "matrixMat.cpp"
00225 #endif
00226 
00227 
00228 
00229 
00230 #endif