From fa515268ee5ac1a7dbdc7f2ff7ae8bb351690353 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 22 Mar 2010 10:15:34 +0000 Subject: [PATCH] Better comments for local_to_global functions. Commented GlobalRowsFromLocal struct. git-svn-id: https://svn.dealii.org/trunk@20874 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/lac/constraint_matrix.templates.h | 437 ++++++++++++------ 1 file changed, 283 insertions(+), 154 deletions(-) diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index d9a2516ad3..62c5937314 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -742,20 +742,26 @@ distribute_local_to_global (const Vector &local_vector, else for (unsigned int i=0; iinhomogeneity; if (val != 0) for (unsigned int j=0; j + // > >, but tuned so that frequent memory + // allocation for each entry is + // avoided. This is not directly used in + // the user code, but accessed via the + // GlobalRowsFromLocal. struct DataCache { DataCache () @@ -1096,77 +1112,162 @@ namespace internals - // collects all the global rows and - // their origin (direct/constraint) - // basically a vector of distributing - // and the data cache. with some - // specialized sort and insert functions. + // collects all the global rows from a + // local contribution (cell) and their + // origin (direct/constraint). this is + // basically a vector consisting of + // "Distributing" structs using access via + // the DataCache. Provides some + // specialized sort and insert functions. + // + // in case there are no constraints, this is + // basically a list of pairs with + // the first index being the global index and + // the second index the local index. The list + // is sorted with respect to the global index. + // + // in case there are constraints, a global dof + // might get a contribution also because it + // gets data from a constrained dof. This + // means that a global dof might also have + // indirect contributions from a local dof via + // a constraint, besides the direct ones. struct GlobalRowsFromLocal { GlobalRowsFromLocal (const unsigned int n_local_rows) : total_row_indices (n_local_rows), n_active_rows (n_local_rows), n_inhomogeneous_rows (0){}; + + // implemented below void insert_index (const unsigned int global_row, const unsigned int local_row, const double constraint_value); void sort (); - const unsigned int & n_additional_dofs (const unsigned int local_dofs); + + // return all kind of information on the + // constraints + + // returns the number of global indices in the + // struct unsigned int size () const { return n_active_rows; }; - unsigned int & global_row (const unsigned int loc_index) - { return total_row_indices[loc_index].global_row; }; - unsigned int size (const unsigned int loc_index) const - { return (total_row_indices[loc_index].constraint_position == + + // returns the global index of the + // counter_index-th entry in the list + unsigned int & global_row (const unsigned int counter_index) + { return total_row_indices[counter_index].global_row; }; + + // returns the number of constraints that are + // associated to the counter_index-th entry in + // the list + unsigned int size (const unsigned int counter_index) const + { return (total_row_indices[counter_index].constraint_position == numbers::invalid_unsigned_int ? 0 : - data_cache.get_size(total_row_indices[loc_index].constraint_position)); }; - const unsigned int & global_row (const unsigned int loc_index) const - { return total_row_indices[loc_index].global_row; }; - const unsigned int & local_row (const unsigned int loc_index) const - { return total_row_indices[loc_index].local_row; }; - unsigned int & local_row (const unsigned int loc_index) - { return total_row_indices[loc_index].local_row; }; - unsigned int local_row (const unsigned int loc_index, + data_cache.get_size(total_row_indices[counter_index]. + constraint_position)); }; + + // returns the global row associated with the + // counter_index-th entry in the list + const unsigned int & global_row (const unsigned int counter_index) const + { return total_row_indices[counter_index].global_row; }; + + // returns the local row in the cell matrix + // associated with the counter_index-th entry + // in the list. Returns invalid_unsigned_int + // for invalid unsigned ints + const unsigned int & local_row (const unsigned int counter_index) const + { return total_row_indices[counter_index].local_row; }; + + // writable index + unsigned int & local_row (const unsigned int counter_index) + { return total_row_indices[counter_index].local_row; }; + + // returns the local row in the cell matrix + // associated with the counter_index-th entry + // in the list in the index_in_constraint-th + // position of constraints + unsigned int local_row (const unsigned int counter_index, const unsigned int index_in_constraint) const - { return (data_cache.get_entry(total_row_indices[loc_index].constraint_position) + { return (data_cache.get_entry(total_row_indices[counter_index].constraint_position) [index_in_constraint]).first; }; - double constraint_value (const unsigned int loc_index, + + // returns the value of the constraint in the + // counter_index-th entry in the list in the + // index_in_constraint-th position of + // constraints + double constraint_value (const unsigned int counter_index, const unsigned int index_in_constraint) const - { return (data_cache.get_entry(total_row_indices[loc_index].constraint_position) + { return (data_cache.get_entry(total_row_indices[counter_index].constraint_position) [index_in_constraint]).second; }; + + // returns whether there is one row with + // indirect contributions (i.e., there has + // been at least one constraint with + // non-trivial ConstraintLine) bool have_indirect_rows () const { return data_cache.element_size; }; + + // append an entry that is constrained. This + // means that there is one less row that data + // needs to be inserted into. void insert_constraint (const unsigned int constrained_local_dof) { --n_active_rows; total_row_indices[n_active_rows].local_row = constrained_local_dof; } + + // last row that was set to be constrained unsigned int last_constrained_local_row () { Assert (total_row_indices.back().global_row == numbers::invalid_unsigned_int, ExcInternalError()); return total_row_indices.back().local_row; }; + + // remove last entry in the list of global + // rows. Some elements at the end of the list + // total_row_indices are used to temporarily + // collect constrained dofs, but they should + // not have a global row index. Check that and + // eliminate the last such entry. void pop_back () { Assert (total_row_indices.back().global_row == numbers::invalid_unsigned_int, ExcInternalError()); total_row_indices.pop_back(); }; + + // store that the constraint_dof at the last + // position is inhomogeneous and needs to be + // considered further. void set_last_row_inhomogeneous () { std::swap (total_row_indices[n_active_rows+n_inhomogeneous_rows].local_row, total_row_indices.back().local_row); ++n_inhomogeneous_rows; }; + + // unsigned int inhomogeneity (unsigned int i) const { return total_row_indices[n_active_rows+i].local_row; }; + // a vector that contains all the global ids + // and the corresponding local ids as well as + // a pointer to that data where we store how + // to resolve constraints. std::vector total_row_indices; + + // holds the actual data from the constraints DataCache data_cache; + + // how many rows there are unsigned int n_active_rows; + + // the number of rows with inhomogeneous + // constraints unsigned int n_inhomogeneous_rows; }; // a function that appends an additional - // row to the list of values, or appends - // a value to an already existing + // row to the list of values, or appends a + // value to an already existing // row. Similar functionality as for - // std::map, - // but here done for a std::vector of - // data type Distributing, and much - // faster. + // std::map, but + // here done for a + // std::vector, much faster + // for short lists as we have them here inline void GlobalRowsFromLocal::insert_index (const unsigned int global_row, @@ -1200,21 +1301,20 @@ namespace internals data_cache.append_index (pos1->constraint_position, constraint); } + // this sort algorithm sorts + // std::vector, but does not + // take the constraints into account. this + // means that in case that constraints are + // already inserted, this function does not + // work as expected. Use shellsort, which + // is very fast in case the indices are + // already sorted (which is the usual case + // with DG elements), and not too slow in + // other cases inline void GlobalRowsFromLocal::sort () { - // this sort algorithm sorts a vector of - // Distributing elements, but does not - // take the constraints into - // account. this means that in case that - // constraints are already inserted, this - // function does not work as expected. - // we use shellsort, which is very fast in - // case the indices are already sorted - // (which is the usual case with DG - // elements), and not too slow in other - // cases unsigned int i, j, j2, temp, templ, istep; unsigned int step; @@ -1222,7 +1322,7 @@ namespace internals // constraints are really empty. const unsigned int length = size(); #ifdef DEBUG - for (unsigned int i=0; i inline void @@ -1288,8 +1394,8 @@ namespace internals block_indices = first_block; } - // transform row indices to local index - // space + // transform row indices to block-local + // index space for (unsigned int i=block_starts[1]; i instead of + // GlobalRowsFromLocal. Used in functions + // for sparsity patterns. template inline void @@ -1336,17 +1444,16 @@ namespace internals - // resolves constraints of one column - // at the innermost loop. goes through - // the origin of each global entry and - // finds out which data we need to collect + // resolves constraints of one column at + // the innermost loop. goes through the + // origin of each global entry and finds + // out which data we need to collect. inline double resolve_matrix_entry (const GlobalRowsFromLocal&global_rows, const unsigned int i, const unsigned int j, const unsigned int loc_row, - const FullMatrix &local_matrix, - const double* matrix_ptr) + const FullMatrix &local_matrix) { const unsigned int loc_col = global_rows.local_row(j); double col_val; @@ -1357,12 +1464,13 @@ namespace internals // set the value to zero. if (loc_row != numbers::invalid_unsigned_int) { - col_val = loc_col != numbers::invalid_unsigned_int ? matrix_ptr[loc_col] : 0; + col_val = ((loc_col != numbers::invalid_unsigned_int) ? + local_matrix(loc_row, loc_col) : 0); // account for indirect contributions by // constraints in column for (unsigned int p=0; p inline void @@ -1435,16 +1547,10 @@ namespace internals // to do some more checks. else { - const double * matrix_ptr = 0; - if (loc_row != numbers::invalid_unsigned_int) - { - Assert (loc_row < local_matrix.m(), ExcInternalError()); - matrix_ptr = &local_matrix(loc_row, 0); - } for (unsigned int j=column_start; j - inline - void add_value (const double value, - const unsigned int row, - const unsigned int column, - const unsigned int * col_ptr, - const bool are_on_diagonal, - unsigned int &counter, - number *val_ptr) + // specialized function that can write into + // the row of a SparseMatrix. + namespace dealiiSparseMatrix { - if (value != 0.) - { - Assert (col_ptr != 0, - typename SparseMatrix::ExcInvalidIndex (row, column)); - if (are_on_diagonal) - { - val_ptr[0] += value; - return; - } - while (col_ptr[counter] < column) - ++counter; - Assert (col_ptr[counter] == column, - typename SparseMatrix::ExcInvalidIndex(row, column)); - val_ptr[counter] += static_cast(value); - } + template + inline + void add_value (const double value, + const unsigned int row, + const unsigned int column, + const unsigned int * col_ptr, + const bool are_on_diagonal, + unsigned int &counter, + number *val_ptr) + { + if (value != 0.) + { + Assert (col_ptr != 0, + typename SparseMatrix::ExcInvalidIndex (row, column)); + if (are_on_diagonal) + { + val_ptr[0] += value; + return; + } + while (col_ptr[counter] < column) + ++counter; + Assert (col_ptr[counter] == column, + typename SparseMatrix::ExcInvalidIndex(row, column)); + val_ptr[counter] += static_cast(value); + } + } } - // similar as before, now with shortcut - // for deal.II sparse matrices. this lets - // use avoid using extra arrays, and does - // all the operations just in place + // similar as before, now with shortcut for + // deal.II sparse matrices. this lets use + // avoid using extra arrays, and does all the + // operations just in place, i.e., in the + // respective matrix row template inline void @@ -1528,12 +1638,13 @@ namespace internals const bool optimize_diagonal = sparsity.optimize_diagonal(); unsigned int counter = optimize_diagonal; - // distinguish three cases about what - // can happen (in order to avoid if() - // at the innermost loop position) - // for checking whether the diagonal is - // the first element of the row - if (!optimize_diagonal) // case 1: no diagonal optimization + // distinguish three cases about what can + // happen for checking whether the diagonal is + // the first element of the row. this avoids + // if statements at the innermost loop + // positions + + if (!optimize_diagonal) // case 1: no diagonal optimization in matrix { if (global_rows.have_indirect_rows() == false) { @@ -1545,20 +1656,21 @@ namespace internals { const unsigned int loc_col = global_rows.local_row(j); const double col_val = matrix_ptr[loc_col]; - add_value(col_val, row, global_rows.global_row(j), col_ptr, - false, counter, val_ptr); + dealiiSparseMatrix::add_value (col_val, row, + global_rows.global_row(j), + col_ptr, false, counter, + val_ptr); } } else { - const double * matrix_ptr = loc_row != numbers::invalid_unsigned_int ? - &local_matrix(loc_row, 0) : 0; for (unsigned int j=column_start; j &local_dof_indices, @@ -2057,6 +2176,8 @@ ConstraintMatrix:: + // Resolve the constraints from the vector and + // apply inhomogeneities. inline double ConstraintMatrix:: @@ -2066,11 +2187,13 @@ ConstraintMatrix:: const std::vector &local_dof_indices, const FullMatrix &local_matrix) const { - // Resolve the constraints from the vector and - // apply inhomogeneities. const unsigned int loc_row = global_rows.local_row(i); const unsigned int n_inhomogeneous_rows = global_rows.n_inhomogeneous_rows; double val = 0; + // has a direct contribution from some local + // entry. If we have inhomogeneous + // constraints, compute the contribution of + // the inhomogeneity in the current row. if (loc_row != numbers::invalid_unsigned_int) { val = local_vector(loc_row); @@ -2081,6 +2204,7 @@ ConstraintMatrix:: local_matrix(loc_row, global_rows.inhomogeneity(i))); } + // go through the indirect contributions for (unsigned int q=0; q void ConstraintMatrix:: @@ -2209,11 +2342,6 @@ distribute_local_to_global (const FullMatrix &local_matrix, VectorType &global_vector, internal::bool2type) const { - // similar function as above, but now - // specialized for block matrices. See - // the other function for additional - // comments. - const bool use_vectors = (local_vector.size() == 0 && global_vector.size() == 0) ? false : true; typedef typename MatrixType::value_type number; @@ -2264,10 +2392,11 @@ distribute_local_to_global (const FullMatrix &local_matrix, vals.resize (n_actual_dofs); } - // the basic difference to the - // non-block variant from now onwards - // is that we go through the blocks - // of the matrix separately. + // the basic difference to the non-block + // variant from now onwards is that we go + // through the blocks of the matrix + // separately, which allows us to set the + // block entries individually for (unsigned int block=0; block