/**
* Standard constructor. Sets
* @p{column} to
- * @p{invalid_entry}.
+ * @p{invalid}.
*/
Entry();
/**
* Comparison operator for finding.
*/
- bool operator==(const Entry&) const;
+// bool operator==(const Entry&) const;
/**
* Less than operator for sorting.
*/
- bool operator < (const Entry&) const;
+// bool operator < (const Entry&) const;
/**
* Non-existent column number.
*/
- static const unsigned int invalid_entry = static_cast<unsigned int>(-1);
+ static const unsigned int invalid = static_cast<unsigned int>(-1);
};
+ /**
+ * Structure for storing
+ * information on a matrix
+ * row. One object for each row
+ * will be stored in the matrix.
+ */
+ struct RowInfo
+ {
+ /**
+ * Constructor.
+ */
+ RowInfo (unsigned int start = Entry::invalid);
+
+ /**
+ * Index of first entry of
+ * the row in the data field.
+ */
+ unsigned int start;
+ /**
+ * Number of entries in this
+ * row.
+ */
+ unsigned short length;
+ /**
+ * Position of the diagonal
+ * element relative tor the
+ * start index.
+ */
+ unsigned short diagonal;
+ /**
+ * Value for non-existing diagonal.
+ */
+ static const unsigned short
+ invalid_diagonal = static_cast<unsigned short>(-1);
+ };
+
+
public:
/**
* Type of matrix entries. In analogy to
/**
* Constructor. Generates a
* matrix of the given size,
- * ready to be filled. The optional parameters
+ * ready to be filled. The
+ * optional parameters
+ * @p{default_row_length} and
+ * @p{default_increment} allow
+ * for preallocating
+ * memory. Providing these
+ * properly is essential for an
+ * efficient assembling of the
+ * matrix.
*/
explicit SparseMatrixEZ (unsigned int n_rows,
unsigned int n_columns = n_rows,
- unsigned int default_row_length = Entry::invalid_entry,
- unsigned int default_increment = Entry::invalid_entry);
+ unsigned int default_row_length = Entry::invalid,
+ unsigned int default_increment = Entry::invalid);
/**
* Destructor. Free all memory, but do not
* Reinitialize the sparse matrix
* to the dimensions provided.
* The matrix will have no
- * entries at this point.
+ * entries at this point. The
+ * optional parameters
+ * @p{default_row_length} and
+ * @p{default_increment} allow
+ * for preallocating
+ * memory. Providing these
+ * properly is essential for an
+ * efficient assembling of the
+ * matrix.
*/
virtual void reinit (unsigned int n_rows,
unsigned int n_columns = n_rows,
- unsigned int default_row_length = Entry::invalid_entry,
- unsigned int default_increment = Entry::invalid_entry);
+ unsigned int default_row_length = Entry::invalid,
+ unsigned int default_increment = Entry::invalid);
/**
* Release all memory and return
const Vector<somenumber> &src,
const number om = 1.) const;
- /**
- * Perform SSOR preconditioning
- * in-place. Apply the
- * preconditioner matrix without
- * copying to a second vector.
- * @p{omega} is the relaxation
- * parameter.
- */
- template <typename somenumber>
- void SSOR (Vector<somenumber> &v,
- const number omega = 1.) const;
-
- /**
- * Perform an SOR preconditioning in-place.
- * The result is $v = (\omega D - L)^{-1} v$.
- * @p{omega} is the damping parameter.
- */
- template <typename somenumber>
- void SOR (Vector<somenumber> &v,
- const number om = 1.) const;
-
- /**
- * Perform a transpose SOR preconditioning in-place.
- * The result is $v = (\omega D - L)^{-1} v$.
- * @p{omega} is the damping parameter.
- */
- template <typename somenumber>
- void TSOR (Vector<somenumber> &v,
- const number om = 1.) const;
-
- /**
- * Do one SOR step on @p{v}.
- * Performs a direct SOR step
- * with right hand side @p{b}.
- */
- template <typename somenumber>
- void SOR_step (Vector<somenumber> &v,
- const Vector<somenumber> &b,
- const number om = 1.) const;
-
- /**
- * Do one adjoint SOR step on
- * @p{v}. Performs a direct TSOR
- * step with right hand side @p{b}.
- */
- template <typename somenumber>
- void TSOR_step (Vector<somenumber> &v,
- const Vector<somenumber> &b,
- const number om = 1.) const;
-
- /**
- * Do one adjoint SSOR step on
- * @p{v}. Performs a direct SSOR
- * step with right hand side @p{b}
- * by performing TSOR after SOR.
- */
- template <typename somenumber>
- void SSOR_step (Vector<somenumber> &v,
- const Vector<somenumber> &b,
- const number om = 1.) const;
-
/**
* Print the matrix to the given
* stream, using the format
* of this object.
*/
unsigned int memory_consumption () const;
-
- /**
- * Exception
- */
- DeclException2 (ExcInvalidIndex,
- int, int,
- << "The entry with index <" << arg1 << ',' << arg2
- << "> does not exist.");
/**
- * Exception
+ * Exception for applying
+ * inverse-type operators to
+ * rectangular matrices.
*/
- DeclException0 (ExcMatrixNotSquare);
-
+ DeclException0(ExcNoSquareMatrix);
+
/**
- * Exception
+ * Exception for missing diagonal entry.
*/
- DeclException2 (ExcIteratorRange,
- int, int,
- << "The iterators denote a range of " << arg1
- << " elements, but the given number of rows was " << arg2);
+ DeclException0(ExcNoDiagonal);
-
private:
/**
*/
const Entry* locate (unsigned int row,
unsigned int col) const;
+
/**
* Find an entry or generate it.
*/
unsigned int n_columns;
/**
- * Start of indices rows. Points
- * into the data field.
+ * Info structure for each row.
*/
- std::vector<unsigned int> row_start;
+ std::vector<RowInfo> row_info;
/**
* Data storage.
inline
SparseMatrixEZ<number>::Entry::Entry()
:
- column(invalid_entry),
+ column(invalid),
value(0)
{}
-template <typename number>
-inline
-bool
-SparseMatrixEZ<number>::Entry::operator==(const Entry& e) const
-{
- return column == e.column;
-}
+// template <typename number>
+// inline
+// bool
+// SparseMatrixEZ<number>::Entry::operator==(const Entry& e) const
+// {
+// return column == e.column;
+// }
+// template <typename number>
+// inline
+// bool
+// SparseMatrixEZ<number>::Entry::operator<(const Entry& e) const
+// {
+// return column < e.column;
+// }
+
template <typename number>
inline
-bool
-SparseMatrixEZ<number>::Entry::operator<(const Entry& e) const
-{
- return column < e.column;
-}
-
+SparseMatrixEZ<number>::RowInfo::RowInfo(unsigned int start)
+ :
+ start(start), length(0), diagonal(invalid_diagonal)
+{}
//----------------------------------------------------------------------//
inline
unsigned int SparseMatrixEZ<number>::m () const
{
- return row_start.size() - 1;
+ return row_info.size();
};
Assert (row<m(), ExcIndexRange(row,0,m()));
Assert (col<n(), ExcIndexRange(col,0,n()));
- const unsigned int end = row_start[row+1];
- for (unsigned int i=row_start[row];i<end;++i)
+ const RowInfo& r = row_info[row];
+ const unsigned int end = r.start + r.length;
+ for (unsigned int i=r.start;i<end;++i)
{
const Entry * const entry = &data[i];
if (entry->column == col)
return entry;
- if (entry->column == Entry::invalid_entry)
+ if (entry->column == Entry::invalid)
return 0;
}
return 0;
Assert (row<m(), ExcIndexRange(row,0,m()));
Assert (col<n(), ExcIndexRange(col,0,n()));
- const unsigned int end = row_start[row+1];
- for (unsigned int i=row_start[row];i<end;++i)
+ RowInfo& r = row_info[row];
+ const unsigned int end = r.start + r.length;
+
+ unsigned int i = r.start;
+ if (r.diagonal != RowInfo::invalid_diagonal && col >= row)
+ i += r.diagonal;
+ // Find position of entry
+ while (i<end && data[i].column < col) ++i;
+
+ Entry* entry = &data[i];
+ // entry found
+ if (entry->column == col)
+ return entry;
+
+ // Now, we must insert the new
+ // entry and move all successive
+ // entries back.
+
+ // If no more space is available
+ // for this row, insert new
+ // elements into the vector.
+ if (row != row_info.size()-1)
{
- Entry* entry = &data[i];
- // entry found
- if (entry->column == col)
- return entry;
- // entry does not exist,
- // create it
- if (entry->column > col)
+ if (end >= row_info[row+1].start)
+ {
+ // Insert new entries
+ data.insert(data.begin()+end, increment, Entry());
+ entry = &data[i];
+ // Update starts of
+ // following rows
+ for (unsigned int rn=row+1;rn<row_info.size();++rn)
+ row_info[rn].start += increment;
+ }
+ } else {
+ if (end >= data.size())
{
- // Save original entry
- Entry temp = data[i];
- // Insert new entry here to
- // make sure all entries
- // are ordered by column
- // index
- entry->column = col;
- entry->value = 0;
-
- // Move all entries in this
- // row up by one
- for (unsigned int j = i+1;j<end;++j)
- {
- Entry temp2 = data[j];
- data[j] = temp;
- temp = temp2;
- if (temp.column == Entry::invalid_entry)
- break;
- }
- // Extend row if there is
- // still a valid entry left
- if (temp.column != Entry::invalid_entry)
- {
- // Insert new entries
- data.insert(data.begin()+end, increment, Entry());
- // Update starts of
- // following rows
- for (unsigned int r=row+1;r<row_start.size();++r)
- row_start[r] += increment;
- // Insert missing entry
- data[end] = temp;
- }
- return entry;
+ // Here, appending a block
+ // does not increase
+ // performance.
+ data.push_back(Entry());
+ entry = &data[i];
}
}
- // Row is full, but entry not
- // found. Therefore, we must
- // extend again
- data.insert(data.begin()+end, increment, Entry());
- // Update starts of
- // following rows
- for (unsigned int r=row+1;r<row_start.size();++r)
- row_start[r] += increment;
- // Insert missing entry
- data[end].column = col;
- return &data[end];
+ // Save original entry
+ Entry temp = *entry;
+ // Insert new entry here to
+ // make sure all entries
+ // are ordered by column
+ // index
+ entry->column = col;
+ entry->value = 0;
+ // Update row_info
+ ++r.length;
+ if (col == row)
+ r.diagonal = i - r.start;
+ else if (col<row && r.diagonal!= RowInfo::invalid_diagonal)
+ ++r.diagonal;
+
+ if (i == end)
+ return entry;
+
+ // Move all entries in this
+ // row up by one
+ for (unsigned int j = i+1; j < end; ++j)
+ {
+ // There should be no invalid
+ // entry below end
+ Assert (data[j].column != Entry::invalid, ExcInternalError());
+ Entry temp2 = data[j];
+ data[j] = temp;
+ temp = temp2;
+ }
+ Assert (data[end].column == Entry::invalid, ExcInternalError());
+ data[end] = temp;
+
+ return entry;
}
unsigned int default_row_length,
unsigned int default_increment)
{
- if (default_row_length == Entry::invalid_entry)
+ clear();
+ if (default_row_length == Entry::invalid)
default_row_length = 5;
- if (default_increment == Entry::invalid_entry)
+ if (default_increment == Entry::invalid)
default_increment = 4;
if (default_increment == 0)
default_increment = 4;
increment = default_increment;
n_columns = n_cols;
- row_start.resize(n_rows+1);
+ row_info.resize(n_rows);
data.reserve(default_row_length * n_rows + n_rows * increment);
data.resize(default_row_length * n_rows);
- for (unsigned int i=0;i<=n_rows;++i)
- row_start[i] = i * default_row_length;
+ for (unsigned int i=0;i<n_rows;++i)
+ row_info[i].start = i * default_row_length;
}
SparseMatrixEZ<number>::clear()
{
n_columns = 0;
- row_start.resize(0);
+ row_info.resize(0);
data.resize(0);
}
Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
- const unsigned int end_row = row_start.size() - 1;
- for (unsigned int i = 0; i < end_row;++i)
+ const unsigned int end_row = row_info.size();
+ for (unsigned int row = 0; row < end_row; ++row)
{
- unsigned int index = row_start[i];
- unsigned int end = row_start[i+1];
+ const RowInfo& ri = row_info[row];
+ typename std::vector<Entry>::const_iterator
+ entry = data.begin() + ri.start;
double s = 0.;
- for (;index != end && data[index].column != Entry::invalid_entry;++index)
+ for (unsigned short i=0;i<ri.length;++i,++entry)
{
- const Entry& entry = data[index];
- s += entry.value * src(entry.column);
+ Assert (entry->column != Entry::invalid,
+ ExcInternalError());
+ s += entry->value * src(entry->column);
}
- dst(i) = s;
+ dst(row) = s;
}
}
Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
- const unsigned int end_row = row_start.size() - 1;
- for (unsigned int i = 0; i < end_row;++i)
+ const unsigned int end_row = row_info.size();
+ for (unsigned int row = 0; row < end_row; ++row)
{
- unsigned int index = row_start[i];
- unsigned int end = row_start[i+1];
+ const RowInfo& ri = row_info[i];
+ typename std::vector<Entry>::const_iterator
+ entry = data.begin() + ri.start;
double s = 0.;
- for (;index != end && data[index].column != Entry::invalid_entry;++index)
+ for (unsigned short i=0;i<ri.length;++i,++entry)
{
- const Entry& entry = data[index];
- s += entry.value * src(entry.column);
+ Assert (entry->column != Entry::invalid,
+ ExcInternalError());
+ s += entry->value * src(entry->column);
}
- dst(i) += s;
+ dst(row) += s;
}
}
Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
- const unsigned int end_row = row_start.size() - 1;
- for (unsigned int i = 0; i < end_row;++i)
+ const unsigned int end_row = row_info.size();
+ for (unsigned int row = 0; row < end_row; ++row)
{
- unsigned int index = row_start[i];
- unsigned int end = row_start[i+1];
- for (;index != end && data[index].column != Entry::invalid_entry;++index)
+ const RowInfo& ri = row_info[row];
+ typename std::vector<Entry>::const_iterator
+ entry = data.begin() + ri.start;
+ for (unsigned short i=0;i<ri.length;++i,++entry)
{
- const Entry& entry = data[index];
- dst(entry.column) += entry.value * src(i);
+ Assert (entry->column != Entry::invalid,
+ ExcInternalError());
+ dst(entry->column) += entry->value * src(row);
}
}
}
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_Jacobi (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om) const
+{
+ Assert (m() == n(), ExcNoSquareMatrix());
+ Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+ Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+ somenumber *dst_ptr = dst.begin();
+ const somenumber *src_ptr = src.begin();
+ typename std::vector<RowInfo>::const_iterator ri = row_info.begin();
+ const typename std::vector<RowInfo>::const_iterator end = row_info.end();
+
+ for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
+ {
+ Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+ *dst_ptr = om * *src_ptr / data[ri->start + ri->diagonal].value;
+ }
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_SOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om) const
+{
+ Assert (m() == n(), ExcNoSquareMatrix());
+ Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+ Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+ somenumber *dst_ptr = dst.begin();
+ const somenumber *src_ptr = src.begin();
+ typename std::vector<RowInfo>::const_iterator ri = row_info.begin();
+ const typename std::vector<RowInfo>::const_iterator end = row_info.end();
+
+ for (; ri != end; ++dst_ptr, ++src_ptr, ++ri)
+ {
+ Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+ number s = *src_ptr;
+ const unsigned int end_row = ri->start + ri->diagonal;
+ for (unsigned int i=ri->start;i<end_row;++i)
+ s -= data[i].value * dst(data[i].column);
+
+ *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
+ }
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_TSOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om) const
+{
+ Assert (m() == n(), ExcNoSquareMatrix());
+ Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+ Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+ somenumber *dst_ptr = dst.begin()+dst.size()-1;
+ const somenumber *src_ptr = src.begin()+src.size()-1;
+ typename std::vector<RowInfo>::const_reverse_iterator
+ ri = row_info.rbegin();
+ const typename std::vector<RowInfo>::const_reverse_iterator
+ end = row_info.rend();
+
+ for (; ri != end; --dst_ptr, --src_ptr, ++ri)
+ {
+ Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+ number s = *src_ptr;
+ const unsigned int end_row = ri->start + ri->length;
+ for (unsigned int i=ri->start+ri->diagonal+1;i<end_row;++i)
+ s -= data[i].value * dst(data[i].column);
+
+ *dst_ptr = om * s / data[ri->start + ri->diagonal].value;
+ }
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrixEZ<number>::precondition_SSOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om) const
+{
+ Assert (m() == n(), ExcNoSquareMatrix());
+ Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+ Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+ somenumber *dst_ptr = dst.begin();
+ const somenumber *src_ptr = src.begin();
+ typename std::vector<RowInfo>::const_iterator ri;
+ const typename std::vector<RowInfo>::const_iterator end = row_info.end();
+
+ // Forward
+ for (ri = row_info.begin(); ri != end; ++dst_ptr, ++src_ptr, ++ri)
+ {
+ Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal());
+ number s = *src_ptr;
+ const unsigned int end_row = ri->start + ri->diagonal;
+ for (unsigned int i=ri->start;i<end_row;++i)
+ s -= om * data[i].value * dst(data[i].column);
+
+ *dst_ptr = s / data[ri->start + ri->diagonal].value;
+ }
+ // Diagonal
+ dst_ptr = dst.begin();
+ for (ri = row_info.begin(); ri != end; ++dst_ptr, ++ri)
+ *dst_ptr *= (2.-om) * data[ri->start + ri->diagonal].value;
+
+ // Backward
+ typename std::vector<RowInfo>::const_reverse_iterator rri;
+ const typename std::vector<RowInfo>::const_reverse_iterator
+ rend = row_info.rend();
+ dst_ptr = dst.begin()+dst.size()-1;
+ for (rri = row_info.rbegin(); rri != rend; --dst_ptr, ++rri)
+ {
+ const unsigned int end_row = rri->start + rri->length;
+ for (unsigned int i=rri->start+rri->diagonal+1;i<end_row;++i)
+ *dst_ptr -= om * data[i].value * dst(data[i].column);
+ *dst_ptr /= data[rri->start + rri->diagonal].value;
+ }
+};
+
+
+
template <typename number>
unsigned int
SparseMatrixEZ<number>::memory_consumption() const
Vector<double> f(A.m());
GrowingVectorMemory<> mem;
- SolverControl control(20, 1.e-3, false);
+ SolverControl control(10, 1.e-13, false);
SolverRichardson<> rich(control, mem, .01);
+ SolverRichardson<> prich(control, mem, 1.);
PreconditionIdentity prec;
-
+ PreconditionJacobi<MATRIX> jacobi;
+ jacobi.initialize(A, .5);
+ PreconditionSOR<MATRIX> sor;
+ sor.initialize(A, 1.2);
+ PreconditionSSOR<MATRIX> ssor;
+ ssor.initialize(A, 1.2);
+
u = 0.;
f = 1.;
- try
- {
- rich.solve(A, u, f, prec);
- }
- catch (...)
- {
- }
+ try { rich.solve(A, u, f, prec); } catch (...) {}
+ try { prich.solve(A, u, f, jacobi); } catch (...) {}
+ try { prich.solve(A, u, f, ssor); } catch (...) {}
+ try { prich.solve(A, u, f, sor); } catch (...) {}
+
+ u = 0.;
deallog << "Transpose" << std::endl;
- try
- {
- rich.Tsolve(A, u, f, prec);
- }
- catch (...)
- {
- }
+ try { rich.Tsolve(A, u, f, prec); } catch (...) {}
+ try { prich.Tsolve(A, u, f, jacobi); } catch (...) {}
+ try { prich.Tsolve(A, u, f, ssor); } catch (...) {}
+ try { prich.Tsolve(A, u, f, sor); } catch (...) {}
deallog.pop();
}
// Switch between regression test
// and benchmark
#ifdef DEBUG
- deallog.depth_console(0);
- const unsigned int size = 10;
+ deallog.depth_console(3);
+ const unsigned int size = 5;
+ const unsigned int row_length = 3;
#else
deallog.depth_console(1000);
deallog.log_execution_time(true);
deallog.log_time_differences(true);
const unsigned int size = 500;
+ const unsigned int row_length = 9;
#endif
FDMatrix testproblem (size, size);
structure.compress();
SparseMatrix<double> A(structure);
deallog << "Assemble" << std::endl;
- testproblem.five_point(A);
+ testproblem.five_point(A, true);
check_vmult_quadratic(A, "5-SparseMatrix<double>");
- SparseMatrixEZ<double> E(dim,dim,5,1);
+ SparseMatrixEZ<double> E(dim,dim,row_length,2);
deallog << "Assemble" << std::endl;
- testproblem.five_point(E);
+ testproblem.five_point(E, true);
check_vmult_quadratic(E, "5-SparseMatrixEZ<double>");
A.clear();
check_vmult_quadratic(A, "9-SparseMatrix<double>");
E.clear();
- E.reinit(dim,dim,9,2);
+ E.reinit(dim,dim,row_length,2);
deallog << "Assemble" << std::endl;
testproblem.nine_point(E);
check_vmult_quadratic(E, "9-SparseMatrixEZ<double>");
template<typename MATRIX>
void
-FDMatrix::nine_point(MATRIX& A) const
+FDMatrix::nine_point(MATRIX& A, bool nonsymmetric) const
{
for(unsigned int i=0;i<=ny-2;i++)
{
template<typename MATRIX>
void
-FDMatrix::five_point(MATRIX& A) const
+FDMatrix::five_point(MATRIX& A, bool nonsymmetric) const
{
for(unsigned int i=0;i<=ny-2;i++)
{
{
// Number of the row to be entered
unsigned int row = j+(nx-1)*i;
- A.set(row, row, 4.);
+ if (nonsymmetric)
+ A.set(row, row, 5.);
+ else
+ A.set(row, row, 4.);
if (j>0)
{
- A.set(row-1, row, -1.);
+ if (nonsymmetric)
+ A.set(row-1, row, -2.);
+ else
+ A.set(row-1, row, -1.);
A.set(row, row-1, -1.);
}
if (j<nx-2)
* Fill the matrix with values.
*/
template <typename MATRIX>
- void five_point(MATRIX&) const;
+ void five_point(MATRIX&, bool nonsymmetric = false) const;
/**
* Fill the matrix with values.
*/
template <typename MATRIX>
- void nine_point(MATRIX&) const;
+ void nine_point(MATRIX&, bool nonsymmetric = false) const;
template <typename number>
void gnuplot_print(std::ostream&, const Vector<number>&) const;
*/
unsigned int ny;
};
-