00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__schur_matrix_h
00013 #define __deal2__schur_matrix_h
00014
00015 #include <deal.II/base/config.h>
00016 #include <deal.II/base/subscriptor.h>
00017 #include <deal.II/base/smartpointer.h>
00018 #include <deal.II/base/logstream.h>
00019 #include <deal.II/lac/vector_memory.h>
00020 #include <deal.II/lac/block_vector.h>
00021 #include <vector>
00022
00023 DEAL_II_NAMESPACE_OPEN
00024
00025
00092 template <class MA_inverse, class MB, class MDt, class MC>
00093 class SchurMatrix : public Subscriptor
00094 {
00095 public:
00096
00110 SchurMatrix(const MA_inverse& Ainv,
00111 const MB& B,
00112 const MDt& Dt,
00113 const MC& C,
00114 VectorMemory<BlockVector<double> >& mem,
00115 const std::vector<unsigned int>& signature = std::vector<unsigned int>(0));
00116
00135 void prepare_rhs (BlockVector<double>& dst,
00136 const BlockVector<double>& src) const;
00137
00142 void vmult (BlockVector<double>& dst,
00143 const BlockVector<double>& src) const;
00144
00145
00146
00151 double residual (BlockVector<double>& dst,
00152 const BlockVector<double>& src,
00153 const BlockVector<double>& rhs) const;
00154
00160 void postprocess (BlockVector<double>& dst,
00161 const BlockVector<double>& src,
00162 const BlockVector<double>& rhs) const;
00163
00172 void debug_level(unsigned int l);
00173 private:
00177 SchurMatrix (const SchurMatrix<MA_inverse, MB, MDt, MC>&);
00181 SchurMatrix& operator = (const SchurMatrix<MA_inverse, MB, MDt, MC>&);
00182
00186 const SmartPointer<const MA_inverse,SchurMatrix<MA_inverse,MB,MDt,MC> > Ainv;
00190 const SmartPointer<const MB,SchurMatrix<MA_inverse,MB,MDt,MC> > B;
00194 const SmartPointer<const MDt,SchurMatrix<MA_inverse,MB,MDt,MC> > Dt;
00198 const SmartPointer<const MC,SchurMatrix<MA_inverse,MB,MDt,MC> > C;
00202 VectorMemory<BlockVector<double> >& mem;
00203
00207 std::vector<unsigned int> signature;
00208
00212 unsigned int debug;
00213 };
00214
00216
00217
00218 template <class MA_inverse, class MB, class MDt, class MC>
00219 SchurMatrix<MA_inverse, MB, MDt, MC>
00220 ::SchurMatrix(const MA_inverse& Ainv,
00221 const MB& B,
00222 const MDt& Dt,
00223 const MC& C,
00224 VectorMemory<BlockVector<double> >& mem,
00225 const std::vector<unsigned int>& signature)
00226 : Ainv(&Ainv), B(&B), Dt(&Dt), C(&C),
00227 mem(mem),
00228 signature(signature),
00229 debug(0)
00230 {
00231 }
00232
00233
00234 template <class MA_inverse, class MB, class MDt, class MC>
00235 void
00236 SchurMatrix<MA_inverse, MB, MDt, MC>
00237 ::debug_level(unsigned int l)
00238 {
00239 debug = l;
00240 }
00241
00242
00243 template <class MA_inverse, class MB, class MDt, class MC>
00244 void SchurMatrix<MA_inverse, MB, MDt, MC>
00245 ::vmult(BlockVector<double>& dst,
00246 const BlockVector<double>& src) const
00247 {
00248 deallog.push("Schur");
00249 if (debug > 0)
00250 deallog << "src:" << src.l2_norm() << std::endl;
00251
00252 C->vmult(dst, src);
00253 if (debug > 0)
00254 deallog << "C:" << dst.l2_norm() << std::endl;
00255
00256 BlockVector<double>* h1 = mem.alloc();
00257 if (signature.size()>0)
00258 h1->reinit(signature);
00259 else
00260 h1->reinit(B->n_block_cols(), src.block(0).size());
00261 Dt->Tvmult(*h1,src);
00262 if (debug > 0)
00263 deallog << "Dt:" << h1->l2_norm() << std::endl;
00264
00265 BlockVector<double>* h2 = mem.alloc();
00266 h2->reinit(*h1);
00267 Ainv->vmult(*h2, *h1);
00268 if (debug > 0)
00269 deallog << "Ainverse:" << h2->l2_norm() << std::endl;
00270
00271 mem.free(h1);
00272 B->vmult_add(dst, *h2);
00273 if (debug > 0)
00274 deallog << "dst:" << dst.l2_norm() << std::endl;
00275
00276 mem.free(h2);
00277 deallog.pop();
00278 }
00279
00280
00281 template <class MA_inverse, class MB, class MDt, class MC>
00282 double SchurMatrix<MA_inverse, MB, MDt, MC>
00283 ::residual(BlockVector<double>& dst,
00284 const BlockVector<double>& src,
00285 const BlockVector<double>& rhs) const
00286 {
00287 vmult(dst, src);
00288 dst.scale(-1.);
00289 dst += rhs;
00290 return dst.l2_norm();
00291 }
00292
00293
00294 template <class MA_inverse, class MB, class MDt, class MC>
00295 void SchurMatrix<MA_inverse, MB, MDt, MC>
00296 ::prepare_rhs(BlockVector<double>& dst,
00297 const BlockVector<double>& src) const
00298 {
00299 Assert (src.n_blocks() == B->n_block_cols(),
00300 ExcDimensionMismatch(src.n_blocks(), B->n_block_cols()));
00301 Assert (dst.n_blocks() == B->n_block_rows(),
00302 ExcDimensionMismatch(dst.n_blocks(), B->n_block_rows()));
00303
00304 deallog.push("Schur-prepare");
00305 if (debug > 0)
00306 deallog << "src:" << src.l2_norm() << std::endl;
00307 BlockVector<double>* h1 = mem.alloc();
00308 if (signature.size()>0)
00309 h1->reinit(signature);
00310 else
00311 h1->reinit(B->n_block_cols(), src.block(0).size());
00312 Ainv->vmult(*h1, src);
00313 if (debug > 0)
00314 deallog << "Ainverse:" << h1->l2_norm() << std::endl;
00315 B->vmult_add(dst, *h1);
00316 if (debug > 0)
00317 deallog << "dst:" << dst.l2_norm() << std::endl;
00318 mem.free(h1);
00319 deallog.pop();
00320 }
00321
00322
00323 template <class MA_inverse, class MB, class MDt, class MC>
00324 void SchurMatrix<MA_inverse, MB, MDt, MC>
00325 ::postprocess(BlockVector<double>& dst,
00326 const BlockVector<double>& src,
00327 const BlockVector<double>& rhs) const
00328 {
00329 Assert (dst.n_blocks() == B->n_block_cols(),
00330 ExcDimensionMismatch(dst.n_blocks(), B->n_block_cols()));
00331 Assert (rhs.n_blocks() == B->n_block_cols(),
00332 ExcDimensionMismatch(rhs.n_blocks(), B->n_block_cols()));
00333 Assert (src.n_blocks() == B->n_block_rows(),
00334 ExcDimensionMismatch(src.n_blocks(), B->n_block_rows()));
00335
00336 deallog.push("Schur-post");
00337 if (debug > 0)
00338 deallog << "src:" << src.l2_norm() << std::endl;
00339 BlockVector<double>* h1 = mem.alloc();
00340 if (signature.size()>0)
00341 h1->reinit(signature);
00342 else
00343 h1->reinit(B->n_block_cols(), src.block(0).size());
00344 Dt->Tvmult(*h1, src);
00345 if (debug > 0)
00346 deallog << "Dt:" << h1->l2_norm() << std::endl;
00347 h1->sadd(-1.,rhs);
00348 Ainv->vmult(dst,*h1);
00349 if (debug > 0)
00350 deallog << "dst:" << dst.l2_norm() << std::endl;
00351 mem.free(h1);
00352 deallog.pop();
00353 }
00354
00355
00356 DEAL_II_NAMESPACE_CLOSE
00357
00358 #endif
00359