Line data Source code
1 : /*
2 : * Copyright (c) 2014-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 : #include "lib_algebra/lib_algebra.h"
34 : #include "lib_algebra/cpu_algebra_types.h"
35 : #include "common/util/histogramm.h"
36 : namespace ug{
37 0 : void checksub(const CPUAlgebra::matrix_type &A)
38 : {
39 0 : if(A.num_rows() == 0)
40 : {
41 : UG_LOG("EMPTY MATRIX!\n");
42 0 : return;
43 : }
44 0 : A.print("A");
45 : UG_LOG(reset_floats);
46 : size_t N = A.num_rows();
47 : // check isolated
48 : size_t iIsolated = 0;
49 0 : for(size_t r=0; r<A.num_rows(); r++)
50 0 : if(A.is_isolated(r)) iIsolated++;
51 :
52 : typedef CPUAlgebra::matrix_type::const_row_iterator row_it;
53 : // typedef CPUAlgebra::matrix_type::value_type value_type;
54 :
55 0 : UG_LOG("Nr of dirichlet nodes: " << iIsolated << " (" << iIsolated*100.0/N << "% )\n");
56 :
57 : // check symmetric
58 :
59 0 : std::vector<double> alpha(N, 0);
60 :
61 : size_t iUnsymmetric=0;
62 : const double unsymmetricEps = 1e-8;
63 0 : for(size_t r=0; r<A.num_rows(); r++)
64 : {
65 : double dUnsymmetric = 0;
66 0 : if(A.is_isolated(r) ) continue;
67 :
68 0 : for(row_it it = A.begin_row(r); it != A.end_row(r); ++it)
69 : {
70 : size_t c = it.index();
71 0 : if(A.is_isolated(c) ) continue;
72 0 : double T = A(c, r);
73 0 : if(A(c, r) != MatrixTranspose(it.value()))
74 : {
75 0 : dUnsymmetric += dabs(T-it.value());
76 : }
77 :
78 : }
79 0 : double diag = A(r, r);
80 0 : if(diag == 0.0) continue;
81 :
82 0 : dUnsymmetric /= diag;
83 :
84 0 : alpha[r] = dUnsymmetric;
85 0 : if(dUnsymmetric > unsymmetricEps)
86 0 : iUnsymmetric++;
87 :
88 : }
89 0 : std::sort(alpha.begin(), alpha.end());
90 :
91 : // check sign condition
92 :
93 0 : if(iUnsymmetric==0)
94 0 : { UG_LOG("Matrix is symmetric! (maximum alpha = " << alpha[N-1] << ")\n"); }
95 : else
96 : {
97 0 : UG_LOG("Matrix is unsymmetric in " << (iUnsymmetric*100)/(N-iIsolated) << "% of the non-dirchlet rows (" << iUnsymmetric << " total)\n");
98 : UG_LOG(" row i alpha-unsymmetric means: sum_{A_{ij} != 0} |A_{ij}-A{ji}| / A_{ii} >= alpha\n")
99 :
100 : UG_LOG("> alpha distribution:\n");
101 0 : UG_LOG(DistributionPercentage(alpha));
102 : }
103 : size_t signConditionMet=0;
104 : size_t zeroDiagonal=0;
105 :
106 0 : double minEW = 1e20;
107 0 : double maxEW = 1e-20;
108 :
109 0 : for(size_t r=0; r<A.num_rows(); r++)
110 : {
111 0 : if(A(r, r)==0.0)
112 : {
113 : zeroDiagonal++;
114 0 : continue;
115 : }
116 : bool bPos = A(r, r) > 0;
117 : double s=0.0;
118 : bool bSignCondMet = true;
119 0 : for(row_it it = A.begin_row(r); it != A.end_row(r); ++it)
120 : {
121 0 : if(it.index() == r) continue;
122 0 : if(bPos)
123 : {
124 0 : if(it.value() > 0)
125 : bSignCondMet = false;
126 : }
127 : else
128 : {
129 0 : if(it.value() < 0)
130 : bSignCondMet = false;
131 : }
132 0 : s += dabs(it.value());
133 : }
134 0 : if(bSignCondMet && A.is_isolated(r) == false)
135 0 : signConditionMet++;
136 :
137 0 : minEW = std::min(minEW, A(r,r)-s);
138 0 : maxEW = std::max(maxEW, A(r,r)+s);
139 : }
140 :
141 0 : if(signConditionMet == N-iIsolated)
142 : { UG_LOG("Sign condition met in all nodes\n"); }
143 : else
144 : {
145 0 : UG_LOG("Sign condition met in " << (signConditionMet*100.0)/(N-iIsolated) << "% of the non-dirichlet rows (" << signConditionMet << " total)\n");
146 : }
147 0 : UG_LOG("Gershgorin Eigenvalues are within [" << minEW << ", " << maxEW << "]\n");
148 :
149 0 : }
150 : }
|
