src/linalg.h

/*
 *   See COPYING file distributed along with the psignifit package for
 *   the copyright and license terms
 */
#ifndef LINALG_H
#define LINALG_H

#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>

/***************************************************************************************
 * This file defines vary basic linear algebra facilities
 *
 * These are not particularly fast or fancy implementations. They are just intended
 * to keep the number of dependencies low.
 */

class MatrixError
{
};

class Matrix
{
        private:
                double *data;
                unsigned int nrows;
                unsigned int ncols;
                // given a decomposition A=LU these two methods help solving the equation Ax=b:
                std::vector<double> forward ( const Matrix* LU, const std::vector<double>& b );    // forward solution of Ly=b
                std::vector<double> backward ( const Matrix*LU, const std::vector<double>& y );    // backward solution of Ux=y
        public:
                Matrix ( const std::vector< std::vector<double> >& A );      
                Matrix ( unsigned int nrows, unsigned int ncols);                   
                Matrix ( const Matrix& A );                                  
                ~Matrix ( void ) { delete [] data; }                         
                double& operator() ( unsigned int i, unsigned int j ) const; 
                void print ( void );                                         
                unsigned int getnrows ( void ) const { return nrows; }       
                unsigned int getncols ( void ) const { return ncols; }       
                Matrix* cholesky_dec ( void ) const;                         
                Matrix* lu_dec ( void ) const;                               
                Matrix * qr_dec ( void ) const;                              
                Matrix * inverse_qr ( void ) const;                          
                Matrix * regularized_inverse ( double alpha ) const;         
                std::vector<double> solve ( const std::vector<double>& b );  
                Matrix* inverse ( void );                                    
                std::vector<double> operator* ( std::vector<double>& x );    
                void scale ( double a );                                     
                bool symmetric ( void );                                     
};

std::vector<double> leastsq ( const Matrix *A, const std::vector<double>& b );  

std::vector<double> leastsq ( const Matrix *M ); 

#endif