Commit 1dd9fb98 authored by Kirill Terekhov's avatar Kirill Terekhov

initial commit

parents
projects/
AsmJit/
compile32_vs2013/
compile64_vs2013/
development/
inmost_autodiff.h
autodiff.cpp
TODO.txt
examples/
examples/Data/
examples/Matrices/
examples/Grids/
tests/
\ No newline at end of file
cmake_minimum_required (VERSION 2.6)
project (INMOST)
set(SOURCE solver.cpp
solver_ani.cpp
solver_petsc.cpp
partitioner.cpp
algorithm.cpp
geometry.cpp
iterator.cpp
storage.cpp
eset.cpp
mesh_file.cpp
timer.cpp
face.cpp
edge.cpp
mesh.cpp
node.cpp
cell.cpp
tag.cpp
element.cpp
mesh_parallel.cpp
modify.cpp
autodiff.cpp)
set(HEADER inmost.h
inmost_options_cmake.h
inmost_common.h
inmost_mesh.h
inmost_solver.h
inmost_partitioner.h
inmost_autodiff.h
container.hpp
io.hpp
solver_ilu2.hpp
# solver_ddpqiluc2.hpp
solver_bcgsl.hpp
solver_prototypes.hpp)
add_library(inmost STATIC ${SOURCE} ${HEADER})
option(USE_MPI "Compile with MPI support" ON)
option(USE_OMP "Compile with OpenMP support (experimental)" OFF)
option(USE_MESH "Compile mesh capabilities" ON)
option(USE_SOLVER "Compile solver capabilities" ON)
option(USE_PARTITIONER "Compile partitioner capabilities" ON)
option(USE_AUTODIFF "Compile automatic differentiation capabilities" ON)
option(TEST_FORTRAN_ANI3D "Test for fortran availibility to compile ANI3D lib" OFF)
option(COMPILE_EXAMPLES "Compile examples" OFF)
option(COMPILE_PROJECTS "Compile projects" OFF)
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_PETSC "Use PETSc solver" OFF)
option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation (under work)" OFF)
option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" ON)
option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake_find/")
if(TEST_FORTRAN_ANI3D)
add_subdirectory(ILU2)
if(ANI3D_AVAILIBLE)
option(USE_SOLVER_ANI "Compile with ANI3D sequential solvers" OFF)
else()
option(USE_SOLVER_ANI "Compile with ANI3D sequential solvers" ON)
endif()
endif()
if(USE_MPI)
find_package(MPI)
if(NOT MPI_FOUND)
option(USE_MPI "Use Message Passing Interface" OFF)
option(USE_MPI2 "Use MPI-2 functions" OFF)
message("MPI NOT FOUND")
else()
include_directories(${MPI_INCLUDE_PATH})
option(USE_MPI "Use Message Passing Interface" ON)
option(USE_MPI2 "Use MPI-2 functions" ON)
set_target_properties(inmost PROPERTIES COMPILE_FLAGS "${MPI_COMPILE_FLAGS}")
message("MPI FOUND")
endif()
endif()
if(USE_OMP)
find_package(OpenMP)
if (OPENMP_FOUND)
set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
else()
option(USE_OMP "Compile with OpenMP support (experimental)" OFF)
endif()
endif()
if(USE_PARTITIONER_PARMETIS)
find_package(ParMETIS)
if(NOT PARMETIS_FOUND)
option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" OFF)
message("PARMETIS NOT FOUND")
else()
include_directories(${PARMETIS_INCLUDE_DIR})
option(USE_PARTITIONER_PARMETIS "Use ParMetis partitioner" ON)
message("PARMETIS FOUND")
endif()
endif()
if(USE_PARTITIONER_ZOLTAN)
find_package(ZOLTAN)
if(NOT ZOLTAN_FOUND)
option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" OFF)
message("ZOLTAN NOT FOUND")
else()
include_directories(${ZOLTAN_INCLUDE_DIRS})
option(USE_PARTITIONER_ZOLTAN "Use Zoltan partitioner" ON)
message("ZOLTAN FOUND")
endif()
endif()
if(USE_SOLVER_PETSC)
find_package(PETSc)
if(NOT PETSC_FOUND)
option(USE_SOLVER_PETSC "Use PETSc solver" OFF)
message("PETSC NOT FOUND")
else()
include_directories(${PETSC_INCLUDES})
option(USE_SOLVER_PETSC "Use PETSc solver" ON)
message("PETSC FOUND")
add_definitions(${PETSC_DEFINITIONS})
#message(${PETSC_LIBRARIES})
endif()
endif()
if(USE_AUTODIFF_OPENCL)
find_package(OpenCL)
if(OPENCL_FOUND)
option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation" ON)
include_directories(${OPENCL_INCLUDE_DIRS})
else()
option(USE_AUTODIFF_OPENCL "Use OpenCL for automatic differentiation" OFF)
message("OpenCL not found")
endif()
endif()
if(USE_AUTODIFF_ASMJIT)
Set(ASMJIT_STATIC TRUE)
include("AsmJit/CMakeLists.txt")
include_directories("${ASMJIT_INC_DIR}/asmjit")
target_link_libraries(inmost asmjit)
endif()
if(USE_AUTODIFF_EXPRESSION_TEMPLATES)
if(MSVC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /bigobj")
endif()
endif()
configure_file("inmost_options_cmake.h" "${PROJECT_BINARY_DIR}/inmost_options.h")
include_directories("${PROJECT_BINARY_DIR}")
if(COMPILE_EXAMPLES)
add_subdirectory(examples)
endif(COMPILE_EXAMPLES)
if(COMPILE_PROJECTS)
add_subdirectory(projects)
endif(COMPILE_PROJECTS)
if(COMPILE_TESTS)
add_subdirectory(tests)
endif(COMPILE_TESTS)
cmake_minimum_required(VERSION 2.8.7)
project(ani3d)
set(SOURCE bcg.f
dsort.f
ilu0.f
iluoo.f)
#if(NOT WIN32 AND NOT CMAKE_GENERATOR STREQUAL Xcode AND NOT CMAKE_GENERATOR STREQUAL Ninja)
#--Work-around for CMake issue 0009220
if(DEFINED CMAKE_Fortran_COMPILER AND CMAKE_Fortran_COMPILER MATCHES "^$")
set(CMAKE_Fortran_COMPILER CMAKE_Fortran_COMPILER-NOTFOUND)
endif()
enable_language(Fortran OPTIONAL)
#endif()
#tutorial for better fortran support:
#http://www.netlib.org/lapack/lawnspdf/lawn270.pdf
if(CMAKE_Fortran_COMPILER_WORKS)
message("Fortran found, adding ani3d")
add_library(ani3d STATIC ${SOURCE})
set(ANI3D_AVAILIBLE 1)
else()
message("Fortran not found")
set(ANI3D_AVAILIBLE 0)
endif()
############################################################
SYSHOME := $(shell /bin/pwd)
include $(SYSHOME)/Rules.make
############################################################
LIBOBJS = bcg.o iluoo.o ilu0.o dsort.o
############################################################
all: help
lib: $(LIBILU) info
clean:
@/bin/rm -f *.o *~
help:
@echo "make {lib|clean|help}"
@echo " "
@echo " lib - create library libilu.a"
@echo " clean - clean this directory"
@echo " help - print help message"
@echo " "
info:
@echo " "
@echo "Libraries are located in lib/"
@echo " "
############################################################
.f.o:
@echo '$(F77) -c $(FFLAGS) ' $*.f
@$(F77) $(FFLAGS) -c $*.f -o $*.o
############################################################
$(LIBILU): $(LIBOBJS)
@echo ''
@rm -f $(LIBILU)
@ar -r $(LIBILU) $(LIBOBJS)
############################################################
Package name: Ani3D-ILU
(Incomplete LU solver)
1. What is aniILU?
The FORTRAN package aniILU solves systems of linear equations
with sparse non-singular matrices.
2. Quick Start
1. make lib
############################################################
# Fortran and C compilers for SEQuential and PARallel codes.
############################################################
F77 = gfortran -m64 #compilers
CC = gcc
CXX = g++
VIEWER = gmv #visualization of GMV files
FLINKER = $(F77) #linkers
CLINKER = $(CC)
LD = $(CXX)
LIBSYS = -lgfortran #for linking C and Fortran
############################################################
FFLAGS = -O5 #-fopenmp #-Wall
CFLAGS = -O5
LDFLAGS =
############################################################
ANIHOME = $(SYSHOME)
ANIILU = $(ANIHOME)/src/aniILU
ANILIB = $(ANIHOME)/lib
ANIINC = $(ANIHOME)/include
ANIBLAS = $(ANIHOME)/blas
ANILAPACK = $(ANIHOME)/src/lapack
############################################################
LIBILU = $(ANILIB)/libilu-2.3.a
#LIBBLAS = -lblas
#LIBLAPACK = -llapack
LIBBLAS = $(ANILIB)/libblas-3.0.a
LIBLAPACK = $(ANILIB)/liblapack-3.0.a
############################################################
INCLUDE =
!DEC$ ATTRIBUTES DLLEXPORT :: DSORTILU
SUBROUTINE DSORTilu (DX, IX, N, KFLAG, IERR)
C***BEGIN PROLOGUE DSORT
C***PURPOSE Sort an array and optionally make the same interchanges in
C an auxiliary array. The array may be sorted in increasing
C or decreasing order. A slightly modified QUICKSORT
C algorithm is used.
C***LIBRARY SLATEC
C***CATEGORY N6A2B
C***TYPE DOUBLE PRECISION (SSORT-S, DSORT-D, ISORT-I)
C***KEYWORDS SINGLETON QUICKSORT, SORT, SORTING
C***AUTHOR Jones, R. E., (SNLA)
C Wisniewski, J. A., (SNLA)
C***MODIFIED by K.Lipnikov (http://humpty.inm.ras.ru/~kos)
C***DESCRIPTION
C
C DSORT sorts array DX and optionally makes the same interchanges in
C array IX. The array DX may be sorted in increasing order or
C decreasing order. A slightly modified quicksort algorithm is used.
C
C Description of Parameters
C DX - array of values to be sorted (usually abscissas)
C IX - array to be (optionally) carried along
C N - number of values in array DX to be sorted
C KFLAG - control parameter
C = 2 means sort DX in increasing order and carry IX along.
C = 1 means sort DX in increasing order (ignoring IX)
C = -1 means sort DX in decreasing order (ignoring IX)
C = -2 means sort DX in decreasing order and carry IX along.
C
C IERR - the error's indicator
C = 0 - Okay
C != 0 - errors
C
C***REFERENCES R. C. Singleton, Algorithm 347, An efficient algorithm
C for sorting with minimal storage, Communications of
C the ACM, 12, 3 (1969), pp. 185-187.
C***ROUTINES CALLED XERMSG
C***REVISION HISTORY (YYMMDD)
C 761101 DATE WRITTEN
C 761118 Modified to use the Singleton quicksort algorithm. (JAW)
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891009 Removed unreferenced statement labels. (WRB)
C 891024 Changed category. (WRB)
C 891024 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
C 901012 Declared all variables; changed X,Y to DX,IX; changed
C code to parallel SSORT. (M. McClain)
C 920501 Reformatted the REFERENCES section. (DWL, WRB)
C 920519 Clarified error messages. (DWL)
C 920801 Declarations section rebuilt and code restructured to use
C IF-THEN-ELSE-ENDIF. (RWC, WRB)
C***END PROLOGUE DSORT
C .. Scalar Arguments ..
INTEGER KFLAG, N, IERR
C .. Array Arguments ..
INTEGER DX(*)
DOUBLE PRECISION IX(*)
C .. Local Scalars ..
DOUBLE PRECISION R
INTEGER T, TT
DOUBLE PRECISION ITX,ITTX
INTEGER I, IJ, J, K, KK, L, M, NN
C .. Local Arrays ..
INTEGER IL(21), IU(21)
C .. Intrinsic Functions ..
INTRINSIC ABS, INT
C***FIRST EXECUTABLE STATEMENT DSORT
IERR = 0
NN = N
IF (NN .LT. 1) THEN
IERR = 1201
RETURN
ENDIF
C
KK = ABS(KFLAG)
IF (KK.NE.1 .AND. KK.NE.2) THEN
IERR = 1202
RETURN
ENDIF
C
C
C Alter array DX to get decreasing order if needed
C
IF (KFLAG .LE. -1) THEN
DO 10 I=1,NN
DX(I) = -DX(I)
10 CONTINUE
ENDIF
C
IF (KK .EQ. 2) GO TO 100
C
C Sort DX only
C
M = 1
I = 1
J = NN
R = 0.375D0
C
20 IF (I .EQ. J) GO TO 60
IF (R .LE. 0.5898437D0) THEN
R = R+3.90625D-2
ELSE
R = R-0.21875D0
ENDIF
C
30 K = I
C
C Select a central element of the array and save it in location T
C
IJ = I + INT((J-I)*R)
T = DX(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
ENDIF
L = J
C
C If last element of array is less than than T, interchange with T
C
IF (DX(J) .LT. T) THEN
DX(IJ) = DX(J)
DX(J) = T
T = DX(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
ENDIF
ENDIF
C
C Find an element in the second half of the array which is smaller
C than T
C
40 L = L-1
IF (DX(L) .GT. T) GO TO 40
C
C Find an element in the first half of the array which is greater
C than T
C
50 K = K+1
IF (DX(K) .LT. T) GO TO 50
C
C Interchange these elements
C
IF (K .LE. L) THEN
TT = DX(L)
DX(L) = DX(K)
DX(K) = TT
GO TO 40
ENDIF
C
C Save upper and lower subscripts of the array yet to be sorted
C
IF (L-I .GT. J-K) THEN
IL(M) = I
IU(M) = L
I = K
M = M+1
ELSE
IL(M) = K
IU(M) = J
J = L
M = M+1
ENDIF
GO TO 70
C
C Begin again on another portion of the unsorted array
C
60 M = M-1
IF (M .EQ. 0) GO TO 190
I = IL(M)
J = IU(M)
C
70 IF (J-I .GE. 1) GO TO 30
IF (I .EQ. 1) GO TO 20
I = I-1
C
80 I = I+1
IF (I .EQ. J) GO TO 60
T = DX(I+1)
IF (DX(I) .LE. T) GO TO 80
K = I
C
90 DX(K+1) = DX(K)
K = K-1
IF (T .LT. DX(K)) GO TO 90
DX(K+1) = T
GO TO 80
C
C Sort DX and carry IX along
C
100 M = 1
I = 1
J = NN
R = 0.375D0
C
110 IF (I .EQ. J) GO TO 150
IF (R .LE. 0.5898437D0) THEN
R = R+3.90625D-2
ELSE
R = R-0.21875D0
ENDIF
C
120 K = I
C
C Select a central element of the array and save it in location T
C
IJ = I + INT((J-I)*R)
T = DX(IJ)
ITX = IX(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
IX(IJ) = IX(I)
IX(I) = ITX
ITX = IX(IJ)
ENDIF
L = J
C
C If last element of array is less than T, interchange with T
C
IF (DX(J) .LT. T) THEN
DX(IJ) = DX(J)
DX(J) = T
T = DX(IJ)
IX(IJ) = IX(J)
IX(J) = ITX
ITX = IX(IJ)
C
C If first element of array is greater than T, interchange with T
C
IF (DX(I) .GT. T) THEN
DX(IJ) = DX(I)
DX(I) = T
T = DX(IJ)
IX(IJ) = IX(I)
IX(I) = ITX
ITX = IX(IJ)
ENDIF
ENDIF
C
C Find an element in the second half of the array which is smaller
C than T
C
130 L = L-1
IF (DX(L) .GT. T) GO TO 130
C
C Find an element in the first half of the array which is greater
C than T
C
140 K = K+1
IF (DX(K) .LT. T) GO TO 140
C
C Interchange these elements
C
IF (K .LE. L) THEN
TT = DX(L)
DX(L) = DX(K)
DX(K) = TT
ITTX = IX(L)
IX(L) = IX(K)
IX(K) = ITTX
GO TO 130
ENDIF
C
C Save upper and lower subscripts of the array yet to be sorted
C
IF (L-I .GT. J-K) THEN
IL(M) = I
IU(M) = L
I = K
M = M+1
ELSE
IL(M) = K
IU(M) = J
J = L
M = M+1
ENDIF
GO TO 160
C
C Begin again on another portion of the unsorted array
C
150 M = M-1
IF (M .EQ. 0) GO TO 190
I = IL(M)
J = IU(M)
C
160 IF (J-I .GE. 1) GO TO 120
IF (I .EQ. 1) GO TO 110
I = I-1
C
170 I = I+1
IF (I .EQ. J) GO TO 150
T = DX(I+1)
ITX = IX(I+1)
IF (DX(I) .LE. T) GO TO 170
K = I
C
180 DX(K+1) = DX(K)
IX(K+1) = IX(K)
K = K-1
IF (T .LT. DX(K)) GO TO 180
DX(K+1) = T
IX(K+1) = ITX
GO TO 170
C
C Clean up
C
190 IF (KFLAG .LE. -1) THEN
DO 200 I=1,NN
DX(I) = -DX(I)
200 CONTINUE
ENDIF
RETURN
END
c-----------------------------------------------------------------------
c This routine solves the system (LU) y = x
c-----------------------------------------------------------------------
SUBROUTINE prevec0( iprevec, dummy, x, y, iwork, dwork)
INTEGER iprevec(*), dummy, iwork(*), n
REAL*8 x(*),y(*),dwork(*)
n = iprevec(1)
call lusol(n, x, y, dwork(1), iwork(n+2), iwork(1))
return
end
subroutine ilu0(n, a, ja, ia, alu, jlu, ju, iw, ierr)
implicit real*8 (a-h,o-z)
real*8 a(*), alu(*)
integer ja(*), ia(*), ju(*), jlu(*), iw(*)
c input:
c---------
c n = dimension of matrix
c a, ja,
c ia = original matrix in compressed sparse row storage.
c
c output:
c alu,jlu = matrix stored in Modified Sparse Row (MSR) format containing
c the L and U factors together. The diagonal (stored in
c alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix
c contains the i-th row of L (excluding the diagonal entry=1)
c followed by the i-th row of U.
c
c ju = pointer to the diagonal elements in alu, jlu.
c
c ierr = integer indicating error code on return
c ierr = 0 --> normal return
c ierr = k --> code encountered a zero pivot at step k.
c work arrays:
c-------------
c iw = integer work array of length n.
c-----------------------------------------------------------------------
integer maxbuf, nbuf
parameter (maxbuf=1000) ! max number of non-zeros in a row
integer jbuf(maxbuf)
real*8 abuf(maxbuf)
c-----------------------------------------------------------------------
ju0 = n+2
jlu(1) = ju0
c
c initialize work vector to zero's
c
do 31 i=1, n
iw(i) = 0
31 continue
c
c main loop
c
do 500 ii = 1, n
js = ju0
c
c sort column indexes
c
nbuf=0
do j=ia(ii),ia(ii+1)-1
nbuf=nbuf+1
abuf(nbuf)=a(j)
jbuf(nbuf)=ja(j)
end do
if (nbuf.gt.maxbuf) stop 'ilu0: maxbuf too small'
call dsortilu (jbuf, abuf, nbuf, 2, ierr)
c generating row number ii of L and U.
c
do 100 j=1,nbuf ! j=ia(ii),ia(ii+1)-1
c
c copy row ii of a, ja, ia into row ii of alu, jlu (L/U) matrix.
c
jcol = jbuf(j) ! ja(j)
if (jcol .eq. ii) then
alu(ii) = abuf(j) ! a(j)
iw(jcol) = ii
ju(ii) = ju0
else