207 STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
212 typedef typename matrix_type::connection connection;
213 m_L.resize_and_clear(A->num_rows(), A->num_cols());
214 m_U.resize_and_clear(A->num_rows(), A->num_cols());
217 std::vector<typename matrix_type::connection> con;
221 static const size_t velocity=0;
222 static const size_t pressure=1;
225 for(
typename matrix_type::row_iterator i_it = A->begin_row(0); i_it != A->end_row(0); ++i_it)
226 con.push_back(connection(i_it.index(), i_it.value()));
227 m_U.set_matrix_row(0, &con[0], con.size());
229 size_t totalentries=0;
232 for(
size_t i=1; i<A->num_rows(); i++)
235 if ((*A)(i,i)!=0) itype=velocity;
else itype=pressure;
244 for(
typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
246 con.push_back(connection(i_it.index(), i_it.value()));
251 typename matrix_type::row_iterator i_it = A->begin_row(i);
252 con.push_back(connection(i_it.index(), i_it.value()));
254 for(; i_it != A->end_row(i); ++i_it)
256 if (i_it.value()==0)
continue;
257 con.push_back(connection(i_it.index(), i_it.value()));
263 for(
size_t i_it = 0; i_it < con.size(); ++i_it)
265 size_t k = con[i_it].iIndex;
272 if(con[i_it].dValue == 0.0)
continue;
273 UG_ASSERT(
m_U.num_connections(k) != 0 &&
m_U.begin_row(k).index() == k,
"");
279 block_type &d = con[i_it].dValue = con[i_it].dValue / ukk;
281 typename matrix_type::row_iterator k_it =
m_U.begin_row(k);
284 while(k_it !=
m_U.end_row(k) && j < con.size())
288 if(k_it.index() == con[j].iIndex)
291 con[j].dValue -= k_it.value() * d;
294 else if(k_it.index() < con[j].iIndex)
299 typename matrix_type::connection
300 c(k_it.index(), k_it.value() * d * -1.0);
305 if ((*A)(k_it.index(),k_it.index())!=0) kitType=velocity;
else kitType=pressure;
307 if ((itype==velocity)&&(kitType==velocity)){
311 con.insert(con.begin()+j, c);
315 if ((itype==velocity)&&(kitType==pressure)){
319 con.insert(con.begin()+j, c);
323 if ((itype==pressure)&&(kitType==velocity)){
327 con.insert(con.begin()+j, c);
331 if ((itype==pressure)&&(kitType==pressure)){
335 con.insert(con.begin()+j, c);
350 if (k_it!=
m_U.end_row(k)){
351 for (;k_it!=
m_U.end_row(k);++k_it){
352 typename matrix_type::connection c(k_it.index(),-k_it.value()*d);
355 if ((*A)(k_it.index(),k_it.index())!=0) kitType=velocity;
else kitType=pressure;
356 if ((itype==velocity)&&(kitType==velocity)){
362 if ((itype==velocity)&&(kitType==pressure)){
368 if ((itype==pressure)&&(kitType==velocity)){
389 totalentries+=con.size();
390 if(maxentries < con.size()) maxentries = con.size();
393 m_L.set_matrix_row(i, &con[0], u_part);
394 m_U.set_matrix_row(i, &con[u_part], con.size()-u_part);
402 UG_LOG(
"\n CRILUT storage information:\n");
403 UG_LOG(
" A nr of connections: " << A->total_num_connections() <<
"\n");
404 UG_LOG(
" L+U nr of connections: " <<
m_L.total_num_connections()+
m_U.total_num_connections() <<
"\n");
405 UG_LOG(
" Increase factor: " << (
float)(
m_L.total_num_connections() +
m_U.total_num_connections() )/A->total_num_connections() <<
"\n");