// $Id$
// Version: $Name$
//
-// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Assert (s.colnums == 0, ExcInvalidConstructorCall());
Assert (s.rows == 0, ExcInvalidConstructorCall());
Assert (s.cols == 0, ExcInvalidConstructorCall());
-
+
reinit (0,0,0, false);
}
SparsityPattern::SparsityPattern (const unsigned int m,
const unsigned int n,
const unsigned int max_per_row,
- const bool optimize_diag)
+ const bool optimize_diag)
:
max_dim(0),
max_vec_len(0),
const unsigned int m,
const unsigned int n,
const std::vector<unsigned int>& row_lengths,
- const bool optimize_diag)
+ const bool optimize_diag)
:
max_dim(0),
max_vec_len(0),
SparsityPattern::SparsityPattern (
const unsigned int m,
const std::vector<unsigned int>& row_lengths,
- const bool optimize_diag)
+ const bool optimize_diag)
:
max_dim(0),
max_vec_len(0),
Assert (original.rows==original.cols, ExcNotQuadratic());
Assert (original.optimize_diagonal(), ExcNotQuadratic());
Assert (original.is_compressed(), ExcNotCompressed());
-
+
reinit (original.rows, original.cols, max_per_row,
original.optimize_diagonal());
for (unsigned int row=0; row<original.rows; ++row)
{
// copy the elements of this row
- // of the other object
- //
+ // of the other object
+ //
// note that the first object actually
// is the main-diagonal element,
// which we need not copy
original_row_end,
row-extra_off_diagonals) :
original_row_start);
-
+
const unsigned int * const
original_first_after_side_diagonals
= (row < rows-extra_off_diagonals-1 ?
const bool optimize_diag)
{
AssertDimension (row_lengths.size(), m);
-
+
rows = m;
cols = n;
delete[] colnums;
colnums = 0;
}
-
+
max_vec_len = vec_len;
colnums = new unsigned int[max_vec_len];
}
-
+
max_row_length = (row_lengths.size() == 0 ?
0 :
std::min (*std::max_element(row_lengths.begin(),
if (diagonal_optimized && (max_row_length==0) && (m!=0))
max_row_length = 1;
-
+
// allocate memory for the rowstart
// values, if necessary. even
// though we re-set the pointers
delete[] colnums;
colnums = 0;
}
-
+
max_vec_len = vec_len;
colnums = new unsigned int[max_vec_len];
}
- // set the rowstart array
+ // set the rowstart array
rowstart[0] = 0;
for (unsigned int i=1; i<=rows; ++i)
rowstart[i] = rowstart[i-1] +
// do nothing if already compressed
if (compressed)
return;
-
+
unsigned int next_free_entry = 0,
next_row_start = 0,
row_length = 0;
// reserve temporary storage to
// store the entries of one row
std::vector<unsigned int> tmp_entries (max_row_length);
-
+
// Traverse all rows
for (unsigned int line=0; line<rows; ++line)
{
// now @p{rowstart} is
// the number of entries in
// this line
-
+
// Sort only beginning at the
// second entry, if optimized
// storage of diagonal entries
// is on.
-
+
// if this line is empty or has
// only one entry, don't sort
if (row_length > 1)
? tmp_entries.begin()+1
: tmp_entries.begin(),
tmp_entries.begin()+row_length);
-
+
// insert column numbers
// into the new field
for (unsigned int j=0; j<row_length; ++j)
// store the size
max_vec_len = nonzero_elements;
-
+
compressed = true;
}
template <typename CSP>
void
SparsityPattern::copy_from (const CSP &csp,
- const bool optimize_diag)
+ const bool optimize_diag)
{
// first determine row lengths for
// each row. if the matrix is
reinit (csp.n_rows(), csp.n_cols(), row_lengths, do_diag_optimize);
// now enter all the elements into
- // the matrix. note that if the
- // matrix is quadratic, then we
- // already have the diagonal
- // element preallocated
- for (unsigned int row = 0; row<csp.n_rows(); ++row)
- {
- unsigned int *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
- typename CSP::row_iterator col_num = csp.row_begin (row),
- end_row = csp.row_end (row);
-
- for (; col_num != end_row; ++col_num)
- {
- const unsigned int col = *col_num;
- if ((col!=row) || !do_diag_optimize)
- *cols++ = col;
- }
- }
+ // the matrix, if there are
+ // any. note that if the matrix is
+ // quadratic, then we already have
+ // the diagonal element
+ // preallocated
+ if (n_rows() * n_cols() > 0)
+ for (unsigned int row = 0; row<csp.n_rows(); ++row)
+ {
+ unsigned int *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
+ typename CSP::row_iterator col_num = csp.row_begin (row),
+ end_row = csp.row_end (row);
+
+ for (; col_num != end_row; ++col_num)
+ {
+ const unsigned int col = *col_num;
+ if ((col!=row) || !do_diag_optimize)
+ *cols++ = col;
+ }
+ }
// do not need to compress the sparsity
// pattern since we already have
(matrix.m() == matrix.n()))
++entries_per_row[row];
}
-
+
reinit (matrix.m(), matrix.n(), entries_per_row, optimize_diag);
// now set entries
unsigned int
-SparsityPattern::max_entries_per_row () const
+SparsityPattern::max_entries_per_row () const
{
// if compress() has not yet been
// called, we can get the maximum
SparsityPattern::operator () (const unsigned int i,
const unsigned int j) const
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (i<rows, ExcIndexRange(i,0,rows));
Assert (j<cols, ExcIndexRange(j,0,cols));
Assert (compressed, ExcNotCompressed());
SparsityPattern::add (const unsigned int i,
const unsigned int j)
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (i<rows, ExcIndexRange(i,0,rows));
Assert (j<cols, ExcIndexRange(j,0,cols));
Assert (compressed==false, ExcMatrixIsCompressed());
SparsityPattern::exists (const unsigned int i,
const unsigned int j) const
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (i<rows, ExcIndexRange(i,0,rows));
Assert (j<cols, ExcIndexRange(j,0,cols));
void
-SparsityPattern::symmetrize ()
+SparsityPattern::symmetrize ()
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
Assert (compressed==false, ExcMatrixIsCompressed());
// Note that we only require a
// quadratic matrix here, no
void
SparsityPattern::print (std::ostream &out) const
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
AssertThrow (out, ExcIO());
out << ',' << colnums[j];
out << ']' << std::endl;
}
-
+
AssertThrow (out, ExcIO());
}
void
SparsityPattern::print_gnuplot (std::ostream &out) const
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
AssertThrow (out, ExcIO());
-
+
for (unsigned int i=0; i<rows; ++i)
for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
if (colnums[j] != invalid_entry)
unsigned int
SparsityPattern::bandwidth () const
{
- Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+ Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
unsigned int b=0;
for (unsigned int i=0; i<rows; ++i)
for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
void
-SparsityPattern::block_write (std::ostream &out) const
+SparsityPattern::block_write (std::ostream &out) const
{
AssertThrow (out, ExcIO());
reinterpret_cast<const char*>(&colnums[max_vec_len])
- reinterpret_cast<const char*>(&colnums[0]));
out << ']';
-
+
AssertThrow (out, ExcIO());
}
rowstart = new std::size_t[max_dim+1];
colnums = new unsigned int[max_vec_len];
-
+
// then read data
in.read (reinterpret_cast<char*>(&rowstart[0]),
reinterpret_cast<char*>(&rowstart[max_dim+1])
const unsigned int*,
const bool);
#ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
-template void SparsityPattern::add_entries<std::vector<unsigned int>::const_iterator>
+template void SparsityPattern::add_entries<std::vector<unsigned int>::const_iterator>
(const unsigned int,
std::vector<unsigned int>::const_iterator,
std::vector<unsigned int>::const_iterator,
const bool);
#endif
-template void SparsityPattern::add_entries<std::vector<unsigned int>::iterator>
+template void SparsityPattern::add_entries<std::vector<unsigned int>::iterator>
(const unsigned int,
std::vector<unsigned int>::iterator,
std::vector<unsigned int>::iterator,
--- /dev/null
+//----------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2006, 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+
+
+// Test reinit(BlockIndices...)
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <lac/block_sparsity_pattern.h>
+
+#include <iomanip>
+#include <fstream>
+
+int main()
+{
+ std::ofstream logfile("block_sparsity_pattern_01/output");
+ logfile.setf(std::ios::fixed);
+ deallog << std::setprecision(2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ BlockSparsityPattern sparsity;
+ std::vector<unsigned int> row_blocks(4);
+ row_blocks[0] = 4;
+ row_blocks[1] = 5;
+ row_blocks[2] = 1;
+ row_blocks[3] = 4;
+ std::vector<unsigned int> col_blocks(3);
+ col_blocks[0] = 2;
+ col_blocks[1] = 3;
+ col_blocks[2] = 2;
+ BlockIndices rows(row_blocks);
+ BlockIndices cols(col_blocks);
+
+ std::vector<std::vector<unsigned int> >
+ row_length(cols.size(),
+ std::vector<unsigned int>(rows.total_size()));
+ for (unsigned int jb=0;jb<row_length.size();++jb)
+ for (unsigned int i=0;i<row_length[jb].size();++i)
+ {
+ const unsigned int d = col_blocks[jb]-1;
+ row_length[jb][i] = (i+1) % d +1;
+ }
+
+ for (unsigned int j=0;j<row_length.size();++j)
+ {
+ for (unsigned int i=0;i<row_length[j].size();++i)
+ deallog << ' ' << row_length[j][i];
+ deallog << std::endl;
+ }
+
+ sparsity.reinit(rows, cols, row_length);
+
+ for (unsigned int ib = 0;ib<rows.size();++ib)
+ for (unsigned int i=0;i<rows.block_size(ib);++i)
+ {
+ const unsigned int ii = rows.local_to_global(ib,i);
+ for (unsigned int jb = 0;jb<cols.size();++jb)
+ for (unsigned int j=0;j<row_length[jb][ii];++j)
+ sparsity.add(ii, cols.local_to_global(jb,j));
+ }
+
+ sparsity.print(logfile);
+
+ return 0;
+}