AFSymMatrix.hpp Source File

00001 #ifndef __AFSymMatrix_H__
00002 #define __AFSymMatrix_H__
00003 
00004 
00005 
00026 #include <stdio.h>
00027 #include "FloatVector.hpp"
00028 
00029 #include "AFloatMatrix.hpp"
00030 
00031 class FDiagonalMatrix;
00032 class AFProductSeqSMatrix;
00033 
00034 class AFSymMatrix : public AFloatMatrix
00035 {
00036 protected:
00037 
00038 int width;
00039 
00040 public:
00041 AFSymMatrix(int w){ width= w; }
00042 virtual ~AFSymMatrix(){}
00043 
00044 virtual AFSymMatrix* copyAnn()
00045                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0;}
00046 virtual AFSymMatrix* copyFull();
00047 
00048 virtual FDiagonalMatrix* copyDiagonal();
00049 // no diagonal
00050 virtual AFSymMatrix* copyL0();                          
00051 virtual AFSymMatrix* copyU0();
00052 // with diagonal
00053 virtual AFSymMatrix* copyL();                           
00054 virtual AFSymMatrix* copyU();
00055 // with unit diagonal
00056 virtual AFSymMatrix* copyL1();                          
00057 virtual AFSymMatrix* copyU1();
00058 
00059 virtual int getWidth() const {return width;}                            // use of abstract is better for dynamic block size,...
00060 virtual int getHeight() const {return width;}
00061 
00062 virtual void setAll(AFloatMatrix* from, bool t=false);
00063 virtual void setAll(float value=0.0);
00064 
00065 // NxN matrix operations  
00066 virtual AFSymMatrix* add_A_B(AFSymMatrix* B, AFSymMatrix* result=0);    // C= A+B
00067 virtual AFSymMatrix* add_B_A(AFSymMatrix* B, AFSymMatrix* result=0)     // C= B+A
00068                 { return B->add_A_B(this, result); }
00069 virtual AFSymMatrix* subst_A_B(AFSymMatrix* B, AFSymMatrix* result=0);  // C= A-B
00070 virtual AFSymMatrix* subst_B_A(AFSymMatrix* B, AFSymMatrix* result=0);  // C= B-A
00071 virtual AFSymMatrix* mult_A_B(AFSymMatrix* B, AFSymMatrix* result=0);   // C= A*B
00072 virtual AFSymMatrix* mult_B_A(AFSymMatrix* B, AFSymMatrix* result=0);   // C= B*A
00073 
00074 // A*x=b <==> 
00075 virtual FloatVector* solve(FloatVector* b, FloatVector* x=0){return 0;} // Ax=b, should be possible if SPD, ... 
00076 virtual AFSymMatrix* invert(){return 0;}                                                // A^-1, is easy if SPD
00077 
00078 
00079 virtual AFSymMatrix* jacobiMatrix();                                            // used to compute SOR omega
00080 virtual FloatVector* jacobiVector(FloatVector* b);                              // xk1= BJ*xk + bJ 
00081 
00082 // direct solvers
00083 virtual FloatVector* linearSolverU0(FloatVector* b, FloatVector* x=0);                  // Ux=b
00084 virtual FloatVector* linearSolverU0(FloatVector* b, int k, FloatVector* x);             // U(k)x(k)=b(k)
00085 virtual FloatVector* linearSolverU1(FloatVector* b, FloatVector* x=0);                  // U1x=b
00086 virtual FloatVector* linearSolverL0(FloatVector* b, FloatVector* x=0);                  // Lx=b
00087 virtual FloatVector* linearSolverL1(FloatVector* b, FloatVector* x=0);                  // L1x=b
00088 
00089 
00090 
00091 
00092 virtual AFSymMatrix* inverseSPD(AFSymMatrix* result=0);                         // SPD inversion
00093 
00094 // Combined solvers
00095 virtual FloatVector* linearSolverL1U0(FloatVector* b, FloatVector* x=0);        // crout LU
00096 virtual FloatVector* linearSolverL0U1(FloatVector* b, FloatVector* x=0);        
00097 
00098 virtual FloatVector* linearSolverLLt(FloatVector* b, FloatVector* x=0);         // cholesky solvers
00099 virtual FloatVector* linearSolverUUt(FloatVector* b, FloatVector* x=0);
00100 
00101 // Cholesky sans racine carree
00102 virtual FloatVector* linearSolverL1DL1t(FloatVector* b, FloatVector* x=0);
00103 virtual FloatVector* linearSolverU1DU1t(FloatVector* b, FloatVector* x=0);
00104 
00105 // factorisation complete
00106 virtual AFSymMatrix* factorisationL1U0();               // L1/U0 Doolittle
00107 virtual AFSymMatrix* factorisationL0U1();               // L0/U1 Crout
00108 
00109 virtual AFSymMatrix* factorisationLLt();                // LLt Cholesky
00110 virtual AFSymMatrix* factorisationUUt();                // UUt Cholesky
00111 virtual AFSymMatrix* factorisationL1DL1t();     // L1DL1t Cholesky sans racine carree
00112 virtual AFSymMatrix* factorisationU1DU1t();     // U1DU1t Cholesky sans racine carree
00113 
00114 virtual AFSymMatrix* choleskyUDUt(AFSymMatrix* result=0);               // in place UDUt Cholesky sans racine carree
00115 virtual AFSymMatrix* inverseUDUt(AFSymMatrix* result=0);                // in place UDUt Cholesky sans racine carree
00116 
00117 // incomplete factorisation & preconditioners
00118 // ------------------------------------------
00119 virtual AFloatMatrix* ilu0()
00120                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0; }
00121 virtual AFloatMatrix* milu0()
00122                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0; }
00123 virtual AFloatMatrix* rilu0(float omega=0.95)
00124                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0; }
00125 virtual AFloatMatrix* iluth(double threshold=0.001)
00126                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0; }
00127 virtual AFloatMatrix* ilut(int p, double eps)
00128                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0;}                // Threshold is related to magitude in row
00129 virtual AFloatMatrix* icf()
00130                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0; }
00131 virtual AFloatMatrix* ildlt()
00132                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__); return 0; }
00133 
00134 virtual FDiagonalMatrix* jacobiPrecond(FDiagonalMatrix* S=0);
00135 
00136 virtual void sor(double omega, AFSymMatrix** R=0, AFSymMatrix**T=0)
00137                         {printf("%s() %s: not Implemented\n",__FUNCTION__, __FILE__);} 
00138 
00139 // new SSOR
00140 virtual AFProductSeqSMatrix* ssorFactorisation();       // allows to select the typed SSOR
00141 
00142 // old SSOR way
00143 virtual void ssor(double omega, AFSymMatrix** R=0, FDiagonalMatrix** S=0, AFSymMatrix**T=0);
00144 virtual void ssorSPD(double omega, AFSymMatrix** R=0, FDiagonalMatrix** S=0);
00145 
00146 // Scaled matrices require diag(A)=I
00147 virtual void ssorScaled(double omega, float* cste, AFSymMatrix** R=0, AFSymMatrix**T=0);
00148 virtual void ssorScaledSPD(double omega, float* cste, AFSymMatrix** R=0);
00149 
00150 // approximate inversion
00151 virtual AFloatMatrix* apinv(int iterOuter, int iterInner){return 0;}            // column based iterative inversion
00152 
00153 // iterative simple linear solvers: Au=b
00154 // -------------------------------------
00155 virtual FloatVector* gaussSeidelLS(int iterMax, FloatVector* b, FloatVector* x=0);
00156 virtual FloatVector* jacobiLS(int iterMax, FloatVector* b, FloatVector* x=0);
00157 virtual FloatVector* sorLS(int iterMax, FloatVector* b, float omega=1.0, FloatVector* x=0);
00158 
00159 // single step iterative linear solvers: Au=b, single iteration
00160 virtual FloatVector* gaussSeidel(FloatVector* b, FloatVector* x=0);
00161 virtual FloatVector* jacobi(FloatVector* b, FloatVector* x=0, FloatVector* dest=0);
00162 virtual FloatVector* sor(FloatVector* b, float omega=1.0, FloatVector* x=0);
00163 
00164 // finding eigenvectors
00165 virtual FloatVector* power(int nmax, float* firstEigenValue, FloatVector* evMax0=0);
00166 virtual FloatVector* inverseIteration(double eigenvalue, int nmax=3, FloatVector* ev=0); 
00167 
00168 // find eigenvalues
00169 virtual float rayleighQuotient(FloatVector* ev=0);
00170 
00171 virtual void output();
00172 virtual void output(FILE* file);
00173 };
00174 
00175 #endif
00176 
00177