ddemidov / amgcl Goto Github PK
View Code? Open in Web Editor NEWC++ library for solving large sparse linear systems with algebraic multigrid method
Home Page: http://amgcl.readthedocs.org/
License: MIT License
C++ library for solving large sparse linear systems with algebraic multigrid method
Home Page: http://amgcl.readthedocs.org/
License: MIT License
Python test fails on my system with the following message (all other tests pass):
(i am using commit 9ac586e)
Number of levels: 3
Operator complexity: 1.33
Grid complexity: 1.18
level unknowns nonzeros
---------------------------------
0 65536 323600 (75.17%)
1 10795 95969 (22.29%)
2 1227 10897 ( 2.53%)
Iterations: 100
Error: 0.999901
F
======================================================================
FAIL: test_solver (__main__.TestPyAMGCL)
----------------------------------------------------------------------
Traceback (most recent call last):
File "./python/test_pyamgcl", line 61, in test_solver
self.assertTrue(np.linalg.norm(res) < 1e-6)
AssertionError: False is not true
----------------------------------------------------------------------
Ran 1 test in 2.720s
FAILED (failures=1)
My system is: Ubuntu 14.04 x64, Python 2.7.6, Boost 1.54, numpy 1.8.1, scipy 0.13.3, GCC 4.8.2
Hi!
How do i use opencl backend in python. Is this possible curently?
Regards,
Marek
Hi,
I'm looking to solve a linear system which is algebraically similar to the Stokes problem, in the sense that I have a matrix in the following block form:
[ A B ]
[ B^T 0 ]
with A symmetric positive definite, and B has full row rank.
I am trying to find an efficient way to solve this system, and I hope to be able to adapt techniques one may use for the Stokes problem, as seen towards the end in this lecture by Wolfgang Bangerth. In this video, Bangerth advocates using the following (right) block preconditioner in conjunction with GMRES:
[ A^-1 -A^-1 B^T S^-1 ]
[ 0 S^-1 ]
where S = B^T A^-1 B. Of course, one does not know the exact inverses of A and S, so one would like to approximate these. For the Stokes problem, one can approximate S^-1 by the inverse of the mass matrix (but for my problem I need to use a different matrix that I hope is somewhat similar), and Bangerth suggests using "one step of multigrid" to approximate the action of applying A^-1.
I've been poking around the docs and the source code, and I think it should be possible to accomplish this kind of block-preconditioning with amgcl
. Is this correct? I see that there's a block_preconditioner.hpp
related to MPI, and also a schur_pressure_correction.hpp
, both of which might be useful for me to understand how to assemble this, but I find it a little hard to decipher what's going on, especially since I'm not familiar with Schur pressure correction.
In particular, can I perform "one step" of AMG for the approximation of A^-1? I'm struggling to see how to assemble the pieces, and any suggestions (or confirmation on whether or not I can do what I want) would be most welcome!
P.S: I realize this is not truly an "issue", but I could not find any other place to ask questions to the community of users of amgcl. I hope this is OK.
Hi,
I recently downloaded AMGCL and attempted to apply it as a solver for a 3D elasticity problem using tetrahedral finite-elements. I've noticed that aggregation coarsening algorithms are very aggressive and typically reduce the number by 2 to 3 orders of magnitude between levels. I was wondering if anyone could offer any guidance to reduce this to approximately a factor of 8?
Thanks,
Nick
Hi Denis,
My research involves modelling biological processes, for example blood vessels growing into a tissue. I use the AMGCL library for solving the diffusion of molecules within the tissue domain. The blood vessels act as a dirichlet BC for the Poisson, acting as a source and a sink for different molecules, so I have a process which roughly like:
A = stiffness_matrix(n,n,n) # stiffness matrix for 3d poisson
for t < total_time:
BC = vessel_locations
A_new = applyBC(A,BC)
precond = pyamgl.amg(A_new)
solver=pyamgcl.solver(precond)
C = agmcl.solver(C)
vessel_locations = move_vessesl(C,vessel_locations)
So at each iteration the BC are changing as the vasculature network grows. I was wondering if you thought it would be feasible for me to "update" the preconditioner each iteration rather than recreating it. I currently do this for the matrix but this doesn't save as much time as if I could also modify the MG.
Cheers,
Duncan
Hi Denis,
Hope you are doing great. Nice to see that your project keeps advancing. There is a minor issue which I cannot really sort out myself. I recently tried to use your QR decomposition on a complex matrix. However, I get a few nan's in the result. When using Mathematica things seem okay. I have added some code in your QR test case to check it:
Input matrix can be found here: X.tar.gz
I don't know if the problem is with the code or with my matrix. Maybe you have time to take a look at it.
Cheers,
Chris
PS: There was a minor bug in the mm.hpp file (see my commit).
Proof-of-concept: https://gist.github.com/ddemidov/a46ff9f58514fbddee59da1d0f3a0bb5
I am trying to speed up CG for (F'F + \lambda*\eye) x = b
by using preconditioning. Matrix F
is large and sparse but computing F'F
explicitly is too costly. Would it be possible to use your software to build AMG preconditioner in this situation? All ideas are welcome.
Just mention that the problem involves solving the same system for several RHSs, i.e. it is no problem if computing the preconditioner takes time.
Thank you.
Hi @ddemidov,
i have been tracking bugs in my code, by compiling with
-std=c++11 -D_DEBUG -D_GLIBCXX_DEBUG
the checks detect a problem in subdomain_deflation when using Eptr within the MPI_Gatherv function
essentially the thing is that Eptr is empty everywhere except on the master_rank, so
&Eptr[0] is not a valid access.
std::vector<int> Eptr;
if (comm.rank == master_rank) {
Eptr.resize(dv_start[cgroup_end] - dv_start[cgroup_beg] + 1, 0);
}
//MY SUGGESTION WOULD BE TO ALWAYS RESIZE Eptr to at least 1 (it will be set to zero)
MPI_Gatherv(
&eptr[1], ndv, MPI_INT, &Eptr[0] + 1,
const_cast<int*>(&ssize[0]), const_cast<int*>(&sstart[0]),
MPI_INT, 0, slaves_comm
);
I am seeing a slightly strange behaviour here. Maybe you could take a quick look at is. I have a complex system I want to solve:
When I use GMRES with solver.maxiter set to 100, the result looks fine. However, when I use solver.maxiter with 25 I get NaN. Other solvers do not show this behaviour.
To me this seems a bit odd... Thanks for looking at it.
Chris
Hi,
I am working on solving a system of FEM equations with a few constraints, formulated as a KKT system
[ K C^T ] [ u ] = [ f ]
[ C 0 ] [ lambda ] [ 0 ]
From what I can tell, this can be done with AMG if separate smoothing/restriction operators are used for the primal variables and for the lagrange multipliers. See
http://crd.lbl.gov/assets/Uploads/ANAG/adamskktamg.pdf
Would there be an easy way to incorporate this into amgcl? I don't know the code well enough yet to know whether I can reuse the existing smoothers/restriction for the system matrix, then just add my own simple one for the constraints.
Thanks!
In general and specifically inside subdomain_deflation
Hello
I am trying to the AMGCL solver in phd thesis to solve the pressure poisson equation
In my solver algorithm, for every time step i do the following
1-construct the matrix
2-setup the preconditioner
3-solve the system with relative tol = 0.01
4-update the rhs
5-solve the system with relative tol = 0.01
6-update the rhs
7-solve the system with relative abstol = 1e-12
During the time step , the system matrix remain unchanged but the rhs changes and iterative solver parameters change
(1) Is it possible to change the iterative solver parameters after setup without having to destroy that setup and redo it again.
(2) Also would it possible to use the same precond for 2 or 3 more time steps even when the matrix system changes
I read the doc but I can't find a way to access the params after setup
I tried the following
boost::property_tree::ptree prm;
prm.put("solver.tol", 1e-2);
prm.put("precond.coarsening.relax", 0.75);
Solver solve(boost::tie(n3, ptr, col, val),prm);
boost::tie(iters, error) = solve(f, x);
solve.prm.solver.abstol = 1e-12; // change the abs tol
boost::tie(iters, error) = solve(f, x);
But the change was ignored
Thank you
Kamra
So to build I am first creating a build directory inside amgcl and then running
[ steven ] [~/Programming/amgcl/build] > cmake .. [master ✔]
-- No build type selected, default to RelWithDebInfo
-- The C compiler identification is GNU 5.3.1
-- The CXX compiler identification is GNU 5.3.1
-- Check for working C compiler: /usr/bin/cc
-- Check for working C compiler: /usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Detecting C compile features
-- Detecting C compile features - done
-- Check for working CXX compiler: /usr/bin/c++
-- Check for working CXX compiler: /usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Try OpenMP C flag = [-fopenmp]
-- Performing Test OpenMP_FLAG_DETECTED
-- Performing Test OpenMP_FLAG_DETECTED - Success
-- Try OpenMP CXX flag = [-fopenmp]
-- Performing Test OpenMP_FLAG_DETECTED
-- Performing Test OpenMP_FLAG_DETECTED - Success
-- Found OpenMP: -fopenmp
-- Boost version: 1.60.0
-- Found the following Boost libraries:
-- program_options
-- system
-- filesystem
-- serialization
-- unit_test_framework
-- Found METIS: /usr/include/metis
-- Could NOT find SCOTCH (missing: SCOTCH_LIBRARIES)
-- Could NOT find PASTIX (missing: PASTIX_INCLUDES PASTIX_LIBRARIES)
-- Looking for pthread.h
-- Looking for pthread.h - found
-- Looking for pthread_create
-- Looking for pthread_create - not found
-- Looking for pthread_create in pthreads
-- Looking for pthread_create in pthreads - not found
-- Looking for pthread_create in pthread
-- Looking for pthread_create in pthread - found
-- Found Threads: TRUE
-- Found PythonInterp: /home/steven/.anaconda2/bin/python (found version "2.7.11")
-- Found PythonLibs: /home/steven/.anaconda2/lib/libpython2.7.so (found suitable version "2.7.11", minimum required is "2.7.11")
-- Found NumPy: version "1.10.4" /home/steven/.anaconda2/lib/python2.7/site-packages/numpy/core/include
-- Boost version: 1.60.0
-- Found the following Boost libraries:
-- system
-- filesystem
-- python
-- Configuring done
-- Generating done
-- Build files have been written to: /home/steven/Programming/amgcl/build
The first thing I notice is that SCOTCH and PASTIX are missing. I have tried installing scotch libraries. I have scotch, libscotchmetis-dev, libscotch-dev
installed. The other thing I notice is that it is calling my boost version 1.60. Running $ dpkg -S /usr/include/boost/version.hpp
I see that the output is libboost1.58-dev:amd64
. Doing $ find ~/.anaconda2 -name boost
(my anaconda install path, see cmake output) I find it listing packages as the 1.60.0 version, as would be suggested by the cmake output.
Since no errors were thrown I decided to try and run make
where it failed at [ 3%] Linking CXX executable test_solver_ns_eigen
on a boost error. Where it couldn't find a library. It is looking in the anaconda/include/boost
folder so I thought I'd reset to the system's version of boost.
$ cmake -DBOOST_ROOT=/usr/include/boost/ -DBOOST_INCLUDEDIR=/usr/include ..
which does show the boost version as 1.58.0
at the end of the cmake. Again it seems to be finding the anaconda version of boost and erroring at the same 3% mark. Here is the error
[ 3%] Linking CXX executable test_solver_ns_eigen
CMakeFiles/test_solver_ns_eigen.dir/test_solver_ns_eigen.cpp.o: In function `__static_initialization_and_destruction_0':
/home/steven/Programming/amgcl/tests/test_solver_ns_eigen.cpp:8: undefined reference to `boost::unit_test::ut_detail::auto_test_unit_registrar::auto_test_unit_registrar(boost::unit_test::basic_cstring<char const>, boost::unit_test::basic_cstring<char const>, unsigned long, boost::unit_test::decorator::collector&)'
CMakeFiles/test_solver_ns_eigen.dir/test_solver_ns_eigen.cpp.o: In function `boost::unit_test::make_test_case(boost::function<void ()> const&, boost::unit_test::basic_cstring<char const>, boost::unit_test::basic_cstring<char const>, unsigned long)':
/home/steven/.anaconda2/include/boost/test/tree/test_unit.hpp:249: undefined reference to `boost::unit_test::test_case::test_case(boost::unit_test::basic_cstring<char const>, boost::unit_test::basic_cstring<char const>, unsigned long, boost::function<void ()> const&)'
CMakeFiles/test_solver_ns_eigen.dir/test_solver_ns_eigen.cpp.o: In function `__static_initialization_and_destruction_0':
/home/steven/Programming/amgcl/tests/test_solver_ns_eigen.cpp:10: undefined reference to `boost::unit_test::ut_detail::auto_test_unit_registrar::auto_test_unit_registrar(boost::unit_test::test_case*, boost::unit_test::decorator::collector&, unsigned long)'
CMakeFiles/test_solver_ns_eigen.dir/test_solver_ns_eigen.cpp.o: In function `bool boost::test_tools::tt_detail::check_frwd<boost::test_tools::check_is_small_t, double, double>(boost::test_tools::check_is_small_t, boost::unit_test::lazy_ostream const&, boost::unit_test::basic_cstring<char const>, unsigned long, boost::test_tools::tt_detail::tool_level, boost::test_tools::tt_detail::check_type, double const&, char const*, double const&, char const*)':
/home/steven/.anaconda2/include/boost/test/tools/old/impl.hpp:92: undefined reference to `boost::test_tools::tt_detail::report_assertion(boost::test_tools::assertion_result const&, boost::unit_test::lazy_ostream const&, boost::unit_test::basic_cstring<char const>, unsigned long, boost::test_tools::tt_detail::tool_level, boost::test_tools::tt_detail::check_type, unsigned long, ...)'
/home/steven/.anaconda2/include/boost/test/tools/old/impl.hpp:92: undefined reference to `boost::test_tools::tt_detail::report_assertion(boost::test_tools::assertion_result const&, boost::unit_test::lazy_ostream const&, boost::unit_test::basic_cstring<char const>, unsigned long, boost::test_tools::tt_detail::tool_level, boost::test_tools::tt_detail::check_type, unsigned long, ...)'
collect2: error: ld returned 1 exit status
tests/CMakeFiles/test_solver_ns_eigen.dir/build.make:99: recipe for target 'tests/test_solver_ns_eigen' failed
make[2]: *** [tests/test_solver_ns_eigen] Error 1
CMakeFiles/Makefile2:105: recipe for target 'tests/CMakeFiles/test_solver_ns_eigen.dir/all' failed
make[1]: *** [tests/CMakeFiles/test_solver_ns_eigen.dir/all] Error 2
Makefile:138: recipe for target 'all' failed
make: *** [all] Error 2
Turning the debug flag for boost on (in cmake) I see no paths leading to anaconda. cmake -DBOOST_ROOT:PATHNAME=/usr/include/boost/ -DBOOST_INCLUDEDIR=/usr/include -DBoost_DEBUG=ON .. | grep ".anaconda2"
shows no output, as expected.
And at this point I'm pretty lost at where the error is arising.
Hi Denis,
I really like your library. I have started using it in place of pyamg due to the higher performance. I had to modify your code a little bit. Updating it to use the latest pybind11 and also I noticed that the pyamgcl did not utilise more than a single core on my PC. I just added the pybind11 commands to acquire and later release the GIL. Was it intentional to use just a single core?
Also if you are interested I can branch and push these small modifications.
Cheers,
Duncan
Hi @ddemidov i tried to make a PR but the system does not allow.
i am attaching here the file with minimal modifications (using numa_vector instead of small_vector)
i should say that i have the feeling that not everything is fixed, since the result still depends on the number of processors, however the PR is still valid
I have problems getting the solver_cuda project to compile, I am using cuda 9.1 and VS 2015, seems to be something with nvcc not handling c++ files correctly, error output:
1>------ Build started: Project: solver_cuda, Configuration: Release x64 ------
1> Building NVCC (Device) object examples/CMakeFiles/solver_cuda.dir/Release/solver_cuda_generated_solver.cu.obj
1> -- Removing D:/eas/gits/amgcl/build_x64/examples/CMakeFiles/solver_cuda.dir//Release/solver_cuda_generated_solver.cu.obj
1> "C:/Program Files/CMake/bin/cmake.exe" -E remove D:/eas/gits/amgcl/build_x64/examples/CMakeFiles/solver_cuda.dir//Release/solver_cuda_generated_solver.cu.obj
1> -- Generating dependency file: D:/eas/gits/amgcl/build_x64/examples/CMakeFiles/solver_cuda.dir//solver_cuda_generated_solver.cu.obj.NVCC-depend
1> "C:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v9.1/bin/nvcc.exe" -M -D__CUDACC__ D:/eas/gits/amgcl/build_x64/examples/solver.cu -o D:/eas/gits/amgcl/build_x64/examples/CMakeFiles/solver_cuda.dir//solver_cuda_generated_solver.cu.obj.NVCC-depend -ccbin "C:/Program Files (x86)/Microsoft Visual Studio 14.0/VC/bin" -m64 -DAMGCL_BLOCK_SIZES=(3)(4) -DNOMINMAX -D_VARIADIC_MAX=10 -D_SCL_SECURE_NO_WARNINGS -DAMGCL_ASYNC_SETUP -Xcompiler ,\"/DWIN32\",\"/D_WINDOWS\",\"/W3\",\"/GR\",\"/EHsc\",\"/MD\",\"/O2\",\"/Ob2\",\"/DNDEBUG\" -Xcompiler -openmp -DSOLVER_BACKEND_CUDA -DNVCC "-IC:/Program Files/NVIDIA GPU Computing Toolkit/CUDA/v9.1/include" -ID:/eas/gits/amgcl/examples/../tests -ID:/eas/boost_1_67_0 -ID:/eas/gits/amgcl
1> solver.cu
1>D:/eas/boost_1_67_0\boost/preprocessor/slot/detail/shared.hpp(27): fatal error C1012: unmatched parenthesis: missing ')'
1> CMake Error at solver_cuda_generated_solver.cu.obj.Release.cmake:219 (message):
1> Error generating
1> D:/eas/gits/amgcl/build_x64/examples/CMakeFiles/solver_cuda.dir//Release/solver_cuda_generated_solver.cu.obj
1>
Any clues to a fix to this problem? The ordinary solver c++ project compiles and build just fine!
So I am reading through the benchmark distributed Poisson equation implementation to figure out how I could have the linear elasticity problem I have with amgcl run on several MPI nodes. But I still have some questions I was hoping you might answer:
Allow dense blocks with overloaded arithmetic operations (e.g. Eigen matrices) to be used as value_type.
Hi Denis
Can you please provide an example of how to wrap a ViennaCL/Eigen sparse matrix to Epetra and vice-versa so that we can use amgcl or use Anasazi that is part of the trilinos project?
I may have found a small bug in the generic GMRES implementation. Please take a look at [1] in Section 6.5.9 on the complex implementation of GMRES. They state that the complex Givens rotation is defined slightly different from the real one, see equation 6.80. When I change it here accordingly, I get much better convergence in the complex case. Can you confirm this? Do you think there will be a similar issue for LGMRES and FGMRES?
Chris
[1] Saad, Yousef. Iterative methods for sparse linear systems. Siam, 2003.
Hi.
I would like to solve a complex-valued problem without rewriting it into a real-valued one. Unfortunately, I wasn't successful in using std::complex instead of double. I get a bunch of errors, for a large part related to implicit conversion where you wrote "1/x" instead of "1./x". Also, some comparisons have to be rewritten as "fabs(eps)>fabs(x)". Both can be fixed straightforwardly.
However, I was wondering whether you see more substantial issues within the algorithms when plugging in a complex type?
Thanks,
Chris
Looking at the code for the Gauss-Seidel relaxation, l.87 of the method apply
reads as:
for(size_t i = n; i-- > 0;) iteration_body(A, rhs, x, i);
Which looks a bit suspicious. I reckon that the iteration should start at n-1 (last row), and the index should be a signed integer (otherwise it will produce an infinite loop).
Am I missing something?
Hi,
Have you considered implementing a parallel version of the Gauss-Seidel relaxation scheme? (so called red-black / multicolor Gauss-Seidel algorithms.) Given the variety of backends I'm not sure this would be easy to do, but I thought I might ask =)
Greetings,
I am trying to compile AMGCL with support to VexCl and I am getting the following error:
/home/vicente/src/vexcl/vexcl/function.hpp:338:1: error: macro "isgreater" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isgreater )
^
/home/vicente/src/vexcl/vexcl/function.hpp:339:1: error: macro "isgreaterequal" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isgreaterequal )
^
/home/vicente/src/vexcl/vexcl/function.hpp:341:1: error: macro "isless" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isless )
^
/home/vicente/src/vexcl/vexcl/function.hpp:342:1: error: macro "islessequal" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, islessequal )
^
/home/vicente/src/vexcl/vexcl/function.hpp:343:1: error: macro "islessgreater" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, islessgreater )
^
/home/vicente/src/vexcl/vexcl/function.hpp:348:1: error: macro "isunordered" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isunordered )
^
In file included from /home/vicente/src/vexcl/vexcl/sparse/ell.hpp:52:0,
from /home/vicente/src/vexcl/vexcl/sparse/matrix.hpp:4,
from /home/vicente/src/vexcl/vexcl/vexcl.hpp:50,
from /home/vicente/src/amgcl/amgcl/backend/vexcl.hpp:40,
from /home/vicente/src/amgcl/examples/vexcl.cpp:5:
/home/vicente/src/vexcl/vexcl/function.hpp:338:1: error: macro "isgreater" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isgreater )
^
/home/vicente/src/vexcl/vexcl/function.hpp:339:1: error: macro "isgreaterequal" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isgreaterequal )
^
/home/vicente/src/vexcl/vexcl/function.hpp:341:1: error: macro "isless" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isless )
^
/home/vicente/src/vexcl/vexcl/function.hpp:342:1: error: macro "islessequal" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, islessequal )
^
/home/vicente/src/vexcl/vexcl/function.hpp:343:1: error: macro "islessgreater" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, islessgreater )
^
/home/vicente/src/vexcl/vexcl/function.hpp:348:1: error: macro "isunordered" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isunordered )
^
In file included from /home/vicente/src/vexcl/vexcl/sparse/ell.hpp:52:0,
from /home/vicente/src/vexcl/vexcl/sparse/matrix.hpp:4,
from /home/vicente/src/vexcl/vexcl/vexcl.hpp:50,
from /home/vicente/src/amgcl/amgcl/backend/vexcl.hpp:40,
from /home/vicente/src/amgcl/examples/solver.cpp:15:
/home/vicente/src/vexcl/vexcl/function.hpp:338:1: error: macro "isgreater" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isgreater )
^
/home/vicente/src/vexcl/vexcl/function.hpp:339:1: error: macro "isgreaterequal" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isgreaterequal )
^
/home/vicente/src/vexcl/vexcl/function.hpp:341:1: error: macro "isless" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isless )
^
/home/vicente/src/vexcl/vexcl/function.hpp:342:1: error: macro "islessequal" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, islessequal )
^
/home/vicente/src/vexcl/vexcl/function.hpp:343:1: error: macro "islessgreater" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, islessgreater )
^
/home/vicente/src/vexcl/vexcl/function.hpp:348:1: error: macro "isunordered" requires 2 arguments, but only 1 given
VEX_BUILTIN_FUNCTION( 2, isunordered )
It continues, but is more of the same error.
I am using boost 1.66 and compiling with c++11 support.
Best regards and thanks in advance
O. Bröker and M. J. Grote. Sparse approximate inverse smoothers for geometric and algebraic multigrid.
Dear Denis
Thanks for your work on amgcl. I'm getting a compile-time error when using the zero-copy backend, and I'm not sure where the problem is coming from. I simply took the Poisson demo straight from the docs, and changed:
Solver solve( boost::tie(n, ptr, col, val) );
to
Solver solve( amgcl::adapter::zero_copy(n, &ptr[0], &col[0], &val[0]) );
which gives me:
In file included from test_amgcl.cpp:7:
../ownCloud/stencils/inst/include/amgcl/adapter/zero_copy.hpp:67:9: error:
implicit instantiation of undefined template
'boost::STATIC_ASSERTION_FAILURE<false>'
BOOST_STATIC_ASSERT(sizeof(Ptr) == sizeof(ptr_type));
^
/usr/include/boost/static_assert.hpp:169:13: note: expanded from macro
'BOOST_STATIC_ASSERT'
sizeof(::boost::STATIC_ASSERTION_FAILURE< BOOST_STATIC_ASSER...
^
/usr/include/boost/smart_ptr/make_shared_object.hpp:819:17: note: in
instantiation of function template specialization
'amgcl::adapter::zero_copy_adapter<double>::zero_copy_adapter<int, int>'
requested here
::new( pv ) T( a1, a2, a3, a4 );
^
../ownCloud/stencils/inst/include/amgcl/adapter/zero_copy.hpp:85:19: note: in
instantiation of function template specialization
'boost::make_shared<amgcl::adapter::zero_copy_adapter<double>, unsigned
long, const int *, const int *, const double *>' requested here
return boost::make_shared< zero_copy_adapter<Val> >(n, ptr, col, val);
^
test_amgcl.cpp:85:32: note: in instantiation of function template specialization
'amgcl::adapter::zero_copy<int, int, double>' requested here
Solver solve(amgcl::adapter::zero_copy(n, ptr.data(), col.data(), val....
^
/usr/include/boost/static_assert.hpp:87:26: note: template is declared here
template <bool x> struct STATIC_ASSERTION_FAILURE;
^
In file included from test_amgcl.cpp:7:
../ownCloud/stencils/inst/include/amgcl/adapter/zero_copy.hpp:68:9: error:
implicit instantiation of undefined template
'boost::STATIC_ASSERTION_FAILURE<false>'
BOOST_STATIC_ASSERT(sizeof(Col) == sizeof(col_type));
^
/usr/include/boost/static_assert.hpp:169:13: note: expanded from macro
'BOOST_STATIC_ASSERT'
sizeof(::boost::STATIC_ASSERTION_FAILURE< BOOST_STATIC_ASSER...
^
/usr/include/boost/static_assert.hpp:87:26: note: template is declared here
template <bool x> struct STATIC_ASSERTION_FAILURE;
Is some extra casting required?
Hi,
When I build program with CUDA, it shows a error below:
1>C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v9.1\include\thrust/system/cuda/detail/internal/copy_cross_system.h(251): error C2664: 'OutputIt thrust::cuda_cub::__copy::cross_system_copy_n<thrust::cuda_cub::tag,Derived,InputIt,Size,OutputIt>(thrust::execution_policy<thrust::cuda_cub::tag> &,thrust::execution_policy<Derived> &,InputIt,Size,OutputIt,thrust::detail::true_type)': cannot convert argument 6 from 'thrust::detail::integral_constant<bool,false>' to 'thrust::detail::true_type'
I'm using solver.cpp and everything is OK without CUDA. I have tried 1.1.0 and the master, all get this error. Do I need to add more settings to run with GPU?
Thanks & Regards,
John
Hi,
I have encountered a few problems while attempting to solve complex systems of equations with a nested value_type, i.e., static_matrixstd::complex<double,int N,int N>. The modifications I needed to perform in static_matrix.hpp are below:
// struct scalar_of< static_matrix<T, N, M> >
typedef typename scalar_of<T>::type type;
// struct adjoint_impl< static_matrix<T, N, M> >
y(j,i) = adjoint(x(i,j));
// struct inner_product_impl< static_matrix<T, N, 1> >
sum += x(i) * adjoint(y(i));
// struct inner_product_impl< static_matrix<T, N, M> >
sum += x(k,i) * adjoint(y(k,j));
// struct norm_impl< static_matrix<T, N, M> >
static typename scalar_of<T>::type
get(const static_matrix<T, N, M> &x) {
T s = math::zero<T>();
for(int i = 0; i < N * M; ++i)
s += x(i) * adjoint(x(i));
return sqrt(std::abs(s));
}
// struct inverse_impl< static_matrix<T, N, N> >
T d = 1. / A(k,k);
Additionally, the compiler was attempting to perform division of an integer by complex in some other parts of the amgcl library due to expressions such as T d = 1 / A(k,k). I added the following work around to my local file.
namespace std
{
inline complex<double>
operator/(const int& x, const complex<double>& y)
{
return static_cast<double>(x) / y;
}
}
I have not tested this solution extensively, but have not encountered further issues. Is it possible to make similar changes in amgcl?
Is there any way of solving singular problems with amgcl? For example, a simple FD approximation of the Poisson problem with pure Neumann bc reads:
[-1 1 .... ]
[1 -2 1 .... ]
[. 1 -2 1 ...]
...
[ ..... 1 -1]
and has a null space (functions with constant values). The standard CG algorithm converges on this problem when the rhs has zero mean. amgcl (with smoothed aggregation, spai0 relaxation) fails with:
Zero sum in skyline_lu factorization
which I'm guessing means the coarse system is singular. Do you know of a way of working around the problem?
Dear Denis,
I've been using amgcl for a small (but fairly important) part of my master's thesis. I was already using Eigen for most of the other numerical work, so using your library with the Eigen backend turned out to be a good deal simpler than integrating with huge libraries like ML or Hypre, as I could cobble something workable together within a very short amount of time. It has certainly worked very well for me, and I greatly appreciate your work, as well as the fact that you've made it open source.
I would like to cite your work, somehow. I couldn't find anything in the docs or in the repo. Do you have a preference for how you would want potential users to cite your work on amgcl in academic papers?
Hi,
I see from the documentation that it is now possible to nest solvers in AmgCL. I was wondering if it is possible to define such nested solvers via the runtime interface? This would be great to allow quick prototyping without resorting to the python bindings.
Dear Denis,
Recently we have been having some problems while compiling AMGCL with Kratos using Visual studio 2015 as it complains that it runs out of memory. We have managed to bypass it by commenting https://github.com/ddemidov/amgcl/blob/master/amgcl/util.hpp#L126.
So our current code looks like:
struct empty_params {
empty_params() {}
empty_params(const boost::property_tree::ptree &p) {
// AMGCL_PARAMS_CHECK(p, BOOST_PP_EMPTY());
}
void get(boost::property_tree::ptree&, const std::string&) const {}
};
But @RiccardoRossi found that this was actually added as a fix for visual studio in 040b9f5.
We are not sure of what this does or if it is safe to comment at all. Any clue on the topic?
hi, i just achived to get the amgcl examples working. I am not sure which libraries are actually used (during setup i configured to path to the system atlas libraries and i also downloaded and installed the blaze-lib).
the testcase that i use is a large non-symmetric matrix-market matrix of size sizes: 631968 x 631968 (nz=23290535). Unfortunately using the solve_mm example reports a solver time of about 11s, whereas our original Incomplete LU BiConjugate Gradient (from the SLATEC library) reports only 0.4s. The results of the ILU solver are not as accurate, but that is no problem since it is only used as preconditioner.
So my question is, what are necessary steps to get the maximum speed out of AMGCL? and how can I check that the correct libraries are used? and whether the parallel execution is working properly? Are there some related benchmark results which compares different libraries, CPU vs GPU, ...?
and secondly, are there other parameters of the multigrid method, which could lead to a better performance? e.g. how can is reduce the accuracy of the solver when it is use as preconditioner? is it perhaps also possible to only solve up to the first sublevel (this would tremendously reduce the number of degrees of freedom)?
thanks for any advice!
greetings
Florian
Dear Denis,
I was looking through the code, because wanted to known, how I could use the AMG as standalone because I need some data about its convergence for my bachelor's thesis. And then I noticed by chance that it should be npost instead of npre in the for-loop in line 514 of amg.hpp. If I'm not mistaken.
By the way great work.
Regards Cedric
Hi @ddemidov
We are receiving this error:
fatal error c1060, the compiler is out of heap space
While trying to compile amgcl alongside with our code (Kratos) in windows using Visual Studio. We tried several methods to avoid this which include among others:
#define AMGCL_RUNTIME_DISABLE_SPAI1
#define AMGCL_RUNTIME_DISABLE_CHEBYSHEV
The error seems to happen always in Win7 and sometimes in Win10 with a machine of similar specs. Visual studio versions used were 2015 update 2 and 2017.
Also let me ping @RiccardoRossi.
Hi,
I was wondering if it is currently possible to extend AMGCL's classes with user-defined schemes. For example if I want to add a custom relaxation/coarsening operator (to do geometric multigrid), or my own smoothing operator, while at the same time keeping the flexibility of the runtime::make_solver
class (to experiment different strategies without recompiling the program).
If it's already possible to do that, could you provide some documentation / example ?
Just a trivial comment. There's an unused parameter in skyline_lu.hpp:207 that generates a long, and hence nasty compiler warning:
[...]
../amgcl/amgcl/solver/skyline_lu.hpp:207:9: warning: unused parameter ‘prm’ [-Wunused-parameter]
skyline_lu(const Matrix &A, const params &prm = params())
Hi again.
This time I have an OpenMP question. I would like to call your code from within a loop which can be trivially parallelized. And I would like to avoid additional parallelization within amgcl. Is it sufficient to put a #pragma omp parallel for
in front of my loop? Or should I remove the omp
statements from amgcl? If this is the case, what do you think about adding a variable for switching off OpenMP in amgcl?
https://msdn.microsoft.com/en-us/library/5187hzke.aspx
Christoph
First of all thank you for developing amgcl. Without it, I would've probably never started this little project.
I am trying to use amgcl to solve the linear elasticity problem using hexahedral elements. I have been following the documentation and the examples in the example folder and roughly speaking I do the following
#include <amgcl/amg.hpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/adapter/ublas.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/solver/bicgstab.hpp>
/*
... Some matrix assembly
*/
typedef amgcl::make_solver<
amgcl::amg<
amgcl::backend::builtin<double>,
amgcl::coarsening::smoothed_aggregation,
amgcl::relaxation::spai0
>,
amgcl::solver::bicgstab<
amgcl::backend::builtin<double>
>
> amgcl_solver;
amgcl_solver my_solver(amgcl::backend::map(my_matrix));
my_solver.precond();
vector x(matrix.size1());
for(size_t i=0; i<x.size(); ++i)
x(i) = 0.0;
size_t iters;
double resid;
boost::tie(iters, resid) = my_solver(b_vector, x);
The matrix is a boost compressed matrix. The vector a boost vector. The system I try to solve is for a single element (24 dof) and I computed by hand the correct values I expect for the result. I use np.linalg.solve
from python to check that the matrix is assembled correctly and python does indeed, when I copied the printed matrix and vector from my C++ code to python yield the correct results. But my C++ does not yield the correct result.
Here is the matrix:
((3.49524,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),(0,3.49524,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),(0,0,3.49524,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),(0.349524,0.611667,-0.174762,3.49524,1.22333,-1.22333,-0.524286,0.087381,-0.087381,0.349524,0.174762,-0.611667,-1.04857,-0.087381,1.22333,-0.699048,-0.174762,0.174762,-0.87381,-0.611667,0.611667,-1.04857,-1.22333,0.087381),(0.611667,0.349524,-0.174762,1.22333,3.49524,-1.22333,-0.087381,-1.04857,1.22333,-0.174762,-0.699048,0.174762,0.087381,-0.524286,-0.087381,0.174762,0.349524,-0.611667,-0.611667,-0.87381,0.611667,-1.22333,-1.04857,0.087381),(0,0,0,0,0,3.49524,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),(0.349524,-0.174762,0.611667,-0.524286,-0.087381,0.087381,3.49524,-1.22333,1.22333,0.349524,-0.611667,0.174762,-1.04857,1.22333,-0.087381,-0.87381,0.611667,-0.611667,-0.699048,0.174762,-0.174762,-1.04857,0.087381,-1.22333),(0.174762,-0.699048,0.174762,0.087381,-1.04857,1.22333,-1.22333,3.49524,-1.22333,-0.611667,0.349524,-0.174762,1.22333,-1.04857,0.087381,0.611667,-0.87381,0.611667,-0.174762,0.349524,-0.611667,-0.087381,-0.524286,-0.087381),(0,0,0,0,0,0,0,0,3.49524,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),(-0.524286,-0.087381,-0.087381,0.349524,-0.174762,-0.611667,0.349524,-0.611667,-0.174762,3.49524,-1.22333,-1.22333,-0.87381,0.611667,0.611667,-1.04857,1.22333,0.087381,-1.04857,0.087381,1.22333,-0.699048,0.174762,0.174762),(0.087381,-1.04857,-1.22333,0.174762,-0.699048,-0.174762,-0.611667,0.349524,0.174762,-1.22333,3.49524,1.22333,0.611667,-0.87381,-0.611667,1.22333,-1.04857,-0.087381,-0.087381,-0.524286,0.087381,-0.174762,0.349524,0.611667),(0,0,0,0,0,0,0,0,0,0,0,3.49524,0,0,0,0,0,0,0,0,0,0,0,0),(-0.699048,0.174762,0.174762,-1.04857,0.087381,1.22333,-1.04857,1.22333,0.087381,-0.87381,0.611667,0.611667,3.49524,-1.22333,-1.22333,0.349524,-0.611667,-0.174762,0.349524,-0.174762,-0.611667,-0.524286,-0.087381,-0.087381),(-0.174762,0.349524,0.611667,-0.087381,-0.524286,0.087381,1.22333,-1.04857,-0.087381,0.611667,-0.87381,-0.611667,-1.22333,3.49524,1.22333,-0.611667,0.349524,0.174762,0.174762,-0.699048,-0.174762,0.087381,-1.04857,-1.22333),(0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.49524,0,0,0,0,0,0,0,0,0),(-1.04857,0.087381,-1.22333,-0.699048,0.174762,-0.174762,-0.87381,0.611667,-0.611667,-1.04857,1.22333,-0.087381,0.349524,-0.611667,0.174762,3.49524,-1.22333,1.22333,-0.524286,-0.087381,0.087381,0.349524,-0.174762,0.611667),(-0.087381,-0.524286,-0.087381,-0.174762,0.349524,-0.611667,0.611667,-0.87381,0.611667,1.22333,-1.04857,0.087381,-0.611667,0.349524,-0.174762,-1.22333,3.49524,-1.22333,0.087381,-1.04857,1.22333,0.174762,-0.699048,0.174762),(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.49524,0,0,0,0,0,0),(-1.04857,-1.22333,0.087381,-0.87381,-0.611667,0.611667,-0.699048,-0.174762,0.174762,-1.04857,-0.087381,1.22333,0.349524,0.174762,-0.611667,-0.524286,0.087381,-0.087381,3.49524,1.22333,-1.22333,0.349524,0.611667,-0.174762),(-1.22333,-1.04857,0.087381,-0.611667,-0.87381,0.611667,0.174762,0.349524,-0.611667,0.087381,-0.524286,-0.087381,-0.174762,-0.699048,0.174762,-0.087381,-1.04857,1.22333,1.22333,3.49524,-1.22333,0.611667,0.349524,-0.174762),(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.49524,0,0,0),(-0.87381,-0.611667,-0.611667,-1.04857,-1.22333,-0.087381,-1.04857,-0.087381,-1.22333,-0.699048,-0.174762,-0.174762,-0.524286,0.087381,0.087381,0.349524,0.174762,0.611667,0.349524,0.611667,0.174762,3.49524,1.22333,1.22333),(-0.611667,-0.87381,-0.611667,-1.22333,-1.04857,-0.087381,0.087381,-0.524286,0.087381,0.174762,0.349524,0.611667,-0.087381,-1.04857,-1.22333,-0.174762,-0.699048,-0.174762,0.611667,0.349524,0.174762,1.22333,3.49524,1.22333),(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.49524))
And the rhs:
(-0.0269257,-0.0269257,0,0,0,-0.179505,0,0,0,0,0,-0.179505,0,0,0,0,0,-0.179505,0,0,0,0,0,-0.179505)
And my C++ code yields as a result for x:
(-0.00770355,-0.00770355,0,-0.00770355,-0.00770355,-0.051357,-0.193071,0.00770355,0,-0.193071,0.00770355,-0.051357,0.00770355,0.177664,0,0.00770355,0.177664,-0.051357,-0.177664,0.193071,0,-0.177664,0.193071,-0.051357)
But the correct result is:
Python result (correct and the same as the pen and paper calculation):
(-0.00770353, -0.00770353, 0. , -0.00770353, -0.00770353, -0.05135699, -0.00770353, 0.00770354, 0. , -0.00770353, 0.00770354, -0.05135699, 0.00770354, -0.00770353, 0. , 0.00770354, -0.00770353, -0.05135699, 0.00770354, 0.00770354, 0. , 0.00770354, 0.00770354, -0.05135699)
So it seems that my C++ setup does compute most entries of the result vector correctly but fails for others? If it wasn't for python I would've assumed that parts of my matrix are incorrect, but this does not seem to be the case.
The C++ solver also outputs the residuum and the number of iterations:
Simulation performed 1 number of iterations and has a final residuum of 1.02456e-15
Do you have any idea if this could be related to:
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.