Skip to content
Snippets Groups Projects
Commit 545d04f0 authored by Ondrej Meca's avatar Ondrej Meca
Browse files

ENH: new version of FETI4I API

parent 63f5f507
Branches master
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 10 additions and 448 deletions
fem/src/DirectSolve.F90
+ 10
− 3
View file @ 545d04f0
...@@ -1883,8 +1883,8 @@ CONTAINS ...@@ -1883,8 +1883,8 @@ CONTAINS
INTEGER :: n_dof_partition INTEGER :: n_dof_partition
INTEGER, ALLOCATABLE :: memb(:), DirichletInds(:), Neighbours(:) INTEGER, ALLOCATABLE :: memb(:), DirichletInds(:), Neighbours(:), IntOptions(:)
REAL(KIND=dp), ALLOCATABLE :: DirichletVals(:) REAL(KIND=dp), ALLOCATABLE :: DirichletVals(:), RealOptions(:)
INTEGER :: Comm_active, Group_active, Group_world INTEGER :: Comm_active, Group_active, Group_world
REAL(KIND=dp), ALLOCATABLE :: dbuf(:) REAL(KIND=dp), ALLOCATABLE :: dbuf(:)
...@@ -1950,10 +1950,17 @@ CONTAINS ...@@ -1950,10 +1950,17 @@ CONTAINS
'inconsistency: A % NumberOfRows /= SIZE(A % ParallelInfo % GlobalDOFs' ) 'inconsistency: A % NumberOfRows /= SIZE(A % ParallelInfo % GlobalDOFs' )
END IF END IF
ALLOCATE(IntOptions(10))
ALLOCATE(RealOptions(1))
CALL FETI4ISetDefaultIntegerOptions(IntOptions)
CALL FETI4ISetDefaultRealOptions(RealOptions)
CALL FETI4ICreateInstance(A % PermonSolverInstance, A % PermonMatrix, & CALL FETI4ICreateInstance(A % PermonSolverInstance, A % PermonMatrix, &
A % NumberOfRows, b, A % ParallelInfo % GlobalDOFs, & A % NumberOfRows, b, A % ParallelInfo % GlobalDOFs, &
n, neighbours, & n, neighbours, &
nd, DirichletInds, DirichletVals) nd, DirichletInds, DirichletVals, &
IntOptions, RealOptions)
END IF END IF
!CALL Permon_Solve( A % PermonSolverInstance, x, b ) !CALL Permon_Solve( A % PermonSolverInstance, x, b )
......
feti4i/feti4i.h deleted 100644 → 0
+ 0
− 134
View file @ 63f5f507
#ifndef FETI4I_H_
#define FETI4I_H_
#include "mpi.h"
/*-----------------------------------------------------------------------------
Set data-types used in FETI4I
Possible values for FETI4I_INT_WIDTH:
32: use 32 bit signed integer
64: use 64 bit signed integer
Possible values for FETI4I_REAL_WIDTH:
64: ESPRESO supports only 64 bit real values
------------------------------------------------------------------------------*/
#ifndef FETI4I_INT_WIDTH
#define FETI4I_INT_WIDTH 32
#endif
#ifndef FETI4I_REAL_WIDTH
#define FETI4I_REAL_WIDTH 64
#endif
#if FETI4I_INT_WIDTH == 32
typedef int FETI4IInt;
#elif FETI4I_INT_WIDTH == 64
typedef long FETI4IInt;
#else
#error "Incorrect user-supplied value of FETI4I_INT_WIDTH"
#endif
/* MPI integer (e.g. rank) are always 32-bit */
typedef int FETI4IMPIInt;
#if FETI4I_REAL_WIDTH == 64
typedef double FETI4IReal;
#else
#error "Incorrect user-supplied value of FETI4I_REAL_WIDTH"
#endif
#ifdef __cplusplus
extern "C" {
#endif
/*-----------------------------------------------------------------------------
Definitions of internal structures used in FETI4I
------------------------------------------------------------------------------*/
typedef struct FETI4IStructMatrix* FETI4IMatrix;
typedef struct FETI4IStructInstance* FETI4IInstance;
/*-----------------------------------------------------------------------------
Functions for manipulating with FETI4I internal structures
------------------------------------------------------------------------------*/
void FETI4ICreateStiffnessMatrix(
FETI4IMatrix *matrix, //TODO size?
FETI4IInt indexBase
);
void FETI4IAddElement( //TODO CRS or CCS ?
FETI4IMatrix matrix,
FETI4IInt size,
FETI4IInt* indices,
FETI4IReal* values
);
/*-----------------------------------------------------------------------------
Functions for creating an instance and solve it
------------------------------------------------------------------------------*/
void FETI4ICreateInstance(
FETI4IInstance *instance,
FETI4IMatrix matrix,
FETI4IInt size,
FETI4IReal* rhs,
FETI4IInt* l2g, /* length of both rhs and l2g is size */
FETI4IMPIInt neighbours_size,
FETI4IMPIInt* neighbours,
FETI4IInt dirichlet_size,
FETI4IInt* dirichlet_indices, //TODO which numbering? we prefer global numbering
FETI4IReal* dirichlet_values
);
void FETI4ISolve(
FETI4IInstance instance,
FETI4IInt solution_size,
FETI4IReal* solution
);
/*-----------------------------------------------------------------------------
Functions for updating a created instance
------------------------------------------------------------------------------*/
void FETI4IUpdateStiffnessMatrix(
FETI4IInstance instance,
FETI4IMatrix stiffnessMatrix
);
void FETI4IUpdateRhs(
FETI4IInstance instance,
FETI4IInt rhs_size,
FETI4IReal* rhs_values
);
//TODO VH: neighbours perhaps should not be passed here and they are FETI4IMPIInt
void FETI4IUpdateDirichlet(
FETI4IInstance instance,
FETI4IInt dirichlet_size,
FETI4IInt* dirichlet_indices,
FETI4IReal* dirichlet_values,
FETI4IReal* l2g,
FETI4IInt neighbour_size,
FETI4IReal* neighbour
);
/*-----------------------------------------------------------------------------
Destroy an arbitrary internal structure
------------------------------------------------------------------------------*/
void FETI4IDestroy(
void* ptr
);
#ifdef __cplusplus
}
#endif
#endif /* FETI4I_H_ */
feti4i/feti4i_fortran_def.h90 deleted 100644 → 0
+ 0
− 50
View file @ 63f5f507
!#ifndef FETI4I_FORTRAN_DEF_H_
!#define FETI4I_FORTRAN_DEF_H_
! ------------------------------------------------------------------------------
! Set data-types used in FETI4I
!
! Possible values for FETI4I_INT_WIDTH:
! 32: use 32 bit signed integer
! 64: use 64 bit signed integer
!
! Possible values for FETI4I_REAL_WIDTH:
! 64: ESPRESO supports only 64 bit real values
! ------------------------------------------------------------------------------
USE, INTRINSIC :: ISO_C_BINDING
#ifndef FETI4I_INT_WIDTH
#define FETI4I_INT_WIDTH 32
#endif
#ifndef FETI4I_REAL_WIDTH
#define FETI4I_REAL_WIDTH 64
#endif
#if FETI4I_INT_WIDTH == 32
INTEGER, PARAMETER :: FETI4IInt = C_INT
#elif FETI4I_INT_WIDTH == 64
INTEGER, PARAMETER :: FETI4IInt = C_LONG
#else
#error "Incorrect user-supplied value of FETI4I_INT_WIDTH"
#endif
! MPI integer (e.g. rank) are always 32-bit
INTEGER, PARAMETER :: FETI4IMPIInt = C_INT
#if FETI4I_REAL_WIDTH == 64
INTEGER, PARAMETER :: FETI4IReal = C_DOUBLE
#else
#error "Incorrect user-supplied value of FETI4I_REAL_WIDTH"
#endif
! -----------------------------------------------------------------------------
! Definitions of internal structures used in FETI4I
! -----------------------------------------------------------------------------
! use TYPE(C_PTR) for opaque pointers
!TODO we could perhaps use Fortran types
!#endif /* FETI4I_FORTRAN_DEF_H_ */
feti4i/feti4i_fortran_interface.h90 deleted 100644 → 0
+ 0
− 52
View file @ 63f5f507
#include "feti4i_fortran_def.h90"
INTERFACE
SUBROUTINE FETI4ICreateStiffnessMatrix(matrix, indexBase) &
BIND(c,NAME='FETI4ICreateStiffnessMatrix')
IMPORT
TYPE(C_PTR) :: matrix ! FETI4IMatrix* type in C (output)
INTEGER(KIND=FETI4IInt), VALUE :: indexBase
END SUBROUTINE FETI4ICreateStiffnessMatrix
SUBROUTINE FETI4IAddElement(matrix, size, indices, values) &
BIND(c,NAME='FETI4IAddElement')
IMPORT
TYPE(C_PTR), VALUE :: matrix ! FETI4IMatrix type in C (input)
INTEGER(KIND=FETI4IInt), VALUE :: size
INTEGER(KIND=FETI4IInt) :: indices(*)
REAL(KIND=FETI4IReal) :: values(*)
END SUBROUTINE FETI4IAddElement
SUBROUTINE FETI4ICreateInstance(instance, matrix, size, rhs, l2g, &
neighbours_size, neighbours, &
dirichlet_size, dirichlet_indices, dirichlet_values) &
BIND(c,NAME='FETI4ICreateInstance')
IMPORT
TYPE(C_PTR) :: instance ! FETI4IInstance* type in C (output)
TYPE(C_PTR), VALUE :: matrix ! FETI4IMatrix type in C (input)
INTEGER(KIND=FETI4IInt), VALUE :: size
REAL(KIND=FETI4IReal) :: rhs(*)
INTEGER(KIND=FETI4IInt) :: l2g(*)
INTEGER(KIND=FETI4IMPIInt), VALUE :: neighbours_size
INTEGER(KIND=FETI4IMPIInt) :: neighbours(*)
INTEGER(KIND=FETI4IInt), VALUE :: dirichlet_size
INTEGER(KIND=FETI4IInt) :: dirichlet_indices(*)
REAL(KIND=FETI4IReal) :: dirichlet_values(*)
END SUBROUTINE FETI4ICreateInstance
SUBROUTINE FETI4ISolve(instance, solution_size, solution) &
BIND(c,NAME='FETI4ISolve')
IMPORT
TYPE(C_PTR), VALUE :: instance ! FETI4IInstance type in C (input)
INTEGER(KIND=FETI4IInt), VALUE :: solution_size
REAL(KIND=FETI4IReal) :: solution(*)
END SUBROUTINE FETI4ISolve
SUBROUTINE FETI4IDestroy(ptr) &
BIND(c,NAME='FETI4IDestroy')
IMPORT
TYPE(C_PTR), VALUE :: ptr ! void* type in C
END SUBROUTINE FETI4IDestroy
END INTERFACE
feti4i/feti4i_fortran_test/Makefile deleted 100644 → 0
+ 0
− 35
View file @ 63f5f507
FC := mpif90
CC := mpicc
MKDIR := mkdir -p
OUTLIBDIR := ../lib
build: test
feti4i_mod.o: ../feti4i_mod.F90
${FC} -c -fPIC $^
feti4i.mod: ../feti4i_mod.F90 feti4i_mod.o
@true
feti4i_dummy.o: feti4i_dummy.c
${CC} -c -fPIC -g $^
${OUTLIBDIR}/libfeti4i.so: feti4i_dummy.o
${MKDIR} ${OUTLIBDIR}
${LINK.F} -o ${OUTLIBDIR}/libfeti4i.so feti4i_dummy.o -shared -fPIC
test.o: test.F90 feti4i.mod
${FC} -c -fPIC -g test.F90 -o $@
test: test.o ${OUTLIBDIR}/libfeti4i.so
${LINK.F} -o test test.o -L${OUTLIBDIR} -Wl,-rpath=${OUTLIBDIR} -lfeti4i
clean:
${RM} test.o feti4i_dummy.o feti4i_mod.o
clean-all: clean
${RM} feti4i.mod test
${RM} -r ${OUTLIBDIR}
.PHONY: build clean
feti4i/feti4i_fortran_test/feti4i_dummy.c deleted 100644 → 0
+ 0
− 95
View file @ 63f5f507
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include "../feti4i.h"
struct FETI4IStructMatrix {
FETI4IInt indexBase;
FETI4IInt size;
FETI4IInt* indices;
FETI4IReal* values;
FETI4IInt nelems;
};
struct FETI4IStructInstance {
FETI4IMatrix matrix;
FETI4IReal* rhs;
FETI4IInt* l2g; /* length of both rhs and l2g is size */
FETI4IMPIInt neighbours_size;
FETI4IMPIInt* neighbours;
FETI4IInt dirichlet_size;
FETI4IInt* dirichlet_indices; //TODO which numbering? we prefer global numbering
FETI4IReal* dirichlet_values;
};
void FETI4ICreateStiffnessMatrix(FETI4IMatrix *matrix, FETI4IInt indexBase)
{
*matrix = (FETI4IMatrix) malloc(sizeof(struct FETI4IStructMatrix));
(*matrix)->indexBase = indexBase;
(*matrix)->size = 1;
(*matrix)->indices = NULL;
(*matrix)->values = NULL;
(*matrix)->nelems = 0;
printf("%s: created matrix %p, indexBase %d\n", __func__, (void*) *matrix, indexBase);
}
void FETI4IAddElement(FETI4IMatrix matrix, FETI4IInt size, FETI4IInt* indices, FETI4IReal* values)
{
matrix->size += (size-1);
matrix->indices = indices;
matrix->values = values;
printf("%s: added element %d to matrix %p\n", __func__, matrix->nelems++, (void*) matrix);
}
void FETI4ICreateInstance(FETI4IInstance *instance, FETI4IMatrix matrix, FETI4IInt size,
FETI4IReal* rhs, FETI4IInt* l2g, FETI4IMPIInt neighbours_size, FETI4IMPIInt* neighbours,
FETI4IInt dirichlet_size, FETI4IInt* dirichlet_indices, FETI4IReal* dirichlet_values)
{
printf("%s: setting rhs %p (%d), l2g %p (%d), "
"neighbours %p (%d), dirichlet_indices %p (%d), "
"dirichlet_values %p (%d)\n",
__func__, (void*) rhs, size, (void*) l2g, size,
(void*) neighbours, neighbours_size, (void*) dirichlet_indices, dirichlet_size,
(void*) dirichlet_indices, dirichlet_size);
*instance = (FETI4IInstance) malloc(sizeof(struct FETI4IStructInstance));
(*instance)->matrix = matrix;
//assert(size == matrix->size);
(*instance)->rhs = rhs;
(*instance)->l2g = l2g;
(*instance)->neighbours_size = neighbours_size;
(*instance)->neighbours = neighbours;
(*instance)->dirichlet_size = dirichlet_size;
(*instance)->dirichlet_indices = dirichlet_indices;
(*instance)->dirichlet_values = dirichlet_values;
}
void FETI4ISolve(FETI4IInstance instance, FETI4IInt solution_size, FETI4IReal* solution)
{
FETI4IInt i;
FETI4IMatrix matrix = instance->matrix;
FETI4IReal *rhs = instance->rhs;
printf("%s: solving with matrix %p, rhs %p\n", __func__, (void*) instance->matrix, (void*) instance->rhs);
//assert(solution_size == matrix->size);
//for (i=0; i<solution_size; i++) {
// printf("#%d %g\n",i,rhs[i]);
//}
//printf("\n");
//for (i=0; i<solution_size; i++) {
// printf("#%d %g ",i,solution[i]);
// solution[i]=-i;
// printf(" %g\n",solution[i]);
//}
}
void FETI4IDestroy(void *ptr)
{
FETI4IInstance instance = (FETI4IInstance) ptr;
free(instance->matrix);
free(instance);
}
feti4i/feti4i_fortran_test/test.F90 deleted 100644 → 0
+ 0
− 74
View file @ 63f5f507
PROGRAM test
USE FETI4I
TYPE(C_PTR) :: instance
TYPE(C_PTR) :: matrix
INTEGER(KIND=FETI4IInt) :: n = 4
INTEGER(KIND=FETI4IInt) :: nelems = 5
INTEGER(KIND=FETI4IInt) :: size = 0 !TODO when this is known
INTEGER(KIND=FETI4IInt), ALLOCATABLE :: inds(:)
REAL(KIND=FETI4IReal), ALLOCATABLE :: vals(:)
REAL(KIND=FETI4IReal), ALLOCATABLE :: b(:)
INTEGER(KIND=FETI4IInt), ALLOCATABLE :: l2g(:)
INTEGER(KIND=FETI4IMPIInt) :: neighbours_size = 3
INTEGER(KIND=FETI4IMPIInt), ALLOCATABLE :: neighbours(:)
INTEGER(KIND=FETI4IInt) :: dirichlet_size = 2
INTEGER(KIND=FETI4IInt), ALLOCATABLE :: dirichlet_indices(:)
REAL(KIND=FETI4IReal), ALLOCATABLE :: dirichlet_values(:)
REAL(KIND=FETI4IReal), ALLOCATABLE :: x(:)
INTEGER :: ALLOC_ERR
CALL FETI4ICreateStiffnessMatrix(matrix, 1)
ALLOCATE(vals(n*n), inds(n))
DO i=1,n
DO j=1,n
vals((i-1)*n+j) = (i-1)*n+j
END DO
inds(i) = i
END DO
size = 1
DO i=1,nelems
CALL FETI4IAddElement(matrix, n, inds, vals)
size = size + (n-1)
END DO
ALLOCATE(b(size), l2g(size))
DO i=1,size
b(i) = -i
l2g(i) = 100+i
END DO
write (*,*) "b: ", b
ALLOCATE(neighbours(neighbours_size))
DO i=1,neighbours_size
neighbours(i) = i-1
END DO
ALLOCATE(dirichlet_indices(dirichlet_size), dirichlet_values(dirichlet_size))
DO i=1,dirichlet_size
dirichlet_indices(i) = 2*(i+1)
dirichlet_values(i) = -2*(i+1)
END DO
CALL FETI4ICreateInstance(instance, matrix, size, b, l2g, &
neighbours_size, neighbours, &
dirichlet_size, dirichlet_indices, dirichlet_values)
ALLOCATE(x(size), STAT = ALLOC_ERR)
write (*,*) "ALLOC_ERR=",ALLOC_ERR
x(:) = 1.0
write (*,*) "x: ", x
CALL FETI4ISolve(instance, size, x)
CALL FETI4IDestroy(instance)
END PROGRAM test
feti4i/feti4i_mod.F90 deleted 100644 → 0
+ 0
− 5
View file @ 63f5f507
MODULE feti4i
#include "feti4i_fortran_interface.h90"
END MODULE feti4i
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment