00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__shifted_matrix_h
00013 #define __deal2__shifted_matrix_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/smartpointer.h>
00018
00019 DEAL_II_NAMESPACE_OPEN
00020
00034 template<class MATRIX>
00035 class ShiftedMatrix
00036 {
00037 public:
00042 ShiftedMatrix (const MATRIX& A, const double sigma);
00043
00047 void shift (const double sigma);
00048
00052 double shift () const;
00053
00057 template <class VECTOR>
00058 void vmult (VECTOR& dst, const VECTOR& src) const;
00059
00063 template <class VECTOR>
00064 double residual (VECTOR& dst, const VECTOR& src, const VECTOR& rhs) const;
00065
00066 private:
00070 SmartPointer<const MATRIX,ShiftedMatrix<MATRIX> > A;
00071
00075
00079 double sigma;
00080 };
00081
00082
00083
00093 template<class MATRIX, class MASSMATRIX, class VECTOR>
00094 class ShiftedMatrixGeneralized
00095 {
00096 public:
00101 ShiftedMatrixGeneralized (const MATRIX& A,
00102 const MASSMATRIX& M,
00103 const double sigma);
00104
00108 void shift (const double sigma);
00109
00113 double shift () const;
00114
00118 void vmult (VECTOR& dst, const VECTOR& src) const;
00119
00123 double residual (VECTOR& dst, const VECTOR& src, const VECTOR& rhs) const;
00124
00125 private:
00129 SmartPointer<const MATRIX,ShiftedMatrixGeneralized<MATRIX,MASSMATRIX,VECTOR> > A;
00133 SmartPointer<const MASSMATRIX,ShiftedMatrixGeneralized<MATRIX,MASSMATRIX,VECTOR> > M;
00134
00138 VECTOR aux;
00139
00143 double sigma;
00144 };
00145
00146
00148
00149
00150 template <class MATRIX>
00151 inline
00152 ShiftedMatrix<MATRIX>::ShiftedMatrix (const MATRIX& A, const double sigma)
00153 :
00154 A(&A), sigma(sigma)
00155 {}
00156
00157
00158
00159 template <class MATRIX>
00160 inline void
00161 ShiftedMatrix<MATRIX>::shift (const double s)
00162 {
00163 sigma = s;
00164 }
00165
00166
00167 template <class MATRIX>
00168 inline double
00169 ShiftedMatrix<MATRIX>::shift () const
00170 {
00171 return sigma;
00172 }
00173
00174
00175
00176 template <class MATRIX>
00177 template <class VECTOR>
00178 inline void
00179 ShiftedMatrix<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
00180 {
00181 A->vmult(dst, src);
00182 if (sigma != 0.)
00183 dst.add(sigma, src);
00184 }
00185
00186
00187 template <class MATRIX>
00188 template <class VECTOR>
00189 inline double
00190 ShiftedMatrix<MATRIX>::residual (VECTOR& dst,
00191 const VECTOR& src,
00192 const VECTOR& rhs) const
00193 {
00194 A->vmult(dst, src);
00195 if (sigma != 0.)
00196 dst.add(sigma, src);
00197 dst.sadd(-1.,1.,rhs);
00198 return dst.l2_norm ();
00199 }
00200
00201
00202
00203 template <class MATRIX, class MASSMATRIX, class VECTOR>
00204 inline
00205 ShiftedMatrixGeneralized<MATRIX, MASSMATRIX, VECTOR>
00206 ::ShiftedMatrixGeneralized (const MATRIX& A,
00207 const MASSMATRIX& M,
00208 const double sigma)
00209 :
00210 A(&A), M(&M), sigma(sigma)
00211 {}
00212
00213
00214 template <class MATRIX, class MASSMATRIX, class VECTOR>
00215 inline void
00216 ShiftedMatrixGeneralized<MATRIX, MASSMATRIX, VECTOR>::shift (const double s)
00217 {
00218 sigma = s;
00219 }
00220
00221 template <class MATRIX, class MASSMATRIX, class VECTOR>
00222 inline double
00223 ShiftedMatrixGeneralized<MATRIX, MASSMATRIX, VECTOR>::shift () const
00224 {
00225 return sigma;
00226 }
00227
00228
00229 template <class MATRIX, class MASSMATRIX, class VECTOR>
00230 inline void
00231 ShiftedMatrixGeneralized<MATRIX, MASSMATRIX, VECTOR>::vmult (VECTOR& dst,
00232 const VECTOR& src) const
00233 {
00234 A->vmult(dst, src);
00235 if (sigma != 0.)
00236 {
00237 aux.reinit(dst);
00238 M.vmult(aux, src);
00239 dst.add(sigma, aux);
00240 }
00241 }
00242
00243
00244 template <class MATRIX, class MASSMATRIX, class VECTOR>
00245 inline double
00246 ShiftedMatrixGeneralized<MATRIX, MASSMATRIX, VECTOR>::residual (VECTOR& dst,
00247 const VECTOR& src,
00248 const VECTOR& rhs) const
00249 {
00250 A->vmult(dst, src);
00251 if (sigma != 0.)
00252 {
00253 aux.reinit(dst);
00254 M.vmult(aux, src);
00255 dst.add(sigma, aux);
00256 }
00257 dst.sadd(-1.,1.,rhs);
00258 return dst.l2_norm ();
00259 }
00260
00261
00262 DEAL_II_NAMESPACE_CLOSE
00263
00264 #endif
00265