Skip to content
GitLab
Commit 86f9bc6f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Explicit Schur primal matrix creation for augmented FETI-DP finished.

parent 2eaf89db
Hide whitespace changes
Inline Side-by-side
AMDiS/src/parallel/PetscHelper.cc
View file @ 86f9bc6f
......@@ -11,6 +11,7 @@
#include "parallel/PetscHelper.h"
#include "Global.h"
namespace AMDiS {
......@@ -70,6 +71,104 @@ namespace AMDiS {
VecAssemblyBegin(vec);
VecAssemblyEnd(vec);
}
void blockMatMatSolve(KSP ksp, Mat mat0, Mat &mat1)
{
// === We have to calculate mat1 = ksp mat0: ===
// === - get all local column vectors from mat0 ===
// === - solve with ksp for this column vector as the rhs vector ===
// === - set the result to the corresponding column vector of ===
// === matrix mat1 ===
PetscInt localRows, localCols, globalRows, globalCols;
MatGetLocalSize(mat0, &localRows, &localCols);
MatGetSize(mat0, &globalRows, &globalCols);
MatCreateAIJ(PETSC_COMM_WORLD,
localRows, localCols, globalRows, globalCols,
150, PETSC_NULL, 150, PETSC_NULL, &mat1);
MatSetOption(mat1, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
// Transform matrix mat0 into a sparse column format.
PetscMatCol mat0_cols;
getMatLocalColumn(mat0, mat0_cols);
Vec tmpVec;
VecCreateSeq(PETSC_COMM_SELF, localRows, &tmpVec);
// Solve for all column vectors of mat A_BPi and create matrix matK
for (PetscMatCol::iterator it = mat0_cols.begin();
it != mat0_cols.end(); ++it) {
getColumnVec(it->second, tmpVec);
KSPSolve(ksp, tmpVec, tmpVec);
setMatLocalColumn(mat1, it->first, tmpVec);
}
VecDestroy(&tmpVec);
MatAssemblyBegin(mat1, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat1, MAT_FINAL_ASSEMBLY);
}
void matNestConvert(Mat matNest, Mat &mat)
{
FUNCNAME("matNestConvert()");
PetscInt nestRows, nestCols;
MatNestGetSize(matNest, &nestRows, &nestCols);
TEST_EXIT(nestRows == 2 && nestCols == 2)
("This funciton is only implemented for 2x2 nested matrices!\n");
Mat mat00, mat01, mat10, mat11;
MatNestGetSubMat(matNest, 0, 0, &mat00);
MatNestGetSubMat(matNest, 0, 1, &mat01);
MatNestGetSubMat(matNest, 1, 0, &mat10);
MatNestGetSubMat(matNest, 1, 1, &mat11);
int nRankRows0, nOverallRows0;
MatGetLocalSize(mat00, &nRankRows0, PETSC_NULL);
MatGetSize(mat00, &nOverallRows0, PETSC_NULL);
int nRankRows1, nOverallRows1;
MatGetLocalSize(mat11, &nRankRows1, PETSC_NULL);
MatGetSize(mat11, &nOverallRows1, PETSC_NULL);
int firstRow0;
MatGetOwnershipRange(mat00, &firstRow0, PETSC_NULL);
int firstRow1;
MatGetOwnershipRange(mat11, &firstRow1, PETSC_NULL);
int nRankRows = nRankRows0 + nRankRows1;
int nOverallRows = nOverallRows0 + nOverallRows1;
int firstRow = firstRow0 + firstRow1;
MatCreateAIJ(PETSC_COMM_WORLD,
nRankRows, nRankRows,
nOverallRows, nOverallRows,
10, PETSC_NULL, 10, PETSC_NULL, &mat);
for (int i = 0; i < nRankRows0; i++) {
PetscInt nCols;
const PetscInt *cols;
const PetscScalar *vals;
MatGetRow(mat00, i + firstRow0, &nCols, &cols, &vals);
for (int j = 0; j < nCols; j++) {
MatSetValue(mat, firstRow + i, cols[j], vals[j], INSERT_VALUES);
}
MatRestoreRow(mat00, i + firstRow0, &nCols, &cols, &vals);
}
MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY);
}
}
}
AMDiS/src/parallel/PetscHelper.h
View file @ 86f9bc6f
......@@ -69,6 +69,26 @@ namespace AMDiS {
* \param[out] vec Vector representing one column of the matrix.
*/
void getColumnVec(const SparseCol &matCol, Vec vec);
/** \brief
* Computes the matrix matrix product inv(A) B = C. Matrices B and C
* are distributed matrices. Matrix A is a local matrix on each MPI
* task. The overall number of rows of local matrices A must be the
* number of distriubted rows in B.
*
* \param[in] ksp inv(A) matrix given by a PETSc solver object.
* \param[in] mat0 matrix B
* \param[out] mat1 resulting matrix C, is created inside the function
*/
void blockMatMatSolve(KSP ksp, Mat mat0, Mat &mat1);
/** \brief
* Converts a 2x2 nested matrix to a MATAIJ matrix (thus not nested).
*
* \param[in] matNest nested input matrix
* \param[out] mat matrix of type MATAIJ, created inside this function.
*/
void matNestConvert(Mat matNest, Mat &mat);
}
}
......
AMDiS/src/parallel/PetscSolverFeti.cc
View file @ 86f9bc6f
......@@ -794,6 +794,7 @@ namespace AMDiS {
(void(*)(void))petscMultMatSchurPrimal);
} else {
schurPrimalAugmentedData.subSolver = subdomain;
schurPrimalAugmentedData.nestedVec = true;
localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
......@@ -824,63 +825,16 @@ namespace AMDiS {
double wtime = MPI::Wtime();
TEST_EXIT(!augmentedLagrange)("Not yet supported!\n");
TEST_EXIT_DBG(meshLevel == 0)
("Does not support for multilevel, check usage of localDofMap.\n");
// === First, calculate matK = inv(A_BB) A_BPi: ===
// === - get all local column vectors from A_BPi ===
// === - solve with A_BB for this column vector as the rhs vector ===
// === - set the result to the corresponding column vector of ===
// === matrix matK ===
// === Create explicit matrix representation of the Schur primal system. ===
int nRowsRankPrimal = primalDofMap.getRankDofs();
int nRowsOverallPrimal = primalDofMap.getOverallDofs();
int nRowsRankB = localDofMap.getRankDofs();
// Transform matrix A_BPi into a sparse column format.
petsc_helper::PetscMatCol mat_b_primal_cols;
petsc_helper::getMatLocalColumn(subdomain->getMatInteriorCoarse(),
mat_b_primal_cols);
TEST_EXIT(static_cast<int>(mat_b_primal_cols.size()) ==
primalDofMap.getLocalDofs())
("Should not happen!\n");
Vec tmpVec;
VecCreateSeq(PETSC_COMM_SELF, nRowsRankB, &tmpVec);
Mat matK;
MatCreateAIJ(mpiCommGlobal,
nRowsRankB, nRowsRankPrimal,
nGlobalOverallInterior, nRowsOverallPrimal,
150, PETSC_NULL, 150, PETSC_NULL, &matK);
MatSetOption(matK, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
// Solve for all column vectors of mat A_BPi and create matrix matK
for (petsc_helper::PetscMatCol::iterator it = mat_b_primal_cols.begin();
it != mat_b_primal_cols.end(); ++it) {
petsc_helper::getColumnVec(it->second, tmpVec);
subdomain->solve(tmpVec, tmpVec);
petsc_helper::setMatLocalColumn(matK, it->first, tmpVec);
}
VecDestroy(&tmpVec);
MatAssemblyBegin(matK, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(matK, MAT_FINAL_ASSEMBLY);
// Calculate: mat_schur_primal = A_PiPi - A_PiB inv(A_BB) ABPi
// = A_PiPi - A_PiB matK
MatDuplicate(subdomain->getMatCoarse(), MAT_COPY_VALUES,
&mat_schur_primal);
Mat matPrimal;
MatMatMult(subdomain->getMatCoarseInterior(), matK, MAT_INITIAL_MATRIX,
PETSC_DEFAULT, &matPrimal);
MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
MatDestroy(&matPrimal);
MatDestroy(&matK);
if (!augmentedLagrange)
createMatExplicitSchurPrimal();
else
createMatExplicitAugmentedSchurPrimal();
// === Create KSP solver object and set appropriate solver options. ====
......@@ -919,6 +873,160 @@ namespace AMDiS {
}
void PetscSolverFeti::createMatExplicitSchurPrimal()
{
FUNCNAME("PetscSolverFeti::createMatExplicitSchurPrimal()");
int creationMode = 0;
Parameters::get("parallel->feti->schur primal creation mode", creationMode);
if (creationMode == 0) {
// matK = inv(A_BB) A_BPi
Mat matK;
petsc_helper::blockMatMatSolve(subdomain->getSolver(),
subdomain->getMatInteriorCoarse(),
matK);
// mat_schur_primal = A_PiPi - A_PiB inv(A_BB) A_BPi
// = A_PiPi - A_PiB matK
MatMatMult(subdomain->getMatCoarseInterior(), matK, MAT_INITIAL_MATRIX,
PETSC_DEFAULT, &mat_schur_primal);
MatAYPX(mat_schur_primal, -1.0, subdomain->getMatCoarse(), DIFFERENT_NONZERO_PATTERN);
MatDestroy(&matK);
} else {
Mat tmp;
schurPrimalData.subSolver = subdomain;
localDofMap.createVec(schurPrimalData.tmp_vec_b, nGlobalOverallInterior);
primalDofMap.createVec(schurPrimalData.tmp_vec_primal);
MatCreateShell(mpiCommGlobal,
primalDofMap.getRankDofs(),
primalDofMap.getRankDofs(),
primalDofMap.getOverallDofs(),
primalDofMap.getOverallDofs(),
&schurPrimalData,
&tmp);
MatShellSetOperation(tmp, MATOP_MULT,
(void(*)(void))petscMultMatSchurPrimal);
MatComputeExplicitOperator(tmp, &mat_schur_primal);
MatDestroy(&tmp);
schurPrimalData.subSolver = NULL;
VecDestroy(&schurPrimalData.tmp_vec_b);
VecDestroy(&schurPrimalData.tmp_vec_primal);
}
}
void PetscSolverFeti::createMatExplicitAugmentedSchurPrimal()
{
FUNCNAME("PetscSolverFeti::createMatExplicitAugmentedSchurPrimal()");
int creationMode = 1;
Parameters::get("parallel->feti->schur primal creation mode", creationMode);
if (creationMode == 0) {
// qj = -Q J
Mat qj;
MatMatMult(mat_augmented_lagrange, mat_lagrange, MAT_INITIAL_MATRIX,
PETSC_DEFAULT, &qj);
MatScale(qj, -1.0);
// matTmp = inv(A_BB) A_BPi
Mat matTmp;
petsc_helper::blockMatMatSolve(subdomain->getSolver(),
subdomain->getMatInteriorCoarse(),
matTmp);
// mat00 = A_PiPi - A_PiB inv(A_BB) A_BPi
// = A_PiPi - A_PiB matTmp
Mat mat00;
MatMatMult(subdomain->getMatCoarseInterior(), matTmp, MAT_INITIAL_MATRIX,
PETSC_DEFAULT, &mat00);
MatAYPX(mat00, -1.0, subdomain->getMatCoarse(), DIFFERENT_NONZERO_PATTERN);
// mat10 = -Q J inv(A_BB) A_BPi
// = qj matTmp
Mat mat10;
MatMatMult(qj, matTmp, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &mat10);
// matTmp = inv(A_BB) trans(J) trans(Q)
Mat qT, jTqT;
MatTranspose(mat_augmented_lagrange, MAT_INITIAL_MATRIX, &qT);
MatTransposeMatMult(mat_lagrange, qT, MAT_INITIAL_MATRIX, PETSC_DEFAULT,
&jTqT);
petsc_helper::blockMatMatSolve(subdomain->getSolver(), jTqT, matTmp);
MatDestroy(&qT);
MatDestroy(&jTqT);
// mat01 = -A_PiB inv(A_BB) trans(J) trans(Q)
// = -A_PiB matTmp
Mat mat01;
MatMatMult(subdomain->getMatCoarseInterior(), matTmp, MAT_INITIAL_MATRIX,
PETSC_DEFAULT, &mat01);
MatScale(mat01, -1.0);
// mat11 = -Q J inv(A_BB) trans(J) trans(Q)
// = qj matTmp
Mat mat11;
MatMatMult(qj, matTmp, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &mat11);
MatDestroy(&matTmp);
MatDestroy(&qj);
Mat nestMat[4] = {mat00, mat01, mat10, mat11};
MatCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, 2, PETSC_NULL, nestMat, &matTmp);
petsc_helper::matNestConvert(matTmp, mat_schur_primal);
MatDestroy(&mat00);
MatDestroy(&mat01);
MatDestroy(&mat10);
MatDestroy(&mat11);
MatDestroy(&matTmp);
} else {
Mat tmp;
schurPrimalAugmentedData.subSolver = subdomain;
schurPrimalAugmentedData.nestedVec = false;
localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b0, nGlobalOverallInterior);
localDofMap.createVec(schurPrimalAugmentedData.tmp_vec_b1, nGlobalOverallInterior);
primalDofMap.createVec(schurPrimalAugmentedData.tmp_vec_primal);
lagrangeMap.createVec(schurPrimalAugmentedData.tmp_vec_lagrange);
schurPrimalAugmentedData.mat_lagrange = &mat_lagrange;
schurPrimalAugmentedData.mat_augmented_lagrange = &mat_augmented_lagrange;
MatCreateShell(mpiCommGlobal,
primalDofMap.getRankDofs() + nRankEdges,
primalDofMap.getRankDofs() + nRankEdges,
primalDofMap.getOverallDofs() + nOverallEdges,
primalDofMap.getOverallDofs() + nOverallEdges,
&schurPrimalAugmentedData,
&tmp);
MatShellSetOperation(tmp, MATOP_MULT,
(void(*)(void))petscMultMatSchurPrimalAugmented);
MatComputeExplicitOperator(tmp, &mat_schur_primal);
MatDestroy(&tmp);
schurPrimalAugmentedData.subSolver = NULL;
schurPrimalAugmentedData.mat_lagrange = NULL;
schurPrimalAugmentedData.mat_augmented_lagrange = NULL;
VecDestroy(&schurPrimalAugmentedData.tmp_vec_b0);
VecDestroy(&schurPrimalAugmentedData.tmp_vec_b1);
VecDestroy(&schurPrimalAugmentedData.tmp_vec_primal);
VecDestroy(&schurPrimalAugmentedData.tmp_vec_lagrange);
}
}
void PetscSolverFeti::destroySchurPrimalKsp()
{
FUNCNAME("PetscSolverFeti::destroySchurPrimal()");
......
AMDiS/src/parallel/PetscSolverFeti.h
View file @ 86f9bc6f
......@@ -140,6 +140,12 @@ namespace AMDiS {
/// system on the primal variables, \ref ksp_schur_primal
void createSchurPrimalKsp(vector<const FiniteElemSpace*> &feSpaces);
///
void createMatExplicitSchurPrimal();
///
void createMatExplicitAugmentedSchurPrimal();
/// Destroys PETSc KSP solver object \ref ksp_schur_primal
void destroySchurPrimalKsp();
......
AMDiS/src/parallel/PetscSolverFetiOperators.cc
View file @ 86f9bc6f
......@@ -41,11 +41,50 @@ namespace AMDiS {
SchurPrimalAugmentedData* data =
static_cast<SchurPrimalAugmentedData*>(ctx);
Vec x_primal, x_mu, y_primal, y_mu;
VecNestGetSubVec(x, 0, &x_primal);
VecNestGetSubVec(x, 1, &x_mu);
VecNestGetSubVec(y, 0, &y_primal);
VecNestGetSubVec(y, 1, &y_mu);
Vec x_primal, x_mu, y_primal, y_mu;
if (data->nestedVec) {
VecNestGetSubVec(x, 0, &x_primal);
VecNestGetSubVec(x, 1, &x_mu);
VecNestGetSubVec(y, 0, &y_primal);
VecNestGetSubVec(y, 1, &y_mu);
} else {
VecDuplicate(data->tmp_vec_primal, &x_primal);
VecDuplicate(data->tmp_vec_primal, &y_primal);
PetscInt allLocalSize, allSize;
VecGetLocalSize(x, &allLocalSize);
VecGetSize(x, &allSize);
PetscInt primalLocalSize, primalSize;
VecGetLocalSize(x_primal, &primalLocalSize);
VecGetSize(x_primal, &primalSize);
TEST_EXIT_DBG(allSize > primalSize)("Should not happen!\n");
TEST_EXIT_DBG(allLocalSize >= primalLocalSize)("Should not happen!\n");
PetscInt muLocalSize = allLocalSize - primalLocalSize;
PetscInt muSize = allSize - primalSize;
VecCreateMPI(PETSC_COMM_WORLD, muLocalSize, muSize, &x_mu);
VecCreateMPI(PETSC_COMM_WORLD, muLocalSize, muSize, &y_mu);
PetscScalar *allValue;
PetscScalar *primalValue;
PetscScalar *muValue;
VecGetArray(x, &allValue);
VecGetArray(x_primal, &primalValue);
VecGetArray(x_mu, &muValue);
for (int i = 0; i < primalLocalSize; i++)
primalValue[i] = allValue[i];
for (int i = 0; i < muLocalSize; i++)
muValue[i] = allValue[primalLocalSize + i];
VecRestoreArray(x, &allValue);
VecRestoreArray(x_primal, &primalValue);
VecRestoreArray(x_mu, &muValue);
}
// inv(K_BB) K_BPi x_Pi
MatMult(data->subSolver->getMatInteriorCoarse(), x_primal, data->tmp_vec_b0);
......@@ -74,6 +113,36 @@ namespace AMDiS {
MatMult(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, y_mu);
VecScale(y_mu, -1.0);
if (!data->nestedVec) {
PetscInt allLocalSize;
VecGetLocalSize(x, &allLocalSize);
PetscInt primalLocalSize;
VecGetLocalSize(x_primal, &primalLocalSize);
PetscInt muLocalSize = allLocalSize - primalLocalSize;
PetscScalar *allValue;
PetscScalar *primalValue;
PetscScalar *muValue;
VecGetArray(y, &allValue);
VecGetArray(y_primal, &primalValue);
VecGetArray(y_mu, &muValue);
for (int i = 0; i < primalLocalSize; i++)
allValue[i] = primalValue[i];
for (int i = 0; i < muLocalSize; i++)
allValue[primalLocalSize + i] = muValue[i];
VecRestoreArray(y, &allValue);
VecRestoreArray(y_primal, &primalValue);
VecRestoreArray(y_mu, &muValue);
VecDestroy(&x_primal);
VecDestroy(&y_primal);
VecDestroy(&x_mu);
VecDestroy(&y_mu);
}
return 0;
}
......
AMDiS/src/parallel/PetscSolverFetiStructs.h
View file @ 86f9bc6f
......@@ -63,6 +63,8 @@ namespace AMDiS {
Mat *mat_augmented_lagrange;
PetscSolver* subSolver;
bool nestedVec;
};
......
Supports Markdown
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