Skip to content
Permalink
Browse files

ILU2: added missing files

Added missing files from ani3D-3.0/aniILU.
Update SOLVER_ANI code.
  • Loading branch information...
a-a-danilov committed Jan 31, 2015
1 parent 6938a77 commit d13bcf941f8457f746ae28225b4329d5cdb491ce
Showing with 274 additions and 5 deletions.
  1. +1 −0 ILU2/CMakeLists.txt
  2. +237 −0 ILU2/bcg.f
  3. +31 −0 ILU2/matvecCSR.f
  4. +1 −1 solver.cpp
  5. +4 −4 solver_ani.h
@@ -1,6 +1,7 @@
cmake_minimum_required(VERSION 2.8.7)
project(ani3d)
set(SOURCE bcg.f
matvecCSR.f
dsort.f
ilu0.f
iluoo.f)
@@ -0,0 +1,237 @@
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---->------------------------------------------------------------------<
@@ -0,0 +1,31 @@
C***********************************************************
C Perform Y:= BETA * Y + ALPHA * A * X
C ALPHA, BETA are real numbers,
C A is the matrix in sparse row format:
C N order of matrix
C A non-zero entries of the matrix in sparse row format
C jA indexes of columns for entries
C iA(1,..,n+1) pointers to jA for the first entry of each row
C***********************************************************
SUBROUTINE matvec( IMATVEC, ALPHA, X, BETA, Y, IA,JA,A )
IMPLICIT NONE

INTEGER IMATVEC(*), N
REAL*8 X(*), Y(*), ALPHA, BETA
Integer iA(*),jA(*)
REAL*8 A(*)

Integer i,j

N = IMATVEC(1)

do i = 1, N
Y(i) = BETA*Y(i)
do j = IA(i),IA(i+1)-1
Y(i) = Y(i) + A(j)*X(JA(j))*ALPHA
end do
end do

return
end

@@ -1937,7 +1937,7 @@ namespace INMOST
VectorPreallocateAni(solution_data,local_size);
{
VectorFillAni(rhs_data,&RHS[vbeg]);
VectorFinalizeAni(rhs_data);s
VectorFinalizeAni(rhs_data);

VectorFillAni(solution_data,&SOL[vbeg]);
VectorFinalizeAni(solution_data);
@@ -6,13 +6,13 @@
#if defined(USE_SOLVER_ANI)
void MatrixCopyDataAni(void ** ppA, void * pB);
void MatrixAssignDataAni(void * pA, void * pB);
void MatrixInitDataAni(void ** ppA, MSPP_MPI_Comm comm, const char * name);
void MatrixInitDataAni(void ** ppA, INMOST_MPI_Comm comm, const char * name);
void MatrixDestroyDataAni(void ** pA);
void MatrixFillAni(void * pA, int size, int * ia, int * ja, double * values);
void MatrixFillValuesAni(void * pA, double * values);
void MatrixFinalizeAni(void * data);

void VectorInitDataAni(void ** ppA, MSPP_MPI_Comm comm, const char * name);
void VectorInitDataAni(void ** ppA, INMOST_MPI_Comm comm, const char * name);
void VectorCopyDataAni(void ** ppA, void * pB);
void VectorAssignDataAni(void * pA, void * pB);
void VectorPreallocateAni(void * pA, int size);
@@ -25,8 +25,8 @@ void SolverInitializeAni(int * argc,char *** argv, const char * file_options);
bool SolverIsFinalizedAni();
void SolverFinalizeAni();
void SolverDestroyDataAni(void ** data);
void SolverInitDataAni(void ** data, MSPP_MPI_Comm comm, const char * name);
void SolverCopyDataAni(void ** data, void * other_data, MSPP_MPI_Comm comm);
void SolverInitDataAni(void ** data, INMOST_MPI_Comm comm, const char * name);
void SolverCopyDataAni(void ** data, void * other_data, INMOST_MPI_Comm comm);
void SolverAssignDataAni(void * data, void * other_data);
void SolverSetMatrixAni(void * data, void * matrix_data, bool same_pattern, bool reuse_preconditioner);
bool SolverSolveAni(void * data, void * rhs_data, void * sol_data);

0 comments on commit d13bcf9

Please sign in to comment.
You can’t perform that action at this time.