---

FBandSMatrix.hpp

#ifndef __FBandSMatrix_H__
#define __FBandSMatrix_H__


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

#include "AFSymMatrix.hpp"

class FBandSMatrix : public AFSymMatrix
{
protected:

int                             n, middle, nmax;                                // n bands

int                             nL0, nD, nU0;

int*                            index;                          // 0=diagonal, <0 L, >0 U
FloatVector**           band;                                   // N-1 size upper diagonal 

public:
        
FBandSMatrix(int asize, int nband, int* iband, float val=0.0);
FBandSMatrix(int asize, int nband, int* iband, float* val);
FBandSMatrix(const FBandSMatrix& v);
FBandSMatrix operator=(const FBandSMatrix& v);

virtual ~FBandSMatrix();

virtual AFloatMatrix* copy();
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();

virtual int numberOfBands(){ return n; }
virtual int numberOfL0Bands(){ return nL0; }
virtual int numberOfDBands(){ return nD;}
virtual int numberOfU0Bands(){ return nU0; }

virtual int getIndexBand(int i){ return index[i]; }

virtual int size()const;

FloatVector* getBand(int i){ return band[i]; }
FloatVector* getUvalues(int i){ return band[middle+i]; }
FloatVector* getDvalues(){ return band[middle]; }
FloatVector* getLvalues(int i){ return band[middle-i]; }

virtual void setBand(int i, FloatVector* v){band[i]->load(v);}

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

virtual void setD0(int i, float value);
virtual float getD0(int i)const;

virtual void setU0(int iband, int i, float value);
virtual float getU0(int iband, int i)const;

virtual void setL0(int iband, int i, float value);
virtual float getL0(int iband, int i)const;

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){ set0(i-1,j-1, value); }
virtual float get(int i, int j){ return get0(i-1,j-1); }

virtual void setD(int i, float value);
virtual float getD(int i)const;

virtual void setU(int iband, int i, float value);
virtual float getU(int iband, int i)const;

virtual void setL(int iband, int i, float value);
virtual float getL(int iband, int i)const;

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

// direct solvers: optimal band Ax=b, O(nband*N)
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

// 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();


// 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

// scalar single operations
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);

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); }


// incomplete factorisation
virtual AFloatMatrix* ilu0();
virtual AFloatMatrix* iluth(double threshold=0.001);

virtual void output(){ printf("Band(%d) ", numberOfBands()); AFSymMatrix::output(); }
};

#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.