Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__LIB_ALGEBRA__LAPACK_ANALYZING_SOLVER__
34 : #define __H__LIB_ALGEBRA__LAPACK_ANALYZING_SOLVER__
35 :
36 : #include "../interface/linear_operator_inverse.h"
37 : #include "../interface/matrix_operator.h"
38 : #include "common/error.h"
39 : #include "common/util/smart_pointer.h"
40 : #include "common/util/histogramm.h"
41 :
42 : namespace ug{
43 :
44 : void checksub(const CPUAlgebra::matrix_type &A);
45 :
46 : template <typename M, typename X, typename Y = X>
47 : class AnalyzingSolver
48 : : public virtual ILinearOperatorInverse<X,Y>
49 : {
50 : public:
51 : /// Domain space
52 : typedef X domain_function_type;
53 :
54 : /// Range space
55 : typedef Y codomain_function_type;
56 :
57 : /// Matrix type
58 : typedef M matrix_type;
59 :
60 : public:
61 0 : virtual bool apply(Y& u, const X& f)
62 : {
63 0 : return m_pLinearOperatorInverse->apply(u, f);
64 : }
65 :
66 0 : virtual bool apply_return_defect(Y& u, X& f)
67 : {
68 0 : return m_pLinearOperatorInverse->apply_return_defect(u, f);
69 : }
70 :
71 0 : AnalyzingSolver(SmartPtr<ILinearOperatorInverse<X,Y> > pLinearOperatorInverse)
72 0 : {
73 0 : m_pLinearOperatorInverse = pLinearOperatorInverse;
74 : }
75 :
76 : /// virtual destructor
77 0 : virtual ~AnalyzingSolver() {};
78 :
79 : /// returns if parallel solving is supported
80 0 : virtual bool supports_parallel() const
81 : {
82 0 : return m_pLinearOperatorInverse->supports_parallel();
83 : }
84 :
85 : public:
86 :
87 :
88 :
89 :
90 0 : void check(const matrix_type &A)
91 : {
92 : UG_LOG("ANALYZING SOLVER:\n");
93 : UG_LOG(" Matrix is of dimension " << A.num_rows() << "\n");
94 0 : if(A.num_rows() != A.num_cols())
95 : { UG_LOG(" Matrix is not quadratic???\n"); }
96 0 : if(GetRows(A(0,0)) == 1)
97 : { UG_LOG(" Submatrices are DOUBLE.\n")}
98 : else
99 0 : {UG_LOG(" Submatrices are Matrices of size " << GetRows(A(0,0)) << " x " << GetCols(A(0,0)) << "\n");}
100 :
101 :
102 : ////
103 : // check symmetry
104 :
105 :
106 : /////
107 :
108 :
109 : const size_t nrOfRows = block_traits<typename matrix_type::value_type>::static_num_rows;
110 0 : size_t m_size = A.num_rows() * nrOfRows;
111 :
112 0 : CPUAlgebra::matrix_type mat;
113 0 : mat.resize_and_clear(m_size, m_size);
114 :
115 0 : for(size_t r=0; r<A.num_rows(); r++)
116 0 : for(typename matrix_type::const_row_iterator it = A.begin_row(r); it != A.end_row(r); ++it)
117 : {
118 0 : size_t rr = r*nrOfRows;
119 0 : size_t cc = it.index()*nrOfRows;
120 0 : for(size_t r2=0; r2<nrOfRows; r2++)
121 0 : for(size_t c2=0; c2<nrOfRows; c2++)
122 : {
123 0 : if(BlockRef(it.value(), r2, c2) != 0.0)
124 0 : mat(rr + r2, cc + c2) = BlockRef(it.value(), r2, c2);
125 : }
126 : }
127 0 : mat.defragment();
128 0 : checksub(mat);
129 0 : }
130 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > A, const Y&u)
131 : {
132 0 : check(A);
133 0 : return m_pLinearOperatorInverse->init(A, u);
134 : }
135 :
136 0 : virtual bool init(SmartPtr<ILinearOperator<Y,X> > A)
137 : {
138 0 : check(A);
139 0 : return m_pLinearOperatorInverse->init(A);
140 : }
141 :
142 0 : void check(SmartPtr<ILinearOperator<Y,X> > A)
143 : {
144 : // cast operator
145 0 : SmartPtr<MatrixOperator<M,Y,X> > op =
146 : A.template cast_dynamic<MatrixOperator<M,Y,X> >();
147 :
148 : // check if correct types are present
149 0 : if(op.invalid())
150 0 : UG_THROW("IMatrixOperatorInverse::init:"
151 : " Passed operator is not matrix-based.");
152 :
153 : // forward request
154 0 : check(*op);
155 0 : }
156 0 : virtual const char *name() const { return m_pLinearOperatorInverse->name(); }
157 :
158 0 : virtual std::string config_string() const
159 : {
160 0 : return "AnalyzingSolver " + m_pLinearOperatorInverse->config_string();
161 : }
162 : private:
163 : SmartPtr<ILinearOperatorInverse<X,Y> > m_pLinearOperatorInverse;
164 : };
165 :
166 :
167 : } // end namespace ug
168 :
169 : #endif /* __H__LIB_ALGEBRA__LAPACK_LU_OPERATOR__ */
|
