Line data Source code
1 : /*
2 : * Copyright (c) 2013-2014: 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 "solve_deficit.h"
34 :
35 : #include "common/log.h"
36 : #include "common/debug_print.h"
37 : #include "lib_algebra/common/operations_vec.h"
38 :
39 : typedef int lapack_int;
40 : typedef float lapack_float;
41 : typedef double lapack_double;
42 : typedef int lapack_ftnlen;
43 :
44 :
45 : inline double dabs(double a) { return a > 0 ? a : -a; }
46 : /*
47 : extern "C"
48 : {
49 :
50 : // factor system *GETRF
51 : void dsytrs_(char *uplo, int *n, int *nrhs, double *a, int *lda,
52 : int *ipivot, double *b, int *ldb, int *info);
53 : }
54 :
55 :
56 : inline lapack_int sytrs(char uplo, int n, int nrhs, double *a, int lda,
57 : int *ipivot, double *b, int ldb)
58 : {
59 : lapack_int info;
60 : dsytrs_(&uplo, &n, &nrhs, a, &lda, ipivot, b, &ldb, &info);
61 : return info;
62 : }
63 :
64 : bool SolveDeficit2(DenseMatrix< VariableArray2<double> > &A,
65 : DenseVector<VariableArray1<double> > &x, DenseVector<VariableArray1<double> > &rhs)
66 : {
67 : DenseMatrix< VariableArray2<double> > A2=A;
68 : DenseVector<VariableArray1<double> > rhs2=rhs;
69 :
70 : std::vector<int> interchange(A.num_rows());
71 : int info = sytrs('U', A.num_rows(), 1, &A(0,0), A.num_rows(), &interchange[0], &rhs[0], A.num_rows());
72 : x=rhs;
73 : DenseVector<VariableArray1<double> > f;
74 : f = A2*x - rhs2;
75 : if(VecNormSquared(f) > 1e-12)
76 : {
77 : UG_LOG("solving was wrong:");
78 : UG_LOGN(CPPString(A2, "Aold"));
79 : rhs2.maple_print("rhs");
80 : UG_LOGN(CPPString(A, "Adecomp"));
81 : rhs.maple_print("rhsDecomp");
82 : x.maple_print("x");
83 : f.maple_print("f");
84 :
85 : }
86 :
87 : }*/
88 :
89 : namespace ug{
90 :
91 : /*
92 : | free variable
93 : xxxxxxxx
94 : 0xxxxxxx
95 : 000xxxxx
96 : 0000xxxx
97 : 00000xxx <- iNonNullRows = 5
98 : 00000000 \
99 : 00000000 | lin. dep.
100 : 00000000 /
101 :
102 : iNonNullRows = number of bounded variable (that is: number of linear independent rows)
103 : returns true if solveable.
104 :
105 : */
106 0 : bool Decomp(DenseMatrix< VariableArray2<double> > &A, DenseVector<VariableArray1<double> > &rhs,
107 : size_t &iNonNullRows, std::vector<size_t> &xinterchange, double deficitTolerance)
108 : {
109 : // LU Decomposition, IKJ Variant
110 : size_t rows = A.num_rows();
111 : size_t cols = A.num_cols();
112 :
113 0 : xinterchange.resize(cols);
114 0 : for(size_t k=0; k<cols; k++)
115 0 : xinterchange[k] = k;
116 :
117 : // UG_LOG("DECOMP " << rows << " x " << cols << "\n");
118 :
119 : std::vector<size_t> colInterchange;
120 : size_t k;
121 0 : for(k=0; k<cols; k++)
122 : {
123 : // UG_LOG("-----------------" << k << " < " << cols << " ------\n");
124 : size_t biggestR=k, biggestC=k;
125 0 : double dBiggest = dabs(A(k,k));
126 :
127 0 : for(size_t r=k; r<rows; r++)
128 0 : for(size_t c=k; c<cols; c++)
129 : {
130 0 : double n = dabs(A(r,c));
131 0 : if(dBiggest < n)
132 : {
133 : dBiggest=n;
134 : biggestR = r;
135 : biggestC = c;
136 : }
137 : }
138 :
139 :
140 : // UG_LOG(CPPString(A, "BeforeSwap"));
141 0 : if(dBiggest < 1e-14)
142 : {
143 : // UG_LOG("k = " << k << " abort");
144 : break;
145 : }
146 : // UG_LOG("k = " << k << " biggest = " << biggestR << ", " << biggestC << "\n");
147 0 : if(biggestR != k)
148 : {
149 0 : for(size_t j=0; j<cols; j++)
150 : std::swap(A(k, j), A(biggestR, j));
151 : std::swap(rhs[k], rhs[biggestR]);
152 : }
153 0 : if(biggestC != k)
154 : {
155 0 : for(size_t j=0; j<rows; j++)
156 : std::swap(A(j, k), A(j, biggestC));
157 : std::swap(xinterchange[k], xinterchange[biggestC]);
158 : }
159 : // UG_LOG(CPPString(A, "AfterSwap"));
160 :
161 0 : for(size_t i=k+1; i<rows; i++)
162 : {
163 0 : double a = A(i,k)/A(k,k);
164 0 : for(size_t j=k+1; j<cols; j++)
165 0 : A(i,j) = A(i,j) - a*A(k,j);
166 0 : rhs[i] -= a*rhs[k];
167 0 : A(i,k) = 0;
168 :
169 : }
170 :
171 :
172 : // UG_LOG(CPPString(A, "MAT"));
173 : // PrintVector(rhs, "rhs");
174 : }
175 0 : iNonNullRows = k;
176 : // every row below iNonNullRows is zero.
177 : // equation system is solveable if every rhs[k] = 0 for k>=iNonNullRows
178 0 : for(; k<rows; k++)
179 0 : if(dabs(rhs[k]) > deficitTolerance )
180 : {
181 : // UG_LOG("row " << k << " > 1e-7" << "\n")
182 : return false;
183 : }
184 : else
185 0 : rhs[k] = 0.0;
186 : // UG_LOGN("iNonNullRows = " << iNonNullRows << " abort");
187 : return true;
188 0 : }
189 :
190 :
191 0 : bool SolveDeficit(DenseMatrix< VariableArray2<double> > &A,
192 : DenseVector<VariableArray1<double> > &x, DenseVector<VariableArray1<double> > &rhs, double deficitTolerance)
193 : {
194 : DenseMatrix< VariableArray2<double> > A2=A;
195 : DenseVector<VariableArray1<double> > rhs2=rhs;
196 :
197 : UG_ASSERT(A.num_rows() == rhs.size(), "");
198 : UG_ASSERT(A.num_cols() == x.size(), "");
199 :
200 : size_t iNonNullRows;
201 0 : x.resize(A.num_cols());
202 0 : for(size_t i=0; i<x.size(); i++)
203 0 : x[i] = 0.0;
204 : std::vector<size_t> interchange;
205 0 : if(Decomp(A, rhs, iNonNullRows, interchange, deficitTolerance) == false) return false;
206 :
207 : // A.maple_print("Adecomp");
208 : // rhs.maple_print("rhs decomp");
209 :
210 0 : for(int i=iNonNullRows-1; i>=0; i--)
211 : {
212 0 : double d=A(i,i);
213 : double s=0;
214 0 : for(size_t k=i+1; k<A.num_cols(); k++)
215 0 : s += A(i,k)*x[interchange[k]];
216 0 : x[interchange[i]] = (rhs[i] - s)/d;
217 : }
218 : DenseVector<VariableArray1<double> > f;
219 0 : f = A2*x - rhs2;
220 0 : if(VecNormSquared(f) > 1e-2)
221 : {
222 0 : UG_LOGN("iNonNullRows = " << iNonNullRows);
223 : UG_LOG("solving was wrong:");
224 0 : UG_LOGN(CPPString(A2, "Aold"));
225 0 : rhs2.maple_print("rhs");
226 0 : UG_LOGN(CPPString(A, "Adecomp"));
227 0 : rhs.maple_print("rhsDecomp");
228 0 : x.maple_print("x");
229 0 : f.maple_print("f");
230 :
231 : }
232 :
233 : return true;
234 0 : }
235 : }
|
