From: kronbichler Date: Tue, 21 Jan 2014 08:30:25 +0000 (+0000) Subject: Test extract_row_copy, fix two bugs X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9bfb7f475b066b102b6eba3a9703376dfc4b38d5;p=dealii-svn.git Test extract_row_copy, fix two bugs git-svn-id: https://svn.dealii.org/trunk@32276 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h index b493ed3f26..188a020e40 100644 --- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h +++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h @@ -605,24 +605,46 @@ ChunkSparseMatrix::extract_row_copy (const size_type row, size_type *column_indices, number *values) const { - row_length = cols->row_length(row); + AssertIndexRange(cols->row_length(row), array_length+1); + AssertIndexRange(row, m()); const unsigned int chunk_size = cols->get_chunk_size(); const size_type reduced_row = row/chunk_size; - AssertIndexRange(cols->sparsity_pattern.row_length(reduced_row)*chunk_size, - array_length+1); SparsityPattern::iterator it = cols->sparsity_pattern.begin(reduced_row), itend = cols->sparsity_pattern.end(reduced_row); const number *val_ptr = &val[(it-cols->sparsity_pattern.begin(0))*chunk_size*chunk_size +(row%chunk_size)*chunk_size]; - for ( ; it < itend; ++it) + + // find out if we did padding and if this row is affected by it + if (cols->n_cols() % chunk_size == 0) + { + for ( ; it < itend; ++it) + { + for (unsigned int c=0; ccolumn()*chunk_size+c; + } + val_ptr += chunk_size*chunk_size; + } + row_length = chunk_size * (cols->sparsity_pattern.row_length(reduced_row)); + } + else { - for (unsigned int c=0; cn_cols() % chunk_size; + row_length = 0; + for ( ; it < itend; ++it) { - *values++ = val_ptr[c]; - *column_indices++ = it->column()*chunk_size+c; + const unsigned int next_chunk_size = + (it->column()==cols->sparsity_pattern.n_cols()-1) ? + last_chunk_size : chunk_size; + for (unsigned int c=0; ccolumn()*chunk_size+c; + } + val_ptr += chunk_size*chunk_size; } - val_ptr += chunk_size*chunk_size; } } diff --git a/deal.II/source/lac/sparsity_pattern.cc b/deal.II/source/lac/sparsity_pattern.cc index 617b820c16..63841efa32 100644 --- a/deal.II/source/lac/sparsity_pattern.cc +++ b/deal.II/source/lac/sparsity_pattern.cc @@ -190,42 +190,29 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original, reinit (original.rows, original.cols, max_per_row); - // now copy the entries from - // the other object + // now copy the entries from the other object for (size_type row=0; row extra_off_diagonals ? @@ -242,9 +229,8 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original, row+extra_off_diagonals) : original_row_end); - // find first free slot. the - // first slot in each row is - // the diagonal element + // find first free slot. the first slot in each row is the diagonal + // element size_type *next_free_slot = &colnums[rowstart[row]] + 1; // copy elements before side-diagonals @@ -265,12 +251,9 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original, original_row_end, next_free_slot); - // this error may happen if the - // sum of previous elements per row - // and those of the new diagonals - // exceeds the maximum number of - // elements per row given to this - // constructor + // this error may happen if the sum of previous elements per row and + // those of the new diagonals exceeds the maximum number of elements per + // row given to this constructor Assert (next_free_slot <= &colnums[rowstart[row+1]], ExcNotEnoughSpace (0,rowstart[row+1]-rowstart[row])); }; @@ -310,8 +293,7 @@ SparsityPattern::reinit (const size_type m, const unsigned int max_per_row, const bool) { - // simply map this function to the - // other @p{reinit} function + // simply map this function to the other @p{reinit} function const std::vector row_lengths (m, max_per_row); reinit (m, n, row_lengths); } @@ -323,8 +305,7 @@ SparsityPattern::reinit (const size_type m, const size_type n, const unsigned int max_per_row) { - // simply map this function to the - // other @p{reinit} function + // simply map this function to the other @p{reinit} function const std::vector row_lengths (m, max_per_row); reinit (m, n, row_lengths); } @@ -360,31 +341,23 @@ SparsityPattern::reinit (const size_type m, rowstart = 0; colnums = 0; max_vec_len = max_dim = rows = cols = 0; - // if dimension is zero: ignore - // max_per_row + // if dimension is zero: ignore max_per_row max_row_length = 0; compressed = false; return; } - // first, if the matrix is - // quadratic, we - // will have to make sure that each - // row has at least one entry for - // the diagonal element. make this - // more obvious by having a - // variable which we can query + // first, if the matrix is quadratic, we will have to make sure that each + // row has at least one entry for the diagonal element. make this more + // obvious by having a variable which we can query store_diagonal_first_in_row = (m == n); - // find out how many entries we - // need in the @p{colnums} array. if - // this number is larger than - // @p{max_vec_len}, then we will need - // to reallocate memory + // find out how many entries we need in the @p{colnums} array. if this + // number is larger than @p{max_vec_len}, then we will need to reallocate + // memory // - // note that the number of elements - // per row is bounded by the number - // of columns + // note that the number of elements per row is bounded by the number of + // columns // std::size_t vec_len = 0; for (size_type i=0; i max_dim) { if (rowstart) @@ -448,8 +409,7 @@ SparsityPattern::reinit (const size_type m, rowstart = new std::size_t[max_dim+1]; } - // allocate memory for the column - // numbers if necessary + // allocate memory for the column numbers if necessary if (vec_len > max_vec_len) { if (colnums) @@ -475,14 +435,11 @@ SparsityPattern::reinit (const size_type m, ((vec_len == 1) && (rowstart[rows] == 0)), ExcInternalError()); - // preset the column numbers by a - // value indicating it is not in - // use + // preset the column numbers by a value indicating it is not in use std::fill_n (&colnums[0], vec_len, invalid_entry); - // if diagonal elements are - // special: let the first entry in - // each row be the diagonal value + // if diagonal elements are special: let the first entry in each row be the + // diagonal value if (store_diagonal_first_in_row) for (size_type i=0; i tmp_entries (max_row_length); // Traverse all rows for (size_type line=0; line 1) std::sort ((store_diagonal_first_in_row) ? tmp_entries.begin()+1 : tmp_entries.begin(), tmp_entries.begin()+row_length); - // insert column numbers - // into the new field + // insert column numbers into the new field for (size_type j=0; j void SparsityPattern::copy_from (const CSP &csp) { - // first determine row lengths for - // each row. if the matrix is - // quadratic, then we might have to - // add an additional entry for the - // diagonal, if that is not yet - // present. as we have to call - // compress anyway later on, don't - // bother to check whether that - // diagonal entry is in a certain - // row or not + // first determine row lengths for each row. if the matrix is quadratic, + // then we might have to add an additional entry for the diagonal, if that + // is not yet present. as we have to call compress anyway later on, don't + // bother to check whether that diagonal entry is in a certain row or not const bool do_diag_optimize = (csp.n_rows() == csp.n_cols()); std::vector row_lengths (csp.n_rows()); for (size_type i=0; i &matrix, template void SparsityPattern::copy_from (const FullMatrix &matrix) { - // first init with the number of entries - // per row. if this matrix is square then we - // also have to allocate memory for the - // diagonal entry, unless we have already - // counted it + // first init with the number of entries per row. if this matrix is square + // then we also have to allocate memory for the diagonal entry, unless we + // have already counted it std::vector entries_per_row (matrix.m(), 0); for (size_type row=0; row(rowstart[i]-rowstart[i-1])); return m; @@ -806,30 +718,21 @@ SparsityPattern::operator () (const size_type i, Assert (j(i) << std::endl; AssertThrow (out, ExcIO()); @@ -1084,8 +968,7 @@ SparsityPattern::bandwidth () const b = std::abs(static_cast(i-colnums[j])); } else - // leave if at the end of - // the entries of this line + // leave if at the end of the entries of this line break; return b; } @@ -1096,8 +979,7 @@ SparsityPattern::block_write (std::ostream &out) const { AssertThrow (out, ExcIO()); - // first the simple objects, - // bracketed in [...] + // first the simple objects, bracketed in [...] out << '[' << max_dim << ' ' << rows << ' ' << cols << ' ' diff --git a/tests/bits/chunk_sparse_matrix_16.cc b/tests/bits/chunk_sparse_matrix_16.cc new file mode 100644 index 0000000000..f9a99c45f7 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_16.cc @@ -0,0 +1,115 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2004 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// compare matrix-vector product in chunk sparse matrix with usual vmult and +// extract_row_copy. + +#include "../tests.h" +#include +#include +#include +#include + + +void test (const unsigned int chunk_size) +{ + deallog << "Chunk size = " << chunk_size << std::endl; + + for (unsigned int n_cols = 4; n_cols<7; ++n_cols) + { + deallog << "n_cols = " << n_cols << std::endl; + ChunkSparsityPattern sp (5,n_cols,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j m(sp); + + // first set a few entries + for (unsigned int i=0; i src(m.n()), dst(m.m()), dst_ref(m.m()); + for (unsigned int i=0; i indices(sp.max_entries_per_row()); + std::vector values(sp.max_entries_per_row()); + for (unsigned int i=0; i