---

FTriDiagonalMatrix.hpp

#ifndef __FTriDiagonalMatrix_H__
#define __FTriDiagonalMatrix_H__



#include <stdio.h>
#include "FloatVector.hpp"

#include "AFSymMatrix.hpp"

class FTriDiagonalMatrix : public AFSymMatrix
{
private:
        
int                             N;

FloatVector*            up;                                     // N-1 size upper diagonal 
FloatVector*            d;                                      // N size diagonal
FloatVector*            low;                            // N-1 size lower diagonal 

public:
        
FTriDiagonalMatrix(int asize, float val=0.0);
FTriDiagonalMatrix(int asize, float lval, float dval, float uval);
FTriDiagonalMatrix(const FTriDiagonalMatrix& v);
FTriDiagonalMatrix operator=(const FTriDiagonalMatrix& v);

virtual ~FTriDiagonalMatrix(){ delete d; delete up; delete low;}

virtual AFloatMatrix* copy(){ return new FTriDiagonalMatrix(*this); }
virtual AFloatMatrix* t(AFloatMatrix* result=0){return (AFloatMatrix*)0;}

// no diagonal
virtual AFSymMatrix* copyL0();                          
virtual AFSymMatrix* copyU0();
// with diagonal
virtual AFSymMatrix* copyL();                           
virtual AFSymMatrix* copyU();
// with unit diagonal
virtual AFSymMatrix* copyL1();                          
virtual AFSymMatrix* copyU1();

void load(FTriDiagonalMatrix* m);
void load(FloatVector* u1, FloatVector* d1, FloatVector* l1);

virtual int size()const { return (3*N-2); }

FloatVector* getUvalues(){ return up; }
FloatVector* getDvalues(){ return d; }
FloatVector* getLvalues(){ return low; }

// base 0
virtual void set0(int i, int j, float value);
virtual float get0(int i, int j);

inline void setD0(int i, float value){ d->set0(i, value); }
inline float getD0(int i)const { return d->get0(i); }

inline void setU0(int i, float value){ up->set0(i, value); }
inline float getU0(int i)const { return up->get0(i); }

inline void setL0(int i, float value){ low->set0(i, value); }
inline float getL0(int i)const { return low->get0(i); }

virtual Simple1DIndexList* indexesInRow0(int k, Simple1DIndexList* oldRow=0);
virtual Simple1DIndexList* indexesInCol0(int k, Simple1DIndexList* oldCol=0);

// base 1
virtual void set(int i, int j, float value);
virtual float get(int i, int j);

inline void setD(int i, float value){ d->set(i, value); }
inline float getD(int i)const { return d->get(i); }

inline void setU(int i, float value){ up->set(i, value); }
inline float getU(int i)const { return up->get(i); }

inline void setL(int i, float value){ low->set(i, value); }
inline float getL(int i)const { return low->get(i); }

virtual Simple1DIndexList* indexesInRow(int k, Simple1DIndexList* oldRow=0);
virtual Simple1DIndexList* indexesInCol(int k, Simple1DIndexList* oldCol=0);

// direct solvers
virtual FloatVector* linearSolverU0(FloatVector* b, FloatVector* x=0);          // Ux=b
virtual FloatVector* linearSolverU1(FloatVector* b, FloatVector* x=0);          // U1x=b
virtual FloatVector* linearSolverL0(FloatVector* b, FloatVector* x=0);          // Lx=b
virtual FloatVector* linearSolverL1(FloatVector* b, FloatVector* x=0);          // L1x=b

// factorisation complete
virtual AFSymMatrix* factorisationL1U0();               // L1/U0 Doolittle
virtual AFSymMatrix* factorisationL0U1();               // L0/U1 Crout

// single step iterative linear solvers: Au=b, single iteration
virtual FloatVector* gaussSeidel(FloatVector* b, FloatVector* x=0);
virtual FloatVector* jacobi(FloatVector* b, FloatVector* x=0, FloatVector* dest=0);
virtual FloatVector* sor(FloatVector* b, float omega=1.0, FloatVector* x=0);

virtual AFSymMatrix* jacobiMatrix();                                    // used to compute SOR omega

// computations

virtual void setAll(float value=0.0);

virtual float det();
virtual float trace();

virtual float norme2();
virtual float sum();

virtual float minimum();
virtual float maximum();
virtual float sigma();

// math vectorial op
FTriDiagonalMatrix sqrt();
FTriDiagonalMatrix pow(float pow);
FTriDiagonalMatrix sqr();
FTriDiagonalMatrix exp();
FTriDiagonalMatrix log();
FTriDiagonalMatrix log(float base);

FTriDiagonalMatrix sin();
FTriDiagonalMatrix cos();
FTriDiagonalMatrix tan();

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

virtual void add(float value);
virtual void subst(float value);
virtual void mult(float value);
virtual void div(float value);

// vector single operations
void operator+=(const FTriDiagonalMatrix& m);
void operator-=(const FTriDiagonalMatrix& m);

virtual void add(AFloatMatrix& m){}
virtual void subst(AFloatMatrix& m){}

// u= A*v including implicit transposition
FloatVector* mult_Av(FloatVector* v, FloatVector* result=0);
FloatVector* mult_ATv(FloatVector* v, FloatVector* result=0);
// u= v*A
FloatVector* mult_vA(FloatVector* v, FloatVector* result=0)     { return mult_ATv(v, result); }
FloatVector* mult_vAT(FloatVector* vt, FloatVector* result=0){ return mult_Av(vt, result); }

// sub vector usage
// u= A*v
FloatVector* mult_Av(FloatVector* v, int col0, int row0, bool incremental=false, FloatVector* result=0);
FloatVector* mult_ATv(FloatVector* v, int col0, int row0, bool incremental=false, FloatVector* result=0);
// u= v*A
FloatVector* mult_vA(FloatVector* v, int row0, int col0, bool incremental=false, FloatVector* result=0)
                        { return mult_ATv(v, col0, row0, incremental, result); }
FloatVector* mult_vAT(FloatVector* vt, int row0, int col0, bool incremental=false, FloatVector* result=0)
                        { return mult_Av(vt, col0, row0, incremental, result); }

// friends
friend FloatVectorT operator*(const FloatVectorT& vt, const FTriDiagonalMatrix& m);     
friend FloatVector operator*(const FTriDiagonalMatrix& m, const FloatVector& 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.