Commit 53f9c8c6 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

petsc : use BAIJ matrices (faster MatSetValues)

parent 68db8d4e
......@@ -19,11 +19,11 @@ struct HXTLinearSystemPETScStruct {
uint32_t nDofs;
int nFields;
uint32_t *nodeMap;
uint32_t *inverseNodeMap;
double *fixed;
KSP ksp;
Mat a;
Vec x;
uint32_t nIdentRows, nIdentRowsAlloc;
PetscInt *identRows;
int assemblyNeeded;
};
......@@ -67,6 +67,12 @@ static HXTStatus assembleIfNeeded(HXTLinearSystemPETSc *system)
HXT_PETSC_CHECK(MatAssemblyEnd(system->a, MAT_FINAL_ASSEMBLY));
system->assemblyNeeded = 0;
}
if (system->nIdentRows > 0){
printf("%i identity rows\n", system->nIdentRows);
HXT_PETSC_CHECK(MatZeroRows(system->a,
system->nIdentRows, system->identRows,1., NULL, NULL));
}
system->nIdentRows = 0;
return HXT_STATUS_OK;
}
......@@ -89,7 +95,7 @@ static HXTStatus linearSystemPreAllocateMatrix(HXTLinearSystemPETSc *system)
}
}
}
int *nByRow = malloc(sizeof(int) * system->nMappedNodes*system->nFields);
int *nByRow = malloc(sizeof(int) * system->nMappedNodes);
for (uint32_t i = 0; i < system->nMappedNodes; ++i){
int n = 0;
SparsityChunk *c = sparsity[i];
......@@ -101,34 +107,26 @@ static HXTStatus linearSystemPreAllocateMatrix(HXTLinearSystemPETSc *system)
if(c->n[j] != -1)
n++;
sparsityChunkFree(sparsity[i]);
for (uint32_t j = 0; j < system->nFields; ++j){
nByRow[i*system->nFields+j] = n*system->nFields;
}
nByRow[i] = n;
}
free(sparsity);
HXT_PETSC_CHECK(MatSeqAIJSetPreallocation(system->a, 0, nByRow));
HXT_PETSC_CHECK(MatSeqBAIJSetPreallocation(system->a, system->nFields,0, nByRow));
free(nByRow);
return HXT_STATUS_OK;
}
HXTStatus hxtLinearSystemPETScAddToMatrix(HXTLinearSystemPETSc *system, int el0, int el1, const double *localMatrix){
int nf = system->nFields;
int nn = system->nNodesByElement;
uint32_t *e0 = system->elements + el0*nn;
uint32_t *e1 = system->elements + el1*nn;
int *map0 = malloc(sizeof(int)*nn*nf);
int *map1 = malloc(sizeof(int)*nn*nf);
int *map0 = malloc(sizeof(int)*nn);
int *map1 = malloc(sizeof(int)*nn);
for(int i=0; i<nn; ++i){
for (int j=0; j<nf; ++j){
//map0[j*nn+i] = system->nodeMap[e0[i]]*nf+j;
//map1[j*nn+i] = system->nodeMap[e1[i]]*nf+j;
map0[i*nf+j] = system->nodeMap[e0[i]]*nf+j;
map1[i*nf+j] = system->nodeMap[e1[i]]*nf+j;
}
map0[i] = system->nodeMap[e0[i]];
map1[i] = system->nodeMap[e1[i]];
}
//HXT_PETSC_CHECK(MatSetValues(system->a,nn*nf,map0,nn*nf,map0,localMatrix,ADD_VALUES));
HXT_PETSC_CHECK(MatSetValues(system->a,nn*nf,map0,nn*nf,map1,localMatrix,ADD_VALUES));
HXT_PETSC_CHECK(MatSetValuesBlocked(system->a,nn,map0,nn,map1,localMatrix,ADD_VALUES));
free(map0);
free(map1);
system->assemblyNeeded = 1;
......@@ -162,18 +160,12 @@ HXTStatus hxtLinearSystemPETScSetMatrixRowIdentity(HXTLinearSystemPETSc *system,
HXT_WARNING("ignoring boundary condition on node %i", node);
return HXT_STATUS_OK;
}
HXT_CHECK(assembleIfNeeded(system));
int row = system->nodeMap[node]*system->nFields + field;
int ncols;
const int *cols;
const double *v;
HXT_PETSC_CHECK(MatGetRow(system->a, row, &ncols, &cols, &v));
for (int i = 0; i < ncols; ++i) {
// I know this is not correct but it works and using MatSetValues or MatZeroRows is much slower
((double*)v)[i] = cols[i] == row ? 1. : 0.;
if(system->nIdentRows == system->nIdentRowsAlloc) {
system->nIdentRowsAlloc *= 2;
system->identRows = realloc(system->identRows, sizeof(PetscInt)*system->nIdentRowsAlloc);
}
HXT_PETSC_CHECK(MatRestoreRow(system->a, row, &ncols, &cols, &v));
//HXT_PETSC_CHECK(MatZeroRows(system->a, 1, &row,1., NULL, NULL));
system->identRows[system->nIdentRows++] = row;
return HXT_STATUS_OK;
}
......@@ -230,11 +222,7 @@ HXTStatus hxtLinearSystemPETScAddRhsEntry(HXTLinearSystemPETSc *system, double *
HXTStatus hxtLinearSystemPETScSolve(HXTLinearSystemPETSc *system, double *rhs, double *solution){
Vec b;
HXT_PETSC_CHECK(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, system->nDofs, rhs, &b));
if(system->assemblyNeeded) {
HXT_PETSC_CHECK(MatAssemblyBegin(system->a, MAT_FINAL_ASSEMBLY));
HXT_PETSC_CHECK(MatAssemblyEnd(system->a, MAT_FINAL_ASSEMBLY));
system->assemblyNeeded = 0;
}
assembleIfNeeded(system);
HXT_PETSC_CHECK(KSPSetOperators(system->ksp, system->a, system->a));
HXT_PETSC_CHECK(KSPSolve(system->ksp, b, system->x));
HXT_PETSC_CHECK(VecDestroy(&b));
......@@ -264,6 +252,9 @@ HXTStatus hxtLinearSystemPETScCreate(HXTLinearSystemPETSc **pSystem, int nElemen
{
HXTLinearSystemPETSc *system = malloc(sizeof(HXTLinearSystemPETSc));
*pSystem = system;
system->nIdentRows = 0;
system->nIdentRowsAlloc = 256;
system->identRows = malloc (sizeof(uint32_t)*system->nIdentRowsAlloc);
system->nElements = nElements;
system->nNodesByElement = nNodesByElement;
system->nFields = nFields;
......@@ -288,6 +279,8 @@ HXTStatus hxtLinearSystemPETScCreate(HXTLinearSystemPETSc **pSystem, int nElemen
uint32_t nDofs = nMappedNodes * nFields;
system->nDofs = nDofs;
HXT_PETSC_CHECK(MatCreate(MPI_COMM_WORLD, &system->a));
HXT_PETSC_CHECK(MatSetType(system->a,MATBAIJ));
HXT_PETSC_CHECK(MatSetSizes(system->a, nDofs, nDofs, PETSC_DETERMINE, PETSC_DETERMINE));
if (prefix != NULL)
HXT_PETSC_CHECK(MatSetOptionsPrefix(system->a, prefix));
......@@ -316,6 +309,7 @@ HXTStatus hxtLinearSystemPETScDelete(HXTLinearSystemPETSc **pSystem)
HXT_PETSC_CHECK(KSPDestroy(&system->ksp));
HXT_PETSC_CHECK(VecDestroy(&system->x));
HXT_PETSC_CHECK(MatDestroy(&system->a));
free(system->identRows);
free(system->nodeMap);
free(system);
*pSystem = NULL;
......
Markdown is supported
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