00001 00012 #ifndef TBCI_SOLVER_DILUPRECOND_H 00013 #define TBCI_SOLVER_DILUPRECOND_H 00014 00015 #include "precond.h" 00016 00017 NAMESPACE_TBCI 00018 00019 // ******************************************************************************** 00020 // ... For general BdMatrices 00021 // ******************************************************************************** 00022 00023 00024 template <typename T> class BdMatrix; 00025 00039 template <typename T> 00040 class DILU_BdMatrixPreconditioner : public Preconditioner_Sig<T, BdMatrix<T> > 00041 { 00042 public: 00043 DILU_BdMatrixPreconditioner () {} 00044 DILU_BdMatrixPreconditioner (const BdMatrix<T> &A) 00045 { update(A); } 00046 ~DILU_BdMatrixPreconditioner () {} 00047 00048 // update function 00049 void update (const BdMatrix<T> &A); 00050 00051 // solve functions 00052 TVector<T> solve (TVector<T> x) const; 00053 00054 inline TVector<T> solve (const Vector<T> &v) const 00055 { 00056 TVector<T> tv(v); 00057 return (solve(tv)); 00058 } 00059 00060 /* template <typename T> */ 00061 inline TVector<T> transSolve (const Vector<T> &v) const 00062 { return solve(v); } 00063 /* template <typename T> */ 00064 inline TVector<T> transSolve ( TVector<T> tv) const 00065 { return solve(tv); } 00066 00067 private: 00068 BVector<T> pivots; 00069 unsigned int dimension; 00070 }; 00071 00072 // See Templates Book, Fig. 3.3 00073 template <typename T> 00074 void DILU_BdMatrixPreconditioner<T>::update (const BdMatrix<T> &A) 00075 { 00076 // copy matrix variables to private vars 00077 dimension = A.dim; 00078 00079 // calculate Pivots 00080 if (pivots.size() != dimension) 00081 pivots.resize(dimension); 00082 BCHK(pivots.size() == 0, VecErr, "Null ptr in DILUBdMPrecCond::update", 0, ); 00083 const unsigned ds = A.diagconf.size(); 00084 unsigned i; 00085 for (i = 0; i < dimension; ++i) 00086 pivots(i) = A.diag[i]; 00087 for (i = 0; i < dimension; ++i) { 00088 register T tmp = pivots(i); 00089 if (UNLIKELY(tmp == T(0.0))) 00090 tmp = T(1.0); 00091 else 00092 tmp = T(1.0)/tmp; 00093 pivots(i) = tmp; 00094 for (unsigned j = 0; j < ds; j++) { 00095 const int df = A.diagconf(j); 00096 if (LIKELY(df < (signed)dimension - (signed)i)) 00097 pivots(i+df) -= A.adiag(df)[i]*tmp*A.bdiag(df)[i]; 00098 } 00099 } 00100 } 00101 00102 template <typename T> 00103 /*inline*/ TVector<T> DILU_BdMatrixPreconditioner<T>::solve (TVector<T> y) const 00104 { 00105 00106 // check dimensions 00107 BCHK(dimension != y.size(), VecErr, "DILU_BdMPreCond::solve: Dimension mismatch", y.size(), ); 00108 00109 // check trivial solution 00110 // Matrix consists of main diagonal only 00111 for (unsigned int i = 0; i < dimension; ++i) 00112 y.setval(i) *= pivots(i); 00113 return y; 00114 } 00115 00116 NAMESPACE_END 00117 00118 #endif /* TBCI_SOLVER_DILUPRECOND_H */
1.5.6