const number& value)
{
// Store end of vector
- const std::vector<Entry>::const_iterator end_col = end();
+ const typename std::vector<Entry>::iterator end_col = end();
// Create Entry for inserting
const Entry e(column, value);
// Find position for inserting
// should return first Entry with
// higher column number.
- std::vector<Entry>::iterator col = lower_bound(begin(), end, e);
+ typename std::vector<Entry>::iterator col = lower_bound(begin(), end_col, e);
// Insert Entry if column did not exist.
// Edit existing entry otherwise.
values.push_back(e);
else
if (col->column == column)
- col->values = value;
+ col->value = value;
else
values.insert(col, e);
}
SparseMatrixEZ<number>::Row::add(unsigned int column,
const number& value)
{
- const std::vector<Entry>::const_iterator end_col = end();
+ const typename std::vector<Entry>::const_iterator end_col = end();
// Create Entry for inserting
- const Entry e(column, value);
+ const typename SparseMatrixEZ<number>::Entry e(column, value);
// Find position for inserting
// should return first Entry with
// higher column number.
- std::vector<Entry>::iterator col = lower_bound(begin(), end, e);
+ typename std::vector<Entry>::iterator col = lower_bound(begin(), end, e);
// Insert Entry if column did not exist.
// Edit existing entry otherwise.
template <typename number>
-SparseMatrixEZ<number>::SparseMatrixEZ(const SparseMatrixEZ&)
+SparseMatrixEZ<number>::SparseMatrixEZ(const SparseMatrixEZ<number>&)
{
Assert(false, ExcNotImplemented());
}
template <typename number>
bool
-SparseMatrixEZ<number>::empty()
+SparseMatrixEZ<number>::empty() const
{
return ((n_columns == 0) && (rows.size()==0));
}
Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
- std::vector<Row>::const_iterator row = rows.begin();
- const std::vector<Row>::const_iterator end_row = rows.end();
+ typename std::vector<Row>::const_iterator row = rows.begin();
+ const typename std::vector<Row>::const_iterator end_row = rows.end();
for (unsigned int i=0; row != end_row; ++i, ++row)
{
- std::vector<Entry>::const_iterator col = row->begin();
- const std::vector<Entry>::const_iterator end_col = row->end();
+ typename std::vector<Entry>::const_iterator col = row->begin();
+ const typename std::vector<Entry>::const_iterator end_col = row->end();
double s = 0.;
for (;col != end_col; ++col)
Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
- std::vector<Row>::const_iterator row = rows.begin();
- const std::vector<Row>::const_iterator end_row = rows.end();
+ typename std::vector<Row>::const_iterator row = rows.begin();
+ const typename std::vector<Row>::const_iterator end_row = rows.end();
for (unsigned int i=0; row != end_row; ++i, ++row)
{
- std::vector<Entry>::const_iterator col = row->begin();
- const std::vector<Entry>::const_iterator end_col = row->end();
+ typename std::vector<Entry>::const_iterator col = row->begin();
+ const typename std::vector<Entry>::const_iterator end_col = row->end();
double s = 0.;
for (;col != end_col; ++col)
Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
- std::vector<Row>::const_iterator row = rows.begin();
- const std::vector<Row>::const_iterator end_row = rows.end();
+ typename std::vector<Row>::const_iterator row = rows.begin();
+ const typename std::vector<Row>::const_iterator end_row = rows.end();
for (unsigned int i=0; row != end_row; ++i, ++row)
{
- std::vector<Entry>::const_iterator col = row->begin();
- const std::vector<Entry>::const_iterator end_col = row->end();
+ typename std::vector<Entry>::const_iterator col = row->begin();
+ const typename std::vector<Entry>::const_iterator end_col = row->end();
- double s = 0.;
for (;col != end_col; ++col)
dst(col->column) += col->value * src(i);
}
}
template <typename number>
-template <typename somenumber>
unsigned int
SparseMatrixEZ<number>::memory_consumption() const
{
unsigned int result = sizeof (*this)
+ sizeof(Row) * sizeof (rows);
- for (std::vector<Row>::const_iterator r = rows.begin();
+ for (typename std::vector<Row>::const_iterator r = rows.begin();
r != rows.end(); ++r)
result += r->size() * sizeof(Entry);
%.go : %.cc Makefile
@echo =====debug========= $<
@$(CXX) $(flags) -c $< -o $@
+
+%.o : %.cc Makefile
+ @echo =====optimized===== $<
+ @$(CXX) $(CXXFLAGS.o) -c $< -o $@
%.exe :
@echo =====linking======= $@
@$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS)
-# Generated automatically from Makefile.in by configure.
-# Generated automatically from Makefile.in by configure.
############################################################
# $Id$
-# Copyright (C) 2000, 2001 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002 by the deal.II authors
############################################################
############################################################
libraries = $(lib-lac.g) \
$(lib-base.g)
+libraries_o = $(lib-lac.o) \
+ $(lib-base.o)
default: run-tests
full_matrix.exe : full_matrix.go $(libraries)
matrix_out.exe : matrix_out.go $(libraries)
solver.exe : solver.go testmatrix.go $(libraries)
+sparse_matrices.exe : sparse_matrices.go testmatrix.go $(libraries)
sparsity_pattern.exe : sparsity_pattern.go testmatrix.go $(libraries)
sparse_ilu.exe : sparse_ilu.go testmatrix.go $(libraries)
vector-vector.exe : vector-vector.go $(libraries)
#mg.exe : mg.go testmatrix.go $(libraries)
-tests = block_matrices block_vector block_vector_iterator full_matrix \
+tests = sparse_matrices block_matrices block_vector block_vector_iterator \
+ full_matrix \
matrix_out solver sparsity_pattern sparse_ilu vector-vector
############################################################
// Make matrix
FDMatrix testproblem(size, size);
SparsityPattern structure(dim, dim, 5);
- testproblem.build_structure(structure);
+ testproblem.five_point_structure(structure);
structure.compress();
SparseMatrix<double> A(structure);
- testproblem.laplacian(A);
+ testproblem.five_point(A);
PreconditionIdentity prec_no;
PreconditionSOR<> prec_sor;
// Make matrix
FDMatrix testproblem(size, size);
SparsityPattern structure(dim, dim, 5);
- testproblem.build_structure(structure);
+ testproblem.five_point_structure(structure);
structure.compress();
SparseMatrix<double> A(structure);
- testproblem.laplacian(A);
+ testproblem.five_point(A);
for (unsigned int test=0; test<2; ++test)
--- /dev/null
+//---------------------------- testmatrix.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2002 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.
+//
+//---------------------------- testmatrix.cc ---------------------------
+
+
+#include "testmatrix.h"
+#include <base/logstream.h>
+#include <lac/sparse_matrix.h>
+#include <lac/sparse_matrix_ez.h>
+#include <lac/vector.h>
+#include <lac/solver_richardson.h>
+#include <lac/precondition.h>
+
+#include <lac/sparse_matrix_ez.templates.h>
+
+#include <fstream>
+
+template<class MATRIX>
+void
+check_vmult_quadratic(const MATRIX& A,
+ const char* prefix)
+{
+ deallog.push(prefix);
+
+ Vector<double> u(A.n());
+ Vector<double> f(A.m());
+ GrowingVectorMemory<> mem;
+
+ SolverControl control(50, 1.e-3, false);
+ SolverRichardson<> rich(control, mem, .01);
+ PreconditionIdentity prec;
+
+ u = 0.;
+ f = 1.;
+
+ try
+ {
+ rich.solve(A, u, f, prec);
+ }
+ catch (...)
+ {
+ }
+ deallog << "Transpose" << std::endl;
+ try
+ {
+ rich.Tsolve(A, u, f, prec);
+ }
+ catch (...)
+ {
+ }
+ deallog.pop();
+}
+
+int main()
+{
+ std::ofstream logfile("sparse_matrices.output");
+ logfile.setf(std::ios::fixed);
+ logfile.precision(2);
+ deallog.attach(logfile);
+
+ // Switch between regression test
+ // and benchmark
+#ifdef DEBUG
+ deallog.depth_console(0);
+ deallog.log_execution_time(true);
+ deallog.log_time_differences(true);
+ const unsigned int size = 1000;
+#else
+ const unsigned int size = 100;
+#endif
+
+ FDMatrix testproblem (size, size);
+ unsigned int dim = (size-1)*(size-1);
+
+ SparsityPattern structure(dim, dim, 5);
+ testproblem.five_point_structure(structure);
+ structure.compress();
+ SparseMatrix<double> A(structure);
+ testproblem.five_point(A);
+ check_vmult_quadratic(A, "5-SparseMatrix<double>");
+
+ SparseMatrixEZ<double> E(dim,dim);
+ testproblem.five_point(E);
+ check_vmult_quadratic(E, "5-SparseMatrixEZ<double>");
+
+ A.clear();
+ structure.reinit(dim, dim, 9);
+ testproblem.nine_point_structure(structure);
+ structure.compress();
+ A.reinit(structure);
+ check_vmult_quadratic(A, "9-SparseMatrix<double>");
+
+ E.clear();
+ E.reinit(dim,dim);
+ testproblem.nine_point(E);
+ check_vmult_quadratic(E, "9-SparseMatrixEZ<double>");
+
+}
--- /dev/null
+
+9.84 :DEAL:5-SparseMatrix<double>:Richardson::Starting
+52.82 :DEAL:5-SparseMatrix<double>:Richardson::Failure step 50
+0.00 :DEAL:5-SparseMatrix<double>::Transpose
+1.56 :DEAL:5-SparseMatrix<double>:RichardsonT::Starting
+80.04 :DEAL:5-SparseMatrix<double>:RichardsonT::Failure step 50
+0.00 :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00 :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+16.61 :DEAL:5-SparseMatrixEZ<double>:Richardson::Starting
+67.84 :DEAL:5-SparseMatrixEZ<double>:Richardson::Failure step 50
+0.00 :DEAL:5-SparseMatrixEZ<double>::Transpose
+1.52 :DEAL:5-SparseMatrixEZ<double>:RichardsonT::Starting
+78.19 :DEAL:5-SparseMatrixEZ<double>:RichardsonT::Failure step 50
+0.00 :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00 :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+13.09 :DEAL:9-SparseMatrix<double>:Richardson::Starting
+73.20 :DEAL:9-SparseMatrix<double>:Richardson::Failure step 50
+0.00 :DEAL:9-SparseMatrix<double>::Transpose
+2.35 :DEAL:9-SparseMatrix<double>:RichardsonT::Starting
+119.03:DEAL:9-SparseMatrix<double>:RichardsonT::Failure step 50
+0.00 :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00 :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+32.49 :DEAL:9-SparseMatrixEZ<double>:Richardson::Starting
+94.45 :DEAL:9-SparseMatrixEZ<double>:Richardson::Failure step 50
+0.00 :DEAL:9-SparseMatrixEZ<double>::Transpose
+2.14 :DEAL:9-SparseMatrixEZ<double>:RichardsonT::Starting
+109.61:DEAL:9-SparseMatrixEZ<double>:RichardsonT::Failure step 50
+0.00 :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00 :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
// generate usual sparsity pattern
const unsigned int N = 15;
SparsityPattern sp1((N-1)*(N-1), (N-1)*(N-1), 5);
- FDMatrix(N,N).build_structure (sp1);
+ FDMatrix(N,N).five_point_structure (sp1);
sp1.compress ();
deallog << sp1.n_rows() << " " << sp1.n_cols() << " "
<< sp1.bandwidth() << " " << sp1.n_nonzero_elements()
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2002 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include "testmatrix.h"
#include <lac/sparse_matrix.h>
+#include <lac/sparse_matrix_ez.h>
#include <lac/vector.h>
+#include <lac/sparse_matrix_ez.templates.h>
+
FDMatrix::FDMatrix(unsigned int nx, unsigned int ny)
:
nx(nx), ny(ny)
{}
void
-FDMatrix::build_structure(SparsityPattern& structure) const
+FDMatrix::five_point_structure(SparsityPattern& structure) const
{
for(unsigned int i=0;i<=ny-2;i++)
{
}
}
-template<typename number>
void
-FDMatrix::laplacian(SparseMatrix<number>& A) const
+FDMatrix::nine_point_structure(SparsityPattern& structure) const
+{
+ for(unsigned int i=0;i<=ny-2;i++)
+ {
+ for(unsigned int j=0;j<=nx-2; j++)
+ {
+ // Number of the row to be entered
+ unsigned int row = j+(nx-1)*i;
+ structure.add(row, row);
+ if (j>0)
+ {
+ structure.add(row-1, row);
+ structure.add(row, row-1);
+ if (i>0)
+ {
+ structure.add(row-1, row-(nx-1));
+ structure.add(row-(nx-1), row-1);
+ }
+ if (i<ny-2)
+ {
+ structure.add(row-1, row+(nx-1));
+ structure.add(row+(nx-1), row-1);
+ }
+ }
+ if (j<nx-2)
+ {
+ structure.add(row+1, row);
+ structure.add(row, row+1);
+ if (i>0)
+ {
+ structure.add(row+1, row-(nx-1));
+ structure.add(row-(nx-1), row+1);
+ }
+ if (i<ny-2)
+ {
+ structure.add(row+1, row+(nx-1));
+ structure.add(row+(nx-1), row+1);
+ }
+ }
+ if (i>0)
+ {
+ structure.add(row-(nx-1), row);
+ structure.add(row, row-(nx-1));
+ }
+ if (i<ny-2)
+ {
+ structure.add(row+(nx-1), row);
+ structure.add(row, row+(nx-1));
+ }
+ }
+ }
+}
+
+template<typename MATRIX>
+void
+FDMatrix::nine_point(MATRIX& A) const
+{
+ for(unsigned int i=0;i<=ny-2;i++)
+ {
+ for(unsigned int j=0;j<=nx-2; j++)
+ {
+ // Number of the row to be entered
+ unsigned int row = j+(nx-1)*i;
+
+ A.set(row, row, 20.);
+ if (j>0)
+ {
+ A.set(row-1, row, -4.);
+ A.set(row, row-1, -4.);
+ if (i>0)
+ {
+ A.set(row-1, row-(nx-1), -1.);
+ A.set(row-(nx-1), row-1, -1.);
+ }
+ if (i<ny-2)
+ {
+ A.set(row-1, row+(nx-1), -1.);
+ A.set(row+(nx-1), row-1, -1.);
+ }
+ }
+ if (j<nx-2)
+ {
+ A.set(row+1, row, -4.);
+ A.set(row, row+1, -4.);
+ if (i>0)
+ {
+ A.set(row+1, row-(nx-1), -1.);
+ A.set(row-(nx-1), row+1, -1.);
+ }
+ if (i<ny-2)
+ {
+ A.set(row+1, row+(nx-1), -1.);
+ A.set(row+(nx-1), row+1, -1.);
+ }
+ }
+ if (i>0)
+ {
+ A.set(row-(nx-1), row, -4.);
+ A.set(row, row-(nx-1), -4.);
+ }
+ if (i<ny-2)
+ {
+ A.set(row+(nx-1), row, -4.);
+ A.set(row, row+(nx-1), -4.);
+ }
+ }
+ }
+}
+
+template<typename MATRIX>
+void
+FDMatrix::five_point(MATRIX& A) const
{
for(unsigned int i=0;i<=ny-2;i++)
{
-template void FDMatrix::laplacian(SparseMatrix<long double>&) const;
-template void FDMatrix::laplacian(SparseMatrix<double>&) const;
-template void FDMatrix::laplacian(SparseMatrix<float>&) const;
+template void FDMatrix::five_point(SparseMatrix<long double>&) const;
+template void FDMatrix::five_point(SparseMatrix<double>&) const;
+template void FDMatrix::five_point(SparseMatrixEZ<double>&) const;
+template void FDMatrix::five_point(SparseMatrix<float>&) const;
+template void FDMatrix::nine_point(SparseMatrix<long double>&) const;
+template void FDMatrix::nine_point(SparseMatrix<double>&) const;
+template void FDMatrix::nine_point(SparseMatrixEZ<double>&) const;
+template void FDMatrix::nine_point(SparseMatrix<float>&) const;
template void FDMatrix::gnuplot_print(std::ostream& s, const Vector<long double>& V) const;
template void FDMatrix::gnuplot_print(std::ostream& s, const Vector<double>& V) const;
template void FDMatrix::gnuplot_print(std::ostream& s, const Vector<float>& V) const;
/**
* Generate the matrix structure.
*/
- void build_structure(SparsityPattern& structure) const;
+ void five_point_structure(SparsityPattern& structure) const;
+
+ /**
+ * Generate the matrix structure.
+ */
+ void nine_point_structure(SparsityPattern& structure) const;
/**
* Fill the matrix with values.
*/
- template <typename number>
- void laplacian(SparseMatrix<number>&) const;
+ template <typename MATRIX>
+ void five_point(MATRIX&) const;
+
+ /**
+ * Fill the matrix with values.
+ */
+ template <typename MATRIX>
+ void nine_point(MATRIX&) const;
template <typename number>
void gnuplot_print(std::ostream&, const Vector<number>&) const;