---

RealSMatrix.hpp

#ifndef __RealSMatrix_H__
#define __RealSMatrix_H__



#include <stdio.h>
#include "RealVector.hpp"
#include "H.hpp"                                // householder reflective matrix

class I;
class RTriDiagonalMatrix;

class RealSMatrix
{
private:

int                     width;
double**                data;

public:
        
RealSMatrix(int n, double value=0.0);
RealSMatrix(I& ident);
RealSMatrix(H& h);
RealSMatrix(const RealSMatrix& m);
RealSMatrix operator=(const RealSMatrix& m);

~RealSMatrix();

RealSMatrix* copy();
RealSMatrix* identity();

RealSMatrix* AtA(RealSMatrix* dest=0);

RealSMatrix t();

int size(){ return width*width; }

int getWidth()const { return width; }
int getHeight()const { return width; }

// base 0
void set0(int i, int j, double value){ data[i][j]= value; }
double get0(int i, int j)const { return data[i][j]; }

// base 1
void set(int i, int j, double value){ data[i-1][j-1]= value; }
double get(int i, int j)const { return data[i-1][j-1]; }

// computations
void setAll(double value=0.0);
void set(RealSMatrix* from);
void setIdentity();

// optimized transform (reduced cons)
RealSMatrix* multA_X(RealSMatrix* dest);                                        // X <-- A*X
 
void addA_kX(RealSMatrix& dest, double k=1.0);                          // A <-- A+k*X 
void addA_k_vvT(RealVector& v, double k=1.0);                           // A <-- A+k*(v*vT) 
double vTAv(RealVector& v);                                                     // quadratic form: vT*A*v 
double rayleighQuotient(RealVector& v){return vTAv(v)/v.vT_v();}        // r= vT*A*v/(vT*v) approximates an eigenvalue 

double det();
double trace();
double cond();

// norm
double normeMin();
double normeInfinite();
double norme2();                        // spectral norme= sqrt( max eigenvalue )
double normeF();                        // frobenius norme      

double sum();

double minimum();
double maximum();
double mean(){ return sum()/size(); }
double sigma();

// matrix generation
RealSMatrix* jacobiEV(RealVector* dEigenValue, RealSMatrix* vEigenVector); // Numerical Recipes p365-6

// math vectorial op
RealSMatrix sqrt();
RealSMatrix pow(double pow);
RealSMatrix sqr();
RealSMatrix exp();
RealSMatrix log();
RealSMatrix log(double base);

RealSMatrix sin();
RealSMatrix cos();
RealSMatrix tan();      

// matrix operations
// =================

void HA_inplace(H& h);          // A= H * A             low cost memory
void HA_inplace(H& h, int nI);          // A= H * A             low cost memory

void AH_inplace(H& h);                                  // A= A * H             low cost memory
void AH_inplace(H& h, int nI);                  // A= A * H             low cost memory
void AH_inplace(H& h, int nI, int ny);          // A= A * H             low cost memory

RealSMatrix* qrHouseholderR(RealSMatrix* r=0);          // QR upper Triangulation

void qrHouseholder(RealSMatrix* q, RealSMatrix* r);
void qrGivens(RealSMatrix* q, RealSMatrix* r);
void qrCGS(RealSMatrix* q, RealSMatrix* r);
void qrMGS(RealSMatrix* q, RealSMatrix* r);
RealSMatrix* upperHessenberg(RealSMatrix* uH=0);

// Penrose pseudoinverse via QR methods
RealSMatrix* pseudoInverse(RealSMatrix* x=0);
RealSMatrix* varianceCovariance(RealSMatrix* aT_aInverse=0);    // X= (At*A)^(-1)

// Neuronic and eigenvectors of a matrix
RealVector* oja(int nmax=100, double beta=0.1, RealVector* oldFirstPCA=0);      // Find First principal component or ev.
RealSMatrix* mOja(double nmax=100, double beta=0.1, RealSMatrix* oldPCA=0);     // Find all principal component/eigenvectors
RealSMatrix* apex(double nmax=100, double beta=0.1, RealSMatrix* oldPCA=0);     // Find all principal component/eigenvectors

// Neuronic and eigenvalues of a matrix: lambdai= E( (wi*x)^2 )= var( wi*x )
double ojaEigenvalue(int nmax=100, int ns=10, double beta=0.1, RealVector* oldFirstPCA=0);      // Evaluates first eigenvalue.
RealVector* apexEigenvalues(int nmax=100, int ns=10, double beta=0.1, RealSMatrix* oldPCA=0);   // Evaluates all eigenvalues.

// Simple PCA: SPCA1 and SPCA2 (M. Partridge, University of Sydney Australia, 1997)
RealVector* firstPcSPCA(int type=2, bool batch=false, int nmax=100, RealVector* oldFirstPCA=0); // Find First principal component or ev.
double firstEigenvalueSPCA(int type=2, bool batch=false, int nmax=100, RealVector* oldFirstPCA=0);      // Find First principal component or ev.

RealVector* kiemePcSPCA(int type=2, int k=0, bool batch=false, int nmax=100, RealVector* oldFirstPCA=0);        // Find kieme principal component or ev.
double kiemeEigenvalueSPCA(int type=2, int k=0, bool batch=false, int nmax=100, RealVector* oldFirstPCA=0);     // Find kieme principal component or ev.

RealSMatrix* fullSPCA(int type=2, bool batch=false, double nmax=100, RealSMatrix* oldPCA=0);    // Find all principal component/eigenvectors
RealVector* eigenvaluesSPCA(int type=2, bool batch=false, double nmax=100, RealSMatrix* oldPCA=0);      // Find all principal component/eigenvectors

// eigenvectors
RealVector* inverseIteration(double eigenvalue, int nmax=3, RealVector* ev=0); 

// eigenvalues
RealSMatrix* bidiagonalisation(RealSMatrix* biDiag=0);
RealSMatrix* svdGolubKahan(RealSMatrix* biDiag=0);

RealSMatrix* qrIteration(int iterMax=3, RealSMatrix* Ak=0);                     // simple QR iteration for eigenvalues O(N^4)
RealSMatrix* qrHessenbergIteration(int iterMax=3, RealSMatrix* r=0);    // Hessenberg QR converge faster and is O(N^3)
RealSMatrix* qrShift(int iterMax=3, RealSMatrix* r=0);
RealSMatrix* qrDoubleShift(int iterMax=3, RealSMatrix* r=0);            // Can support complexes via double iterations
//RealSMatrix* qrImplicit(int iterMax=3, RealSMatrix* r=0);

RealSMatrix* qrSymmetric(int iterMax=3, RealSMatrix* r=0);      // A symmetric=>Hu is tridiagonal & eigenvalues are all real
RTriDiagonalMatrix* lanczosSymmetric(RTriDiagonalMatrix* t=0);  // Generates cheap large tridiagonals
RTriDiagonalMatrix* lanczosMGS(RTriDiagonalMatrix* t=0);        // Tridiagonalise easier with modified Gram-Schmidt Lanczos

// Factorisation
void luDoolittle(RealSMatrix* l1, RealSMatrix* u);              // A -> L1, U
void luCrout(RealSMatrix* l, RealSMatrix* u1){}                 // A -> L, U1

// incomplete factorisation
RealSMatrix* ilu0();
RealSMatrix* iluth(double threshold=0.001);
//RealSMatrix* iluk();

// linear solvers A*x=b
RealVector* backwardSubstitution(RealVector& b, RealVector* x=0);               // U*x=b
RealVector* backwardSubstitutionU1(RealVector& b, RealVector* x=0);     // U1*x=b
RealVector* forwardElimination(RealVector& b, RealVector* x=0);         // L*x=b
RealVector* forwardEliminationL1(RealVector& b, RealVector* x=0);               // L1*x=b

RealVector* linearSolverU(RealVector& b, RealVector* x=0);
RealVector* linearSolverU1(RealVector& b, RealVector* x=0);
RealVector* linearSolverL(RealVector& b, RealVector* x=0);
RealVector* linearSolverL1(RealVector& b, RealVector* x=0);

// U, L
RealSMatrix* inverseOfU(RealSMatrix* uinv=0);           // Uinv = U^-1
RealSMatrix* inverseOfU1(RealSMatrix* uinv=0);          // U1inv = U1^-1

// scalar single operations
void operator+=(double value);
void operator-=(double value);
void operator*=(double value);
void operator/=(double value);

// matrix single operations
void operator+=(const RealSMatrix& m);
void operator-=(const RealSMatrix& m);
void operator*=(const RealSMatrix& m);

// friends
friend RealVectorT operator*(const RealVectorT& vt, const RealSMatrix& m);      
friend RealVector operator*(const RealSMatrix& m, const RealVector& v); 
        
void output();
void output(FILE* file);
};

#endif


         

Restoreinpaint sourceforge project `C++/Java Image Processing, Restoration, Inpainting Project'.

Bernard De Cuyper: Open Project Leader: Concept, design and development. Bernard De Cuyper & Eddy Fraiha 2002, 2003. Bernard De Cuyper 2004. Open and free, for friendly usage only. Modifications on Belgium ground of this piece of artistic work, by governement institutions or companies, must be notified to Bernard De Cuyper.