Roland Kaufmann
8c6302f123
Every program that relies on manual inspection has been moved to a new (hopefully short-lived) directory called not-unit/; every remaining file has been given the prefix test_ to indicate that this is the executable test to be run.
127 lines
3.3 KiB
C++
127 lines
3.3 KiB
C++
/*
|
|
Copyright 2012 SINTEF ICT, Applied Mathematics.
|
|
Copyright 2012 Statoil ASA.
|
|
|
|
This file is part of the Open Porous Media Project (OPM).
|
|
|
|
OPM is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OPM is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
|
*/
|
|
|
|
#include <cassert>
|
|
#include <cmath>
|
|
#include <cstddef>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <vector>
|
|
|
|
#include <boost/shared_ptr.hpp>
|
|
|
|
#include <opm/core/linalg/sparse_sys.h>
|
|
#include <opm/core/linalg/LinearSolverAGMG.hpp>
|
|
|
|
namespace {
|
|
std::size_t
|
|
compute_nnz(const std::size_t m)
|
|
{
|
|
assert (m > 0);
|
|
|
|
std::size_t nnz = m; // A(i,i)
|
|
|
|
nnz += (m > 1) ? 2 : 0; // A(0,1), A(m-1,m-2)
|
|
nnz += (m > 2) ? 2 * (m - 2) : 0; // A(i,i-1), A(i,i+1)
|
|
|
|
return nnz;
|
|
}
|
|
|
|
|
|
boost::shared_ptr<CSRMatrix>
|
|
build_laplace_1d(const std::size_t m)
|
|
{
|
|
assert (m >= 2);
|
|
|
|
const std::size_t nnz = compute_nnz(m);
|
|
|
|
boost::shared_ptr<CSRMatrix>
|
|
A(csrmatrix_new_known_nnz(m, nnz), csrmatrix_delete);
|
|
|
|
A->ia[ 0 ] = 0;
|
|
|
|
// First row
|
|
A->ia[ 0 + 1 ] = A->ia[ 0 ];
|
|
A->ja[ A->ia[0 + 1] ] = 0 + 0;
|
|
A->sa[ A->ia[0 + 1] ++ ] = 2.0;
|
|
A->ja[ A->ia[0 + 1] ] = 0 + 1;
|
|
A->sa[ A->ia[0 + 1] ++ ] = - 1.0;
|
|
|
|
// General rows
|
|
for (std::size_t i = 1; i < m - 1; ++i) {
|
|
A->ia[i + 1] = A->ia[i];
|
|
|
|
A->ja[ A->ia[i + 1] ] = int(i) - 1;
|
|
A->sa[ A->ia[i + 1] ++ ] = - 1.0;
|
|
|
|
A->ja[ A->ia[i + 1] ] = int(i) ;
|
|
A->sa[ A->ia[i + 1] ++ ] = 2.0;
|
|
|
|
A->ja[ A->ia[i + 1] ] = int(i) + 1;
|
|
A->sa[ A->ia[i + 1] ++ ] = - 1.0;
|
|
}
|
|
|
|
// Last row
|
|
A->ia[ (m - 1) + 1 ] = A->ia[ m - 1 ];
|
|
|
|
A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) - 1;
|
|
A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = - 1.0;
|
|
|
|
A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) ;
|
|
A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = 2.0;
|
|
|
|
return A;
|
|
}
|
|
}
|
|
|
|
|
|
int main()
|
|
{
|
|
const std::size_t m = 10;
|
|
|
|
boost::shared_ptr<CSRMatrix> A = build_laplace_1d(m);
|
|
|
|
// Form right-hand side [1, 0, 0, ...., 0, 1]
|
|
std::vector<double> b(m, 0.0);
|
|
b[0] = 1.0; b.back() = 1.0;
|
|
|
|
// Allocate solution vector
|
|
std::vector<double> x(m);
|
|
|
|
// Create solver for SPD system.
|
|
Opm::LinearSolverAGMG linsolve(100, 1e-9, true);
|
|
|
|
Opm::LinearSolverInterface::LinearSolverReport
|
|
rpt = linsolve.solve(A.get(), & b[0], & x[0]);
|
|
|
|
double e = 0.0;
|
|
for (std::size_t i = 0; i < m; ++i) {
|
|
const double d = x[i] - 1.0;
|
|
e += d * d;
|
|
}
|
|
|
|
std::cerr << "|| e ||_2 = "
|
|
<< std::scientific
|
|
<< std::setprecision(5)
|
|
<< std::sqrt(e) / double(m) << '\n';
|
|
|
|
return 0;
|
|
}
|