114 mat.set_as_copy_of(*pOp);
115 STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
119 typedef typename matrix_type::connection connection;
120 m_L.resize_and_clear(A->num_rows(), A->num_cols());
121 m_U.resize_and_clear(A->num_rows(), A->num_cols());
126 std::vector<typename matrix_type::connection> con;
131 for(
typename matrix_type::row_iterator i_it = A->begin_row(0); i_it != A->end_row(0); ++i_it)
132 con.push_back(connection(i_it.index(), i_it.value()));
133 m_U.set_matrix_row(0, &con[0], con.size());
135 size_t totalentries=0;
140 "Using ILUT(" <<
m_eps <<
") on " << A->num_rows() <<
" x " << A->num_rows() <<
" matrix...");
143 for(
size_t i=1; i<A->num_rows()/2; i++)
151 for(
typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
153 con.push_back(connection(i_it.index(), i_it.value()));
160 totalentries+=con.size();
161 if(maxentries < con.size()) maxentries = con.size();
163 for(
size_t i_it=0; i_it<con.size(); i_it++)
165 if(con[i_it].iIndex < i)
m_L(i, con[i_it].iIndex) = con[i_it].dValue;
166 else m_U(i, con[i_it].iIndex) = con[i_it].dValue;
171 for(
size_t i=A->num_rows()/2; i<A->num_rows(); i++)
178 for(
typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
180 con.push_back(connection(i_it.index(), i_it.value()));
187 totalentries+=con.size();
188 if(maxentries < con.size()) maxentries = con.size();
190 A->set_matrix_row(i, &con[0], con.size());
194 for(
typename matrix_type::row_iterator i_it = A->begin_row(A->num_rows()/2); i_it != A->end_row(A->num_rows()/2); ++i_it)
196 if(i_it.index() <= A->num_rows()/2)
198 con.push_back(connection(i_it.index(), i_it.value()));
200 m_L.set_matrix_row(A->num_rows()/2, &con[0], u_part);
201 m_U.set_matrix_row(A->num_rows()/2, &con[u_part], con.size()-u_part);
204 for(
size_t i=A->num_rows()/2+1; i<A->num_rows(); i++)
212 for(
typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
214 con.push_back(connection(i_it.index(), i_it.value()));
221 totalentries+=con.size();
222 if(maxentries < con.size()) maxentries = con.size();
224 for(
size_t i_it=0; i_it<con.size(); i_it++)
226 if(con[i_it].iIndex < i)
m_L(i, con[i_it].iIndex) = con[i_it].dValue;
227 else m_U(i, con[i_it].iIndex) = con[i_it].dValue;
240 UG_LOG(
"\n ILUT storage information:\n");
241 UG_LOG(
" A nr of connections: " << A->total_num_connections() <<
"\n");
242 UG_LOG(
" L+U nr of connections: " <<
m_L.total_num_connections()+
m_U.total_num_connections() <<
"\n");
243 UG_LOG(
" Increase factor: " << (
float)(
m_L.total_num_connections() +
m_U.total_num_connections() )/A->total_num_connections() <<
"\n");
244 UG_LOG(
reset_floats <<
"Total entries: " << totalentries <<
" (" << ((
double)totalentries) / (A->num_rows()*A->num_rows()) <<
"% of dense)");
250 size_t eliminate_row(std::vector<typename matrix_type::connection> &con,
size_t start,
size_t stop,
double dmax)
254 for(
size_t i_it = 0; i_it < con.size(); ++i_it)
256 size_t k = con[i_it].iIndex;
257 if(k < start)
continue;
264 if(con[i_it].dValue == 0.0)
continue;
265 UG_ASSERT(
m_U.num_connections(k) != 0 &&
m_U.begin_row(k).index() == k,
"");
271 block_type &d = con[i_it].dValue = con[i_it].dValue / ukk;
273 typename matrix_type::row_iterator k_it =
m_U.begin_row(k);
276 while(k_it !=
m_U.end_row(k) && j < con.size())
279 if(k_it.index() == con[j].iIndex)
282 con[j].dValue -= k_it.value() * d;
285 else if(k_it.index() < con[j].iIndex)
290 typename matrix_type::connection
291 c(k_it.index(), k_it.value() * d * -1.0);
296 con.insert(con.begin()+j, c);
309 if (k_it!=
m_U.end_row(k)){
310 for (;k_it!=
m_U.end_row(k);++k_it){
311 typename matrix_type::connection c(k_it.index(),-k_it.value()*d);