Skip to content

I
INMOST
Commit 330cef96 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files
Options

Sync changes 

Synchronizing various changes
parent 0fde7a26
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 4226 additions and 3150 deletions
.gitignore
View file @ 330cef96
.svn
tests/solver_test001/matrices
solver_mtiluc2.cpp
solver_mtiluc2.hpp
CMakeLists.txt
View file @ 330cef96
......@@ -22,7 +22,8 @@ set(SOURCE solver.cpp
modify.cpp
earray.cpp
comparator.cpp
autodiff.cpp)
autodiff.cpp
solver_ddpqiluc2.cpp)
set(HEADER inmost.h
inmost_options_cmake.h
inmost_common.h
......@@ -42,6 +43,15 @@ set(INMOST_MINOR_VERSION 1)
set(INMOST_PATCH_VERSION 0)
set(INMOST_VERSION "${INMOST_MAJOR_VERSION}.${INMOST_MINOR_VERSION}.${INMOST_PATCH_VERSION}")
if( EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_mtiluc2.hpp" AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/solver_mtiluc2.cpp" )
add_definitions(-DHAVE_SOLVER_MPTILUC2)
set(HAVE_SOLVER_MPTILUC2 TRUE)
list(APPEND HEADER solver_mtiluc2.hpp)
list(APPEND SOURCE solver_mtiluc2.cpp)
else()
set(HAVE_SOLVER_MPTILUC2 FALSE)
endif()
add_library(inmost STATIC ${SOURCE} ${HEADER})
......@@ -62,6 +72,7 @@ option(COMPILE_TESTS "Compile some tests" OFF)
option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" OFF)
option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" OFF)
option(USE_SOLVER_METIS "Use METIS for matrix reordering" OFF)
option(USE_SOLVER_PETSC "Use PETSc solvers" OFF)
option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF)
#option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF)
......@@ -123,6 +134,17 @@ if(USE_PARTITIONER_PARMETIS)
endif()
endif()
if(USE_SOLVER_METIS)
find_package(METIS)
if(NOT METIS_FOUND)
set(USE_SOLVER_METIS OFF)
message("METIS NOT FOUND")
else()
include_directories(${METIS_INCLUDE_DIR})
set(USE_SOLVER_METIS ON)
message("METIS FOUND")
endif()
endif()
if(USE_PARTITIONER_ZOLTAN)
find_package(ZOLTAN)
......
Documentation/Doxyfile_pages
View file @ 330cef96
......@@ -1027,7 +1027,7 @@ HTML_FILE_EXTENSION = .html
# of the possible markers and block names see the documentation.
# This tag requires that the tag GENERATE_HTML is set to YES.
HTML_HEADER =
HTML_HEADER = header.html
# The HTML_FOOTER tag can be used to specify a user-defined HTML footer for each
# generated HTML page. If the tag is left blank doxygen will generate a standard
......@@ -1037,7 +1037,7 @@ HTML_HEADER =
# that doxygen normally uses.
# This tag requires that the tag GENERATE_HTML is set to YES.
HTML_FOOTER =
HTML_FOOTER = footer.html
# The HTML_STYLESHEET tag can be used to specify a user-defined cascading style
# sheet that is used by each HTML page. It can be used to fine-tune the look of
......
Documentation/README
View file @ 330cef96
Please use Doxygen (http://www.doxygen.org/) to generate html or pdf documentation.
Online documentation is available at http://dodo.inm.ras.ru/inmost/docs/latest/html/index.html
Online documentation is available at http://doxy.inmost.org
Documentation/footer.html 0 → 100644
View file @ 330cef96
<!-- HTML footer for doxygen 1.8.5-->
<!-- start footer part -->
<!--BEGIN GENERATE_TREEVIEW-->
<div id="nav-path" class="navpath"><!-- id is needed for treeview function! -->
<ul>
$navpath
<li class="footer">$generatedby
<a href="http://www.doxygen.org/index.html">
<img class="footer" src="$relpath^doxygen.png" alt="doxygen"/></a> $doxygenversion </li>
</ul>
</div>
<!--END GENERATE_TREEVIEW-->
<!--BEGIN !GENERATE_TREEVIEW-->
<hr class="footer"/><address class="footer"><small>
$generatedby &#160;<a href="http://www.doxygen.org/index.html">
<img class="footer" src="$relpath^doxygen.png" alt="doxygen"/>
</a> $doxygenversion
</small></address>
<!--END !GENERATE_TREEVIEW-->
</body>
</html>
Documentation/header.html 0 → 100644
View file @ 330cef96
<!-- HTML header for doxygen 1.8.5-->
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<script>
(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
})(window,document,'script','//www.google-analytics.com/analytics.js','ga');
ga('create', 'UA-59408561-3', 'auto');
ga('send', 'pageview');
</script>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen $doxygenversion"/>
<!--BEGIN PROJECT_NAME--><title>$projectname: $title</title><!--END PROJECT_NAME-->
<!--BEGIN !PROJECT_NAME--><title>$title</title><!--END !PROJECT_NAME-->
<link href="$relpath^tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="$relpath^jquery.js"></script>
<script type="text/javascript" src="$relpath^dynsections.js"></script>
$treeview
$search
$mathjax
<link href="$relpath^$stylesheet" rel="stylesheet" type="text/css" />
$extrastylesheet
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<!--BEGIN TITLEAREA-->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<!--BEGIN PROJECT_LOGO-->
<td id="projectlogo"><img alt="Logo" src="$relpath^$projectlogo"/></td>
<!--END PROJECT_LOGO-->
<!--BEGIN PROJECT_NAME-->
<td style="padding-left: 0.5em;">
<div id="projectname">$projectname
<!--BEGIN PROJECT_NUMBER-->&#160;<span id="projectnumber">$projectnumber</span><!--END PROJECT_NUMBER-->
</div>
<!--BEGIN PROJECT_BRIEF--><div id="projectbrief">$projectbrief</div><!--END PROJECT_BRIEF-->
</td>
<!--END PROJECT_NAME-->
<!--BEGIN !PROJECT_NAME-->
<!--BEGIN PROJECT_BRIEF-->
<td style="padding-left: 0.5em;">
<div id="projectbrief">$projectbrief</div>
</td>
<!--END PROJECT_BRIEF-->
<!--END !PROJECT_NAME-->
<!--BEGIN DISABLE_INDEX-->
<!--BEGIN SEARCHENGINE-->
<td>$searchbox</td>
<!--END SEARCHENGINE-->
<!--END DISABLE_INDEX-->
</tr>
</tbody>
</table>
</div>
<!--END TITLEAREA-->
<!-- end header part -->
TODO.txt
View file @ 330cef96
mesh + !partitioner + !domain + solvers + +/-nonlinear solvers + autodiff + !system handler
:
0) () , data()
1) () erase array , resize'
......
examples/DrawMatrix/CMakeLists.txt
View file @ 330cef96
......@@ -30,7 +30,13 @@ if(OPENGL_FOUND)
message("linking DrawMatrix with Trilinos")
target_link_libraries(DrawMatrix ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_METIS)
message("linking DrawMatrix with Metis")
target_link_libraries(DrawMatrix ${METIS_LIBRARIES})
endif()
endif()
install(TARGETS DrawMatrix EXPORT inmost-targets RUNTIME DESTINATION bin)
else(GLUT_FOUND)
message("GLUT not found, not building DrawMatrix")
......
examples/DrawMatrix/main.cpp
View file @ 330cef96
......@@ -387,10 +387,10 @@ void keyboard(unsigned char key, int x, int y)
void DrawEntry(int i, int j)//, Storage::real r)
{
//~ glColor3f(r,0.0,1.0-r);
glVertex2i(i-(zoom-1),j-(zoom-1));
glVertex2i(i+zoom,j-(zoom-1));
glVertex2i(i+zoom,j+zoom);
glVertex2i(i-(zoom-1),j+zoom);
glVertex2i(i-(zoom-1),j-(zoom-1)-1);
glVertex2i(i+zoom,j-(zoom-1)-1);
glVertex2i(i+zoom,j+zoom-1);
glVertex2i(i-(zoom-1),j+zoom-1);
}
......
examples/FVDiscr/CMakeLists.txt
View file @ 330cef96
......@@ -21,6 +21,10 @@ if(USE_SOLVER)
message("linking FVDiscr with Trilinos")
target_link_libraries(FVDiscr ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_METIS)
message("linking FVDiscr with Metis")
target_link_libraries(FVDiscr ${METIS_LIBRARIES})
endif()
endif()
......
examples/MatSolve/CMakeLists.txt
View file @ 330cef96
......@@ -29,6 +29,23 @@ if(USE_SOLVER)
message("linking MatSolve with Trilinos")
target_link_libraries(MatSolve ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_METIS)
message("linking MatSolve with Metis")
target_link_libraries(MatSolve ${METIS_LIBRARIES})
endif()
endif()
if(USE_PARTITIONER)
if(USE_PARTITIONER_ZOLTAN)
message("linking MatSolve with Zoltan")
target_link_libraries(MatSolve ${ZOLTAN_LIBRARIES})
endif()
if(USE_PARTITIONER_PARMETIS)
message("linking MatSolve with ParMETIS")
target_link_libraries(MatSolve ${PARMETIS_LIBRARIES})
endif()
endif()
install(TARGETS MatSolve EXPORT inmost-targets RUNTIME DESTINATION bin)
examples/MatSolve/main.cpp
View file @ 330cef96
......@@ -22,13 +22,14 @@ int main(int argc, char ** argv)
switch(atoi(argv[1]))
{
case 0: type = Solver::INNER_ILU2; break;
case 1: type = Solver::INNER_MLILUC; break;
case 1: type = Solver::INNER_DDPQILUC; break;
case 2: type = Solver::PETSc; break;
case 3: type = Solver::Trilinos_Aztec; break;
case 4: type = Solver::Trilinos_Belos; break;
case 5: type = Solver::Trilinos_Ifpack; break;
case 6: type = Solver::Trilinos_ML; break;
case 7: type = Solver::ANI; break;
case 8: type = Solver::INNER_MPTILUC; break;
}
Solver::Initialize(&argc,&argv,argc > 4 ? argv[4] : NULL); // Initialize the linear solver in accordance with args
{
......@@ -71,15 +72,19 @@ int main(int argc, char ** argv)
{
Solver s(type); // Declare the linear solver by specified type
s.SetParameterEnum("gmres_substeps",3);
s.SetParameterEnum("gmres_substeps",4);
s.SetParameterReal("relative_tolerance",1.0e-9);
s.SetParameterReal("absolute_tolerance",1.0e-16);
s.SetParameterEnum("reorder_nonzeros",0);
s.SetParameterEnum("rescale_iterations",8);
s.SetParameterEnum("adapt_ddpq_tolerance",0);
s.SetParameterReal("drop_tolerance",0.001);
s.SetParameterReal("reuse_tolerance",0.00001);
s.SetParameterReal("ddpq_tolerance",0.7);
s.SetParameterReal("drop_tolerance",1.0e-4);
s.SetParameterReal("reuse_tolerance",1.0e-8);
s.SetParameterReal("ddpq_tolerance",0.8);
s.SetParameterEnum("condition_estimation",1);
t = Timer();
......
examples/Solver/CMakeLists.txt
View file @ 330cef96
......@@ -5,21 +5,28 @@ add_executable(Solver ${SOURCE})
target_link_libraries(Solver inmost)
if(USE_SOLVER_ANI)
message("linking Solver with ani3d and BLAS")
target_link_libraries(Solver ani3d ${BLAS_LIBRARIES})
if(BLAS_LINKER_FLAGS)
set_target_properties(Solver PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}")
if(USE_SOLVER)
if(USE_SOLVER_ANI)
message("linking Solver with ani3d and BLAS")
target_link_libraries(Solver ani3d ${BLAS_LIBRARIES})
if(BLAS_LINKER_FLAGS)
set_target_properties(Solver PROPERTIES LINK_FLAGS "${BLAS_LINKER_FLAGS}")
endif()
endif()
endif()
if(USE_SOLVER_PETSC)
message("linking Solver with PETSc")
target_link_libraries(Solver ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking Solver with Trilinos")
target_link_libraries(Solver ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_PETSC)
message("linking Solver with PETSc")
target_link_libraries(Solver ${PETSC_LIBRARIES})
endif()
if(USE_SOLVER_TRILINOS)
message("linking Solver with Trilinos")
target_link_libraries(Solver ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES})
endif()
if(USE_SOLVER_METIS)
message("linking Solver with Metis")
target_link_libraries(Solver ${METIS_LIBRARIES})
endif()
endif(USE_SOLVER)
if(USE_PARTITIONER)
if(USE_PARTITIONER_ZOLTAN)
message("linking Solver with Zoltan")
......@@ -29,7 +36,7 @@ if(USE_PARTITIONER)
message("linking Solver with ParMETIS")
target_link_libraries(Solver ${PARMETIS_LIBRARIES})
endif()
endif()
endif(USE_PARTITIONER)
if(USE_MPI)
message("linking Solver with MPI")
......
inmost_common.h
View file @ 330cef96
......@@ -15,6 +15,7 @@
//#define USE_PARTITIONER_PARMETIS
#define USE_SOLVER
//#define USE_SOLVER_METIS
//#define USE_SOLVER_PETSC
//#define USE_SOLVER_TRILINOS
//#define USE_SOLVER_ANI
......
inmost_mesh.h
View file @ 330cef96
......@@ -834,6 +834,27 @@ namespace INMOST
__INLINE const Edge & InvalidEdge() {static Edge ret(NULL,InvalidHandle()); return ret;}
/// \brief An interface for elements of type FACE.
///
/// This interface carry the link to the mesh and the element's handle that
/// represents position of the element's data in the mesh.
///
/// Interface provides some operations that can be done uniquely on the object of class
/// Element for which Element::GetElementType retruns FACE.
///
/// For the basic set of operations on all of the elements check class Element.
///
/// For the basic set of operations on the data of the element check class Storage.
///
/// You can obtain object of class Face from Mesh::iteratorFace,
/// in this case obtained object is always valid.
/// Also you can get it through Mesh::FaceByLocalID, check with Element::isValid to see whether you
/// have obtained a valid object.
/// You can convert an object of class Element into an object of class Face by Element::getAsFace. In debug
/// mode it will internally check that the element is of type FACE.
/// You may compose an object of class Face by using constructor and specifing pointer to the mesh
/// and providing element's handle. You can make a handle with ComposeHandle(FACE,local_id) function.
class Face : public Element //implemented in face.cpp
{
public:
......@@ -861,15 +882,17 @@ namespace INMOST
Node getEnd () const;
/// \brief Retrieve the cell for which the normal points outwards.
///
/// Depending on the grid construction the normal may incorrectly point inwards.
/// \warning Depending on the grid construction the normal may incorrectly point inwards.
/// You can resolve this situation by Face::FixNormalOrientation.
///
/// @return the cell for which normal points outwards.
/// @see Face::FixNormalOrientation
Cell BackCell () const;
/// \brief Retrieve the cell for which the normal points inwards.
///
/// Depending on the grid construction the normal may incorrectly point outwards.
/// \warning Depending on the grid construction the normal may incorrectly point outwards.
/// You can resolve this situation by Face::FixNormalOrientation.
///
/// @return the cell for which normal points inwards.
/// @see Face::FixNormalOrientation
Cell FrontCell () const;
......@@ -902,8 +925,8 @@ namespace INMOST
/// This interface carry the link to the mesh and the element's handle that
/// represents position of the element's data in the mesh.
///
/// Interface provides some operations that can be done uniquely on the element
/// of type CELL.
/// Interface provides some operations that can be done uniquely on the object of class
/// Element for which Element::GetElementType retruns CELL.
///
/// For the basic set of operations on all of the elements check class Element.
///
......@@ -913,10 +936,10 @@ namespace INMOST
/// in this case obtained object is always valid.
/// Also you can get it through Mesh::CellByLocalID, check with Element::isValid to see whether you
/// have obtained valid object.
/// You can convert object of class Element into object of class Cell by Element::getAsCell. In debug
/// mode it will internally check that element is of type CELL.
/// You may compose object of class Cell by using constructor and specifing pointer to the mesh
/// and providing element's handle. You can make handle with ComposeHandle(CELL,local_id) function.
/// You can convert an object of class Element into an object of class Cell by Element::getAsCell. In debug
/// mode it will internally check that the element is of type CELL.
/// You may compose an object of class Cell by using constructor and specifing pointer to the mesh
/// and providing element's handle. You can make a handle with ComposeHandle(CELL,local_id) function.
class Cell : public Element
{
public:
......@@ -1000,13 +1023,6 @@ namespace INMOST
///
/// The order of the edges is unspecified.
///
/// WARNING: Note that this function uses markers to check for duplication of edges
/// in output array. This allows for faster algorithm, but will lead to penalties
/// in shared parallel execution, since operations on markers are declared atomic.
/// For atomic declaration to work you have to define USE_OMP during CMake configuration
/// or in inmost_common.h. If you attempt to use this function in shared parallel execution
/// without USE_OMP you may expect side effects.
///
/// @return Set of edges that compose current cell.
///
/// \todo
......@@ -1017,6 +1033,13 @@ namespace INMOST
/// 1. Use of markers (current wariant).
/// 2. Put all elements into array with duplications, then run std::sort and std::unique.
/// 3. Put all elements into array, check for duplication by running through array.
///
/// \warning Note that this function uses markers to check for duplication of edges
/// in output array. This allows for faster algorithm, but will lead to penalties
/// in shared parallel execution, since operations on markers are declared atomic.
/// For atomic declaration to work you have to define USE_OMP during CMake configuration
/// or in inmost_common.h. If you attempt to use this function in shared parallel execution
/// without USE_OMP you may expect side effects.
ElementArray<Edge> getEdges () const;
/// \brief Get all the faces of the current cell.
///
......@@ -1037,13 +1060,13 @@ namespace INMOST
/// the order is unspecified. When suggest_nodes_order was provided into Mesh::CreateCell
/// then the order follows provided order.
///
/// WARNING: To work correctly in shared parallel environment this function requires
/// USE_OMP to be defined in CMake or in inmost_common.h.
///
/// @param mask Marker that should be used to filter elements.
/// @param invert_mask If false then those elements that are marked will be taken, otherwise
/// elements that are not marked will be taken.
/// @return Set of the nodes that compose current cell and (not) marked by marker.
///
/// \warning To work correctly in shared parallel environment this function requires
/// USE_OMP to be defined in CMake or in inmost_common.h.
ElementArray<Node> getNodes (MarkerType mask,bool invert_mask = false) const;
/// \brief Get the subset of the edges of the current cell that are (not) marked by provided marker.
///
......@@ -1051,13 +1074,13 @@ namespace INMOST
///
/// The order of the edges is unspecified.
///
/// WARNING: To work correctly in shared parallel environment this function requires
/// USE_OMP to be defined in CMake or in inmost_common.h.
///
/// @param mask Marker that should be used to filter elements.
/// @param invert_mask If false then those elements that are marked will be taken, otherwise
/// elements that are not marked will be taken.
/// @return Set of the edges that compose current cell and (not) marked by marker.
///
/// \warning To work correctly in shared parallel environment this function requires
/// USE_OMP to be defined in CMake or in inmost_common.h.
ElementArray<Edge> getEdges (MarkerType mask,bool invert_mask = false) const;
/// \brief Get the subset of the faces of the current cell that are (not) marked by provided marker.
///
......@@ -1066,13 +1089,13 @@ namespace INMOST
/// The order of the faces returned by this function is preserved from
/// the moment of the construction of the cell.
///
/// WARNING: To work correctly in shared parallel environment this function requires
/// USE_OMP to be defined in CMake or in inmost_common.h.
///
/// @param mask Marker that should be used to filter elements.
/// @param invert_mask If false then those elements that are marked will be taken, otherwise
/// elements that are not marked will be taken.
/// @return Set of the faces that compose current cell and (not) marked by marker.
///
/// \warning To work correctly in shared parallel environment this function requires
/// USE_OMP to be defined in CMake or in inmost_common.h.
ElementArray<Face> getFaces (MarkerType mask,bool invert_mask = false) const;
/// \brief Check that sequence of edges form a closed loop and each edge have a node that matches one of the nodes at the next edge.
///
......@@ -1185,6 +1208,7 @@ namespace INMOST
bool Closure () const; // test integrity of cell
};
///
__INLINE const Cell & InvalidCell() {static Cell ret(NULL,InvalidHandle()); return ret;}
class ElementSet : public Element //implemented in eset.cpp
......@@ -1530,6 +1554,9 @@ namespace INMOST
__INLINE const void * MGetLink (HandleType h, const Tag & t) const {if( !t.isSparseByDim(GetHandleElementNum(h)) ) return MGetDenseLink(h,t); else return MGetSparseLink(h,t);}
__INLINE void * MGetLink (HandleType h, const Tag & t) {if( !t.isSparseByDim(GetHandleElementNum(h)) ) return MGetDenseLink(h,t); else {void * & q = MGetSparseLink(h,t); if( q == NULL ) q = calloc(1,t.GetRecordSize()); return q;}}
public:
/// Remove all data and all elements from the mesh
/// Reset geometry service and topology check flags
void Clear();
/// For debug purposes
integer HandleDataPos (HandleType h) {return links[GetHandleElementNum(h)][GetHandleID(h)];}
/// For parmetis
......
inmost_options_cmake.h
View file @ 330cef96
......@@ -15,6 +15,7 @@
#cmakedefine USE_PARTITIONER_PARMETIS
#cmakedefine USE_SOLVER
#cmakedefine USE_SOLVER_METIS
#cmakedefine USE_SOLVER_PETSC
#cmakedefine USE_SOLVER_TRILINOS
#cmakedefine USE_SOLVER_ANI
......
inmost_solver.h
View file @ 330cef96
......@@ -42,7 +42,8 @@ namespace INMOST
enum Type
{
INNER_ILU2, ///< inner Solver based on BiCGStab(L) solver with second order ILU factorization as preconditioner.
INNER_MLILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and unsymmetric reordering for diagonal dominance as preconditioner.
INNER_DDPQILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and unsymmetric reordering for diagonal dominance as preconditioner.
INNER_MPTILUC, ///< inner Solver based on BiCGStab(L) solver with second order Crout-ILU with inversed-based condition estimation and maximum product transversal reordering as preconditioner.
Trilinos_Aztec, ///< external Solver AztecOO from Trilinos package
Trilinos_Belos, ///< external Solver Belos from Trilinos package, currently without preconditioner
Trilinos_ML, ///< external Solver AztecOO with ML preconditioner
......
mesh.cpp
View file @ 330cef96
......@@ -336,6 +336,49 @@ namespace INMOST
//AssignGlobalID(other.have_global_id);
return *this;
}
void Mesh::Clear()
{
for(ElementType etype = NODE; etype <= MESH; etype = NextElementType(etype))
{
for(tag_array_type::size_type i = 0; i < tags.size(); ++i)
{
if( tags[i].isDefined(etype) )
{
if( tags[i].isSparse(etype) )
{
for(integer lid = 0; lid < LastLocalID(etype); ++lid)
if( isValidElement(etype,lid) )
DelSparseData(ComposeHandle(etype,lid),tags[i]);
}
else if( tags[i].GetSize() == ENUMUNDEF )
{
for(integer lid = 0; lid < LastLocalID(etype); ++lid)
if( isValidElement(etype,lid) )
DelDenseData(ComposeHandle(etype,lid),tags[i]);
}
}
}
}
memset(remember,0,sizeof(bool)*15);
tags.clear();
//clear links
dense_data.clear();
for(int i = 0; i < 5; i++)
{
links[i].clear();
empty_links[i].clear();
empty_space[i].clear();
}
for(int i = 0; i < 6; i++)
{
sparse_data[i].clear();
back_links[i].clear();
}
RemTopologyCheck(ENUMUNDEF);
SetTopologyCheck(DEFAULT_CHECK);
}
Mesh::~Mesh()
{
......
solver.cpp
View file @ 330cef96
......@@ -4,6 +4,9 @@
#include "solver_ani.h"
#include "solver_ilu2.hpp"
#include "solver_ddpqiluc2.hpp"
#if defined(HAVE_SOLVER_MPTILUC2)
#include "solver_mtiluc2.hpp"
#endif
#include "solver_bcgsl.hpp"
#include <fstream>
#include <sstream>
......@@ -1354,7 +1357,7 @@ namespace INMOST
solver_data = static_cast<void *>(problem);
}
#endif
if( _pack == INNER_ILU2 || _pack == INNER_MLILUC )
if( _pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC)
{
Method * prec;
if (_pack == INNER_ILU2)
......@@ -1366,13 +1369,22 @@ namespace INMOST
prec = NULL;
#endif
}
else
else if( _pack == INNER_DDPQILUC )
{
#if defined(__SOLVER_DDPQILUC2__)
prec = new ILUC_preconditioner(info);
#else
std::cout << "Sorry, multilevel condition estimation crout-ilu preconditioner with reordering for diagonal dominance is not included in this release" << std::endl;
prec = NULL;
#endif
}
else if( _pack == INNER_MPTILUC )
{
#if defined(__SOLVER_MTILUC2__)
prec = new MTILUC_preconditioner(info);
#else
std::cout << "Sorry, maximum product transverse condition estimation crout-ilu preconditioner is not included in this release" << std::endl;
prec = NULL;
#endif
}
solver_data = new KSOLVER(prec, info);
......@@ -1420,7 +1432,7 @@ namespace INMOST
throw - 1; //You should not really want to copy solver's information
}
#endif
if (_pack == INNER_ILU2 || _pack == INNER_MLILUC)
if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC)
{
throw - 1; //You should not really want to copy solver's information
}
......@@ -1471,7 +1483,7 @@ namespace INMOST
}
}
#endif
if (_pack == INNER_ILU2 || _pack == INNER_MLILUC)
if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC)
{
if( matrix_data != NULL )
{
......@@ -1559,7 +1571,7 @@ namespace INMOST
throw - 1; //You should not really want to copy solver's information
}
#endif
if (_pack == INNER_ILU2 || _pack == INNER_MLILUC)
if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC)
{
throw - 1; //You should not really want to copy solver's information
}
......@@ -1815,7 +1827,7 @@ namespace INMOST
ok = true;
}
#endif
if (_pack == INNER_ILU2 || _pack == INNER_MLILUC)
if (_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC)
{
Solver::Matrix * mat = new Solver::Matrix(A);
......@@ -1829,13 +1841,17 @@ namespace INMOST
sol->RealParameter(":tau2") = preconditioner_reuse_tolerance;
sol->EnumParameter(":scale_iters") = preconditioner_rescale_iterations;
if( _pack == INNER_MLILUC )
if( _pack == INNER_DDPQILUC )
{
sol->RealParameter(":ddpq_tau") = preconditioner_ddpq_tolerance;
sol->EnumParameter(":reorder_nnz") = preconditioner_reorder_nonzero;
sol->EnumParameter(":estimator") = preconditioner_condition_estimation;
sol->EnumParameter(":ddpq_tau_adapt") = preconditioner_adapt_ddpq_tolerance;
}
else if( _pack == INNER_MPTILUC )
{
sol->EnumParameter(":estimator") = preconditioner_condition_estimation;
}
else sol->EnumParameter(":fill") = static_cast<INMOST_DATA_ENUM_TYPE>(preconditioner_fill_level);
if( sizeof(KSOLVER) == sizeof(BCGSL_solver) )
......@@ -2224,7 +2240,7 @@ namespace INMOST
}
#endif
if(_pack == INNER_ILU2 || _pack == INNER_MLILUC)
if(_pack == INNER_ILU2 || _pack == INNER_DDPQILUC || _pack == INNER_MPTILUC)
{
IterativeMethod * sol = static_cast<IterativeMethod *>(solver_data);
......
solver_bcgsl.hpp
View file @ 330cef96
......@@ -14,6 +14,92 @@
namespace INMOST
{
int solvenxn(INMOST_DATA_REAL_TYPE * A, INMOST_DATA_REAL_TYPE * x, INMOST_DATA_REAL_TYPE * b, int n, int * order)
{
INMOST_DATA_REAL_TYPE temp, max;
int temp2;
for(int i = 0; i < n; i++) order[i] = i;
for(int i = 0; i < n; i++)
{
int maxk = i, maxq = i;
max = fabs(A[maxk*n+maxq]);
//Find best pivot
for(int q = i; q < n; q++) // over columns
{
for(int k = i; k < n; k++) // over rows
{
if( fabs(A[k*n+q]) > max )
{
max = fabs(A[k*n+q]);
maxk = k;
maxq = q;
}
}
}
//Exchange rows
if( maxk != i )
{
for(int q = 0; q < n; q++)
{
temp = A[maxk*n+q];
A[maxk*n+q] = A[i*n+q];