00001
00002
00003
00004 #ifndef compositions_hh
00005 #define compositions_hh
00006
00007 #include "basics/outputOperator.hh"
00008 #include "basics/exceptions.hh"
00009 #include "space/space.hh"
00010 #include "function/basis.hh"
00011 #include "function/vector.hh"
00012 #include "basics/debug.hh"
00013
00014
00015 #define ComposeConstr_D 0
00016 #define ComposeAppl_D 0
00017 #define SchurComplContr_D 0
00018
00019 namespace concepts {
00020
00021
00022 template<class F>
00023 class SparseMatrix;
00024
00025 template<class F>
00026 class DenseMatrix;
00027
00028 template<class F>
00029 class Matrix;
00030
00031 template<class F>
00032 class HashedSparseMatrix;
00033
00034
00035
00039 template<class F>
00040 class Operator : public OutputOperator {
00041 public:
00043 typedef F type;
00045 typedef typename Realtype<F>::type r_type;
00047 typedef typename Cmplxtype<F>::type c_type;
00048
00049 Operator(uint dimX, uint dimY) : dimX_(dimX), dimY_(dimY) {}
00050 virtual ~Operator() {}
00051
00067 virtual void operator()(const Function<r_type>& fncY, Function<F>& fncX);
00082 virtual void operator()(const Function<c_type>& fncY, Function<c_type>& fncX);
00083
00087 virtual inline const uint dimX() const { return dimX_; }
00088
00092 virtual inline const uint dimY() const { return dimY_; }
00093 protected:
00094 virtual std::ostream& info(std::ostream& os) const;
00096 uint dimX_, dimY_;
00097 };
00098
00099
00100
00106 template<class F>
00107 class VecOperator : public Operator<F> {
00108 public:
00110 typedef typename Realtype<F>::type r_type;
00112 typedef typename Cmplxtype<F>::type c_type;
00113
00114 VecOperator(uint dimX, uint dimY) : Operator<F>(dimY, dimX) {}
00115
00116 virtual void operator()(const Function<r_type>& fncY, Function<F>& fncX);
00117 virtual void operator()(const Function<c_type>& fncY,
00118 Function<c_type>& fncX);
00119
00135 virtual void operator()(const Vector<r_type>& fncY, Vector<F>& fncX);
00150 virtual void operator()(const Vector<c_type>& fncY, Vector<c_type>& fncX);
00151
00153 void operator()(const Matrix<r_type>& mX, Matrix<F>& mY);
00155 void operator()(const Matrix<c_type>& mX, Matrix<c_type>& mY);
00156 protected:
00157 virtual std::ostream& info(std::ostream& os) const;
00161 virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX) = 0;
00162 };
00163
00164
00165
00176 template<class F, class H = F>
00177 class Compose : public Operator<F> {
00178 public:
00180 inline Compose(Operator<F>& A, Operator<H>& B);
00181
00184 inline void operator()(const Function<F>& fncY, Function<F>& fncX);
00185
00194 bool collapse(Matrix<F>& dest) const;
00195 protected:
00196 virtual std::ostream& info(std::ostream& os) const;
00197 private:
00199 Operator<F>& A_;
00201 Operator<H>& B_;
00203 Vector<F> f_;
00204
00205 void collapse_(DenseMatrix<F>& A, DenseMatrix<H>& B, Matrix<F>& dest) const;
00206 };
00207
00208 template<class F, class H>
00209 Compose<F,H>::Compose(Operator<F>& A, Operator<H>& B)
00210 : Operator<F>(A.dimX(), B.dimY()), A_(A), B_(B), f_(B.dimX()) {
00211 DEBUGL(ComposeConstr_D,
00212 "A is (" << A.dimX() << "x" << A.dimY() << ") *" <<
00213 "B is (" << B.dimX() << "x" << B.dimY() << ")");
00214 conceptsAssert(A.dimY() == B.dimX(), Assertion());
00215 }
00216
00217 template<class F, class H>
00218 void Compose<F,H>::operator()(const Function<F>& fncY, Function<F>& fncX)
00219 {
00220 DEBUGL(ComposeAppl_D, *this);
00221 conceptsAssert3(fncY.dim() == this->dimY_, Assertion(),
00222 "fncY.dim = " << fncY.dim() << ", dimY = " << this->dimY());
00223 conceptsAssert3(fncX.dim() == this->dimX(), Assertion(),
00224 "fncX.dim = " << fncX.dim() << ", dimX = " << this->dimX());
00225 DEBUGL(ComposeAppl_D, "y = " << fncY);
00226 B_(fncY, f_);
00227 DEBUGL(ComposeAppl_D, "f = " << f_);
00228 A_(f_, fncX);
00229 DEBUGL(ComposeAppl_D, "x = " << fncX);
00230 }
00231
00232
00233
00237 template<class F>
00238 class Multiple : public Operator<F> {
00239 public:
00241 typedef typename Realtype<F>::type r_type;
00243 typedef typename Cmplxtype<F>::type c_type;
00244
00245 Multiple(Operator<F>& A, F a)
00246 : Operator<F>(A.dimX(), A.dimY()), A_(A), a_(a) {}
00247
00255 inline void operator()(const Function<r_type>& fncY, Function<F>& fncX)
00256 {
00257 apply_(fncY, fncX);
00258 }
00259 inline void operator()(const Function<c_type>& fncY,
00260 Function<c_type>& fncX) {
00261 apply_(fncY, fncX);
00262 }
00263
00264 protected:
00265 virtual std::ostream& info(std::ostream& os) const;
00266 private:
00268 Operator<F>& A_;
00270 F a_;
00271
00272 template<class H, class I>
00273 inline void apply_(const Function<H>& fncY, Function<I>& fncX);
00274 };
00275
00276 template<class F>
00277 template<class H, class I>
00278 void Multiple<F>:: apply_(const Function<H>& fncY, Function<I>& fncX) {
00279 conceptsAssert(fncY.dim() == this->dimY(), Assertion());
00280 conceptsAssert(fncX.dim() == this->dimX(), Assertion());
00281 A_(fncY, fncX);
00282 if (a_ != 1.0)
00283 fncX *= a_;
00284 }
00285
00286
00287
00292 template<class F>
00293 class LiCoI : public Operator<F> {
00294 public:
00296 typedef typename Realtype<F>::type r_type;
00298 typedef typename Cmplxtype<F>::type c_type;
00299
00300 LiCoI(Operator<F>& A, F a = 1.0, F b = 1.0)
00301 : Operator<F>(A.dimX(), A.dimY()), A_(A), a_(a), b_(b) {}
00302
00310 virtual void operator()(const Function<r_type>& fncY, Function<F>& fncX);
00311 virtual void operator()(const Function<c_type>& fncY,
00312 Function<c_type>& fncX);
00313
00330 bool collapse(Matrix<F>& dest, const F fact = 1.0) const;
00331 protected:
00332 virtual std::ostream& info(std::ostream& os) const;
00333 private:
00335 Operator<F>& A_;
00337 F a_;
00339 F b_;
00340 };
00341
00342
00343
00350 template<class F>
00351 class LiCo : public Operator<F> {
00352 public:
00354 typedef typename Realtype<F>::type r_type;
00356 typedef typename Cmplxtype<F>::type c_type;
00357
00366 LiCo(Operator<F>& A, Operator<F>& B, F a = 1.0, F b = 1.0) :
00367 Operator<F>(A.dimX(), A.dimY()), A_(A), B_(B), a_(a), b_(b) {
00368 conceptsAssert(A.dimX() == B.dimX(), Assertion());
00369 conceptsAssert(A.dimY() == B.dimY(), Assertion());
00370 }
00371
00380 inline void operator()(const Function<r_type>& fncY, Function<F>& fncX)
00381 {
00382 apply_(fncY, fncX);
00383 }
00384 inline void operator()(const Function<c_type>& fncY, Function<c_type>& fncX)
00385 {
00386 apply_(fncY, fncX);
00387 }
00388
00404 bool collapse(Matrix<F>& dest, const F fact = 1.0) const;
00405 protected:
00406 virtual std::ostream& info(std::ostream& os) const;
00407 private:
00409 Operator<F>& A_;
00411 Operator<F>& B_;
00413 F a_;
00415 F b_;
00416
00417 template<class H, class I>
00418 inline void apply_(const Function<H>& fncY, Function<I>& fncX);
00419 };
00420
00421 template<class F>
00422 template<class H, class I>
00423 void LiCo<F>::apply_(const Function<H>& fncY, Function<I>& fncX) {
00424 conceptsAssert(this->dimX() == fncX.dim(), Assertion());
00425 conceptsAssert(this->dimY() == fncY.dim(), Assertion());
00426 A_(fncY, fncX);
00427 fncX *= a_;
00428 Vector<I> f(this->dimX());
00429 B_(fncY, f);
00430 fncX.add(f, b_);
00431 }
00432
00433
00434
00443 template<class F>
00444 class SchurCompl : public LiCo<F> {
00445 public:
00447 SchurCompl(Operator<F>& A_II_inv, Operator<F>& A_IB,
00448 Operator<F>& A_BI, Operator<F>& A_BB);
00449 virtual ~SchurCompl() {
00450 delete D_; delete C_;
00451 }
00452 private:
00454 Compose<F> *C_, *D_;
00455 };
00456
00457 template<class F>
00458 SchurCompl<F>::SchurCompl(Operator<F>& A_II_inv, Operator<F>& A_IB,
00459 Operator<F>& A_BI, Operator<F>& A_BB)
00460 : LiCo<F>(A_BB,
00461 *(D_ = new Compose<F>(A_BI,
00462 *(C_ = new Compose<F>(A_II_inv, A_IB)))
00463 ), 1.0, -1.0)
00464 {
00465 DEBUGL(SchurComplContr_D, "A_II_inv = " <<
00466 (DenseMatrix<F>(SparseMatrix<F>(A_II_inv, true))));
00467 DEBUGL(SchurComplContr_D, "A_IB = " <<
00468 (DenseMatrix<F>(SparseMatrix<F>(A_IB, true))));
00469 DEBUGL(SchurComplContr_D, "A_BI = " <<
00470 (DenseMatrix<F>(SparseMatrix<F>(A_BI, true))));
00471 DEBUGL(SchurComplContr_D, "A_BB = " <<
00472 (DenseMatrix<F>(SparseMatrix<F>(A_BB, true))));
00473 DEBUGL(SchurComplContr_D, "C_ = AII^-1 * AIB = " <<
00474 (DenseMatrix<F>(SparseMatrix<F>(*C_, true))));
00475 DEBUGL(SchurComplContr_D, "D_ = ABI * AII^-1 * AIB = " <<
00476 (DenseMatrix<F>(SparseMatrix<F>(*D_, true))));
00477 DEBUGL(SchurComplContr_D, "S_ = A_BB - A_BI * A_II^-1 * A_IB = " <<
00478 (DenseMatrix<F>(SparseMatrix<F>(*this, true))));
00479 }
00480
00481 }
00482
00483 #endif // compositions_hh