Skip to content

GitLab

Commit b53160fb authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

remove support of old fortran ilu2 solver

parent a03370ea
Expand all Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ b53160fb
......@@ -26,7 +26,6 @@ option(USE_SOLVER "Compile solver capabilities" ON)
option(USE_PARTITIONER "Compile partitioner capabilities" ON)
option(USE_AUTODIFF "Compile automatic differentiation capabilities" ON)
option(USE_NONLINEAR "Compile nonlinear solver capabilities" ON)
option(TEST_FORTRAN_ANI3D "Test for fortran availibility to compile ANI3D lib" OFF)
option(COMPILE_EXAMPLES "Compile examples" OFF)
option(COMPILE_TESTS "Compile some tests" OFF)
......@@ -69,21 +68,6 @@ foreach(CompilerFlag ${CompilerFlags})
endforeach()
endif( USE_MT )
if(TEST_FORTRAN_ANI3D)
include(CheckLanguage)
check_language(Fortran)
find_package(BLAS)
if(NOT CMAKE_Fortran_COMPILER OR NOT BLAS_FOUND)
message(STATUS "No Fortran or BLAS support")
option(USE_SOLVER_ANI "Compile with ANI3D sequential solvers" OFF)
else()
enable_language(Fortran)
add_subdirectory(ILU2)
option(USE_SOLVER_ANI "Compile with ANI3D sequential solvers" ON)
endif()
endif()
if(USE_OPTIMIZER_BAYESIAN)
include_directories(/Users/bvdmitri/Projects/CXX/limbo/src)
target_link_libraries(inmost /Users/bvdmitri/Projects/CXX/limbo/build/src/liblimbo.a)
......
ILU2/CMakeLists.txt deleted 100644 → 0
View file @ a03370ea
cmake_minimum_required(VERSION 2.8.7)
project(ani3d)
set(SOURCE bcg.f
matvecCSR.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()
ILU2/Makefile deleted 100644 → 0
View file @ a03370ea
############################################################
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)
############################################################
ILU2/README deleted 100644 → 0
View file @ a03370ea
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
ILU2/README.md deleted 100644 → 0
View file @ a03370ea
ANI3D solvers
=======
This folder contains solvers from Ani3D package from:
http://sourceforge.net/projects/ani3d/
ILU2/Rules.make deleted 100644 → 0
View file @ a03370ea
############################################################
# 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 =
ILU2/bcg.f deleted 100644 → 0
View file @ a03370ea
c---->------------------------------------------------------------------<
c Preconditioned bi-conjugate gradient stabilized method
c Templates for the solution if linear Systems...
c http://www.netlib.org
c---->------------------------------------------------------------------<
SUBROUTINE slpbcgs(
> prevec, IPREVEC, iW,rW,
> matvec, IMATVEC, IA,JA,A,
> WORK, MW, NW,
> N, RHS, SOL,
> ITER, RESID,
> INFO, NUNIT )
c---->
IMPLICIT NONE
c---->------------------------------------------------------------------<
c Argument types:
c
EXTERNAL matvec, prevec
INTEGER IMATVEC(*), IPREVEC(*), IA(*), JA(*), iW(*)
INTEGER N, MW, NW, ITER, INFO, NUNIT
REAL*8 RESID, A(*)
REAL*8 RHS(*), SOL(*), WORK(MW,NW), rW(*)
c---->
c Argument Descriptions:
c
c prevec : extern : Precondition-vector routine
c IPREVEC : input : Configuration data for 'prevec'
c matvec : extern : Matrix-vector multiply routine
c IMATVEC : input : Configuration data for 'matvec'
c
c WORK : work : Workspace (MW,NW)
c MW : input : leading dimension of workspace >= N
c NW : input : trailing dimension of workspace >= 8
c
c N : input : Length of vectors
c RHS : input : RHS vector
c SOL : in/out : Initial guess / iterated solution
c ITER : in/out : Maximum iterations / actual iterations
c RESID : in/out : Convergence target / Norm of final residual
c INFO : output : = 0, converged
c : = 1, did not converge
c : = 2, could not continue since OMEGA = ZERO
c : = 3, could not continue since || S || is too small
c : = 4, could not continue since RHO = ZERO
c : < 0, error with input
c---->
c External routine specifications:
c
c matvec( IMATVEC, A, X, B, Y ) <=> Y = A * Mat * X + B * Y
c prevec( IPREVEC, i, X, Y ) <=> Y = (MatP_{i})^{-1} * X
c where MatP is the approximation of Mat
c---->------------------------------------------------------------------<
c Local Parameters
c
REAL*8 ZERO,ONE
PARAMETER ( ZERO = 0.0 , ONE = 1.0 )
c---->------------------------------------------------------------------<
c Local Variables:
c
INTEGER MAXIT
INTEGER JR, JP, JS, JR1, JP1, JS1, JV, JT, i
REAL*8 RHO, RHOPREV, OMEGA, ALPHA, BETA, TOL, TMP(1), TMP2
& ,minSol,maxSol
c
c---->------------------------------------------------------------------<
c External BLAS, etc.:
c
EXTERNAL ddot,daxpy,dcopy,dscal
REAL*8 ddot
INTRINSIC sqrt, min
c---->------------------------------------------------------------------<
c
c Test the input parameters.
c
INFO = 0
c
if ( N .eq. 0 ) then
return
else if ( N .lt. 0 ) then
INFO = -10
else if ( MW .lt. N ) then
INFO = -20
else if ( NW .lt. 8 ) then
INFO = -30
else if ( ITER .le. 0 ) then
INFO = -40
endif
c
if ( INFO .ne. 0 ) return
c---->------------------------------------------------------------------<
c Save input iteration limit and convergence tolerance
c
MAXIT = ITER
TOL = RESID
c---->
c Alias workspace columns.
c
JV = 1
JT = JV + 1
JR = JT + 1
JP = JR + 1
JS = JP + 1
JR1 = JS + 1
JP1 = JR1+ 1
JS1 = JP1+ 1
c---->
c Set initial residual
c
call dcopy( N, RHS, 1, WORK(1,JR), 1 )
c
TMP2 = ddot( N, SOL, 1, SOL, 1 )
if ( TMP2 .ne. ZERO ) then
call matvec( IMATVEC, -ONE, SOL, ONE, WORK(1,JR), IA,JA,A )
c call matvec( IMATVEC, -ONE, SOL, ZERO, WORK(1,JR) )
c call daxpy( N, ONE, RHS, 1, WORK(1,JR), 1 )
endif
c---->
TMP2 = ddot( N, WORK(1,JR), 1, WORK(1,JR), 1 )
RESID = sqrt( TMP2 )
c---->
ITER = 0
RHO = ZERO
if ( RESID .lt. TOL ) GOTO 20
call dcopy( N, WORK(1,JR), 1, WORK(1,JR1), 1 )
c---->------------------------------------------------------------------<
c PBCGS iteration point
c---->--
10 continue
c
ITER = ITER + 1
c---->----
RHOPREV = RHO
RHO = ddot( N, WORK(1,JR), 1, WORK(1,JR1), 1 )
if ( RHO .eq. ZERO ) then
if ( NUNIT .gt. 0 ) then
write(NUNIT,*) "PBCGS: Bad rho_tilde: method fails"
end if
INFO = 4
goto 20
end if
IF (ITER.eq.1) THEN
call dcopy( N, WORK(1,JR), 1, WORK(1,JP), 1 )
TMP(1) = ZERO
call dcopy( N, TMP, 0, WORK(1,JV), 1 )
call dcopy( N, TMP, 0, WORK(1,JT), 1 )
ELSE
BETA = ( RHO / RHOPREV ) * ( ALPHA / OMEGA )
call daxpy( N, -OMEGA, WORK(1,JV), 1, WORK(1,JP), 1 )
call dscal( N, BETA, WORK(1,JP), 1 )
call daxpy( N, ONE, WORK(1,JR), 1, WORK(1,JP), 1 )
END IF
call prevec( IPREVEC, 0, WORK(1,JP), WORK(1,JP1), iW,rW )
call matvec( IMATVEC, ONE, WORK(1,JP1), ZERO, WORK(1,JV),
& IA,JA,A )
TMP2 = ddot( N, WORK(1,JV), 1, WORK(1,JR1), 1 )
ALPHA = RHO / TMP2
call dcopy( N, WORK(1,JR), 1, WORK(1,JS), 1 )
call daxpy( N, -ALPHA, WORK(1,JV), 1, WORK(1,JS), 1 )
TMP2 = ddot( N, WORK(1,JS), 1, WORK(1,JS), 1 )
if ( sqrt( TMP2 ) .lt. 1d-16 ) then
call daxpy( N, ALPHA, WORK(1,JP1), 1, SOL, 1 )
call matvec( IMATVEC, -ONE, SOL, ZERO, WORK(1,JR), IA,JA,A )
call daxpy( N, ONE, RHS, 1, WORK(1,JR), 1 )
TMP2 = ddot( N, WORK(1,JR), 1, WORK(1,JR), 1 )
RESID = sqrt( TMP2 )
INFO = 3
goto 20
end if
call prevec( IPREVEC, 0, WORK(1,JS), WORK(1,JS1), iW,rW )
call matvec( IMATVEC, ONE, WORK(1,JS1), ZERO, WORK(1,JT),
& IA,JA,A )
OMEGA = ddot( N, WORK(1,JT), 1, WORK(1,JS), 1 )
TMP2 = ddot( N, WORK(1,JT), 1, WORK(1,JT), 1 )
OMEGA = OMEGA / TMP2
call daxpy( N, ALPHA, WORK(1,JP1), 1, SOL, 1 )
call daxpy( N, OMEGA, WORK(1,JS1), 1, SOL, 1 )
call dcopy( N, WORK(1,JS), 1, WORK(1,JR), 1 )
call daxpy( N, -OMEGA, WORK(1,JT), 1, WORK(1,JR), 1 )
c---->----
c Check convergence
TMP2 = ddot( N, WORK(1,JR), 1, WORK(1,JR), 1 )
RESID = sqrt( TMP2 )
c---->------------------------------------------------------------------<
c Continue BPCGS loop while:
c 1) Less than maximum iteration
c 2) Have not converged
c print*,'pbcgs: ',ITER, RESID
if (NUNIT.gt.0) write(NUNIT,*) ITER, RESID , ''
c For continuation it is necessary that OMEGA .ne. ZERO
c
if ( ITER .lt. MAXIT .and. RESID .ge. TOL
> .and. OMEGA .ne. ZERO ) go to 10
c---->--
c
c Convergence failure?
c
if ( ITER .ge. MAXIT .and. RESID .ge. TOL ) INFO = 1
if ( OMEGA .eq. ZERO ) INFO = 2
c---->------------------------------------------------------------------<
c Output
c
20 continue
TMP2 = ddot( N , SOL, 1, SOL, 1 )
TMP2 = sqrt( TMP2 )
if ( NUNIT .gt. 0 ) then
minSol=1d30
maxSol=-1d30
do i=1,n
minSol=min(SOL(i),minSol)
maxSol=max(SOL(i),maxSol)
end do
WRITE(NUNIT,9000) ITER,RESID,TMP2,minSol,maxSol
9000 FORMAT('ILU: iters:',I4,', residual=',E9.3,' (|SOL|= ',E10.4,
& ', ',E11.4,' < SOL < ',E11.4,')')
end if
c---->------------------------------------------------------------------<
return
end
c---->------------------------------------------------------------------<
ILU2/dsort.f deleted 100644 → 0
View file @ a03370ea
!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)