// // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology // All rights reserved. // Authors: Simon Vey, Thomas Witkowski et al. // // This file is part of AMDiS // // See also license.opensource.txt in the distribution. #include "AMDiS.h" #include "parallel/PetscSolverGlobalMatrix.h" #include "parallel/StdMpi.h" #include "parallel/MpiHelper.h" namespace AMDiS { void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix *seqMat) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()"); TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n"); TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n"); TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n"); double wtime = MPI::Wtime(); createMatVec(*seqMat); if (coarseSpaceMap.size()) { fillPetscMatrixWithCoarseSpace(seqMat); return; } // === Create PETSc vector (solution and a temporary vector). === #if (DEBUG != 0) MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif // === Transfer values from DOF matrices to the PETSc matrix. === int nComponents = seqMat->getNumRows(); for (int i = 0; i < nComponents; i++) for (int j = 0; j < nComponents; j++) if ((*seqMat)[i][j]) setDofMatrix((*seqMat)[i][j], i, j); #if (DEBUG != 0) MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif matAssembly(); if (printMatInfo) { MatInfo matInfo; MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo); MSG("Matrix info:\n"); MSG(" memory usage: %e MB\n", matInfo.memory / (1024.0 * 1024.0)); MSG(" mallocs: %d\n", static_cast(matInfo.mallocs)); MSG(" nz allocated: %d\n", static_cast(matInfo.nz_allocated)); MSG(" nz used: %d\n", static_cast(matInfo.nz_used)); MSG(" nz unneeded: %d\n", static_cast(matInfo.nz_unneeded)); } // === Init PETSc solver. === KSPCreate(mpiCommGlobal, &kspInterior); KSPGetPC(kspInterior, &pcInterior); KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT); KSPSetType(kspInterior, KSPBCGS); KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str()); KSPSetFromOptions(kspInterior); initPreconditioner(pcInterior); // Do not delete the solution vector, use it for the initial guess. if (!zeroStartVector) KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE); removeDirichletRows(seqMat, dofMap, getMatInterior()); #if (DEBUG != 0) MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif // === For debugging allow to write the matrix to a file. === bool dbgWriteMatrix = false; Parameters::get("parallel->debug->write matrix", dbgWriteMatrix); if (dbgWriteMatrix) { PetscViewer matView; PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.mat", FILE_MODE_WRITE, &matView); MatView(getMatInterior(), matView); PetscViewerDestroy(&matView); } } void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix *seqMat) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()"); TEST_EXIT_DBG(interiorMap)("Should not happen!\n"); TEST_EXIT_DBG(coarseSpaceMap.size() == seqMat->getSize()) ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize()); // === Prepare traverse of sequentially created matrices. === using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; vector cols, colsOther; vector values, valuesOther; cols.reserve(300); colsOther.reserve(300); values.reserve(300); valuesOther.reserve(300); // === Traverse all sequentially created matrices and add the values to === // === the global PETSc matrices. === int nComponents = seqMat->getSize(); for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) { for (int colComponent = 0; colComponent < nComponents; colComponent++) { DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent]; if (!dofMat) continue; ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent]; ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent]; traits::col::type col(dofMat->getBaseMatrix()); traits::const_value::type value(dofMat->getBaseMatrix()); // Traverse all rows. for (cursor_type cursor = begin(dofMat->getBaseMatrix()), cend = end(dofMat->getBaseMatrix()); cursor != cend; ++cursor) { bool isRowCoarse = isCoarseSpace(rowComponent, *cursor); cols.clear(); colsOther.clear(); values.clear(); valuesOther.clear(); // Traverse all columns. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { bool isColCoarse = isCoarseSpace(colComponent, col(*icursor)); if (isColCoarse == false) if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false) continue; if (isColCoarse == isRowCoarse) { cols.push_back(col(*icursor)); values.push_back(value(*icursor)); } else { colsOther.push_back(col(*icursor)); valuesOther.push_back(value(*icursor)); } } // for each nnz in row // === Set matrix values. === if (isRowCoarse) { for (unsigned int i = 0; i < cols.size(); i++) cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]); int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor); MatSetValues(getMatCoarseByComponent(rowComponent, colComponent), 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) { for (unsigned int i = 0; i < colsOther.size(); i++) colsOther[i] = interiorMap->getMatIndex(colComponent, colsOther[i]) + rStartInterior; MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 1, &rowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } } else { if ((*interiorMap)[rowComponent].isSet(*cursor) == false) continue; int localRowIndex = (subdomainLevel == 0 ? interiorMap->getLocalMatIndex(rowComponent, *cursor) : interiorMap->getMatIndex(rowComponent, *cursor)); for (unsigned int i = 0; i < cols.size(); i++) { if (subdomainLevel == 0) cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]); else cols[i] = interiorMap->getMatIndex(colComponent, cols[i]); } MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) { int globalRowIndex = interiorMap->getMatIndex(rowComponent, *cursor) + rStartInterior; for (unsigned int i = 0; i < colsOther.size(); i++) colsOther[i] = colCoarseSpace->getMatIndex(colComponent, colsOther[i]); MatSetValues(getMatInteriorCoarseByComponent(colComponent), 1, &globalRowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } } } } } matAssembly(); // === Create solver for the non primal (thus local) variables. === KSPCreate(mpiCommLocal, &kspInterior); KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(kspInterior, "interior_"); KSPSetType(kspInterior, KSPPREONLY); KSPGetPC(kspInterior, &pcInterior); if (isSymmetric) { PCSetType(pcInterior, PCCHOLESKY); PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS); } else { PCSetType(pcInterior, PCLU); if (subdomainLevel == 0) PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK); else PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS); } KSPSetFromOptions(kspInterior); } void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()"); TEST_EXIT_DBG(vec)("No DOF vector defined!\n"); TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n"); // === Transfer values from DOF vector to the PETSc vector. === if (coarseSpaceMap.size()) { for (int i = 0; i < vec->getSize(); i++) setDofVector(getVecRhsInterior(), getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i); } else { for (int i = 0; i < vec->getSize(); i++) setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i); } vecRhsAssembly(); removeDirichletRows(vec, dofMap, getVecRhsInterior()); // === For debugging allow to write the rhs vector to a file. === bool dbgWriteRhs = false; Parameters::get("parallel->debug->write rhs", dbgWriteRhs); if (dbgWriteRhs) { PetscViewer vecView; PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.vec", FILE_MODE_WRITE, &vecView); VecView(getVecRhsInterior(), vecView); PetscViewerDestroy(&vecView); } } void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) { FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()"); int nComponents = vec.getSize(); // === Set old solution to be initiual guess for PETSc solver. === if (!zeroStartVector) { TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n"); VecSet(getVecSolInterior(), 0.0); for (int i = 0; i < nComponents; i++) setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true); vecSolAssembly(); } MatNullSpace matNullspace; Vec nullspaceBasis; if (nullspace.size() > 0 || hasConstantNullspace || constNullspaceComponent.size() > 0) { TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n"); if (constNullspaceComponent.size() > 0) { nullspace.clear(); SystemVector *basisVec = new SystemVector(vec); basisVec->set(0.0); for (unsigned int i = 0; i < constNullspaceComponent.size(); i++) basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0); nullspace.push_back(basisVec); } if (nullspace.size() > 0) { VecDuplicate(getVecSolInterior(), &nullspaceBasis); for (int i = 0; i < nComponents; i++) setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true); VecAssemblyBegin(nullspaceBasis); VecAssemblyEnd(nullspaceBasis); VecNormalize(nullspaceBasis, PETSC_NULL); MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 1, &nullspaceBasis, &matNullspace); MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior()); PetscReal n; VecNorm(getVecSolInterior(), NORM_2, &n); MSG("NORM IS: %e\n", n); } else { MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace); } MSG("NULLSPACE IS NOT REMOVED!\n"); // MatSetNullSpace(getMatInterior(), matNullspace); // KSPSetNullSpace(kspInterior, matNullspace); // === Remove null space, if requested. === if (removeRhsNullspace) { TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n"); MSG("Remove nullspace from rhs vector.\n"); MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL); } } else { TEST_EXIT(removeRhsNullspace == false) ("No nullspace provided that should be removed from rhs!\n"); } // PETSc. solve(getVecRhsInterior(), getVecSolInterior()); if (nullspace.size() > 0) { MatNullSpaceDestroy(&matNullspace); VecDestroy(&nullspaceBasis); } // === Transfere values from PETSc's solution vectors to the DOF vectors. === PetscScalar *vecPointer; VecGetArray(getVecSolInterior(), &vecPointer); int c = 0; for (int component = 0; component < nComponents; component++) { DOFVector &dv = *(vec.getDOFVector(component)); DofMap& d = (*interiorMap)[component].getMap(); for (DofMap::iterator it = d.begin(); it != d.end(); ++it) if (it->second.local != -1) dv[it->first] = vecPointer[c++]; } VecRestoreArray(getVecSolInterior(), &vecPointer); // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to === // === more than one partition. === meshDistributor->synchVector(vec); } void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol) { FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()"); double wtime = MPI::Wtime(); double t0 = 0.0, t1 = 0.0; Vec tmp; if (mpiCommLocal.Get_size() == 1) interiorMap->createLocalVec(tmp); else interiorMap->createVec(tmp); PetscScalar *tmpValues, *rhsValues; VecGetArray(tmp, &tmpValues); VecGetArray(rhs, &rhsValues); for (int i = 0; i < interiorMap->getRankDofs(); i++) tmpValues[i] = rhsValues[i]; VecRestoreArray(rhs, &rhsValues); VecRestoreArray(tmp, &tmpValues); t0 = MPI::Wtime() - wtime; wtime = MPI::Wtime(); KSPSolve(kspInterior, tmp, tmp); t1 = MPI::Wtime() - wtime; wtime = MPI::Wtime(); VecGetArray(tmp, &tmpValues); VecGetArray(sol, &rhsValues); for (int i = 0; i < interiorMap->getRankDofs(); i++) rhsValues[i] = tmpValues[i]; VecRestoreArray(sol, &rhsValues); VecRestoreArray(tmp, &tmpValues); VecDestroy(&tmp); t0 += MPI::Wtime() - wtime; // MSG("TIMEING: %.5f %.5f\n", t0, t1); } void PetscSolverGlobalMatrix::destroyMatrixData() { FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()"); exitPreconditioner(pcInterior); matDestroy(); KSPDestroy(&kspInterior); } void PetscSolverGlobalMatrix::destroyVectorData() { FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()"); vecDestroy(); } void PetscSolverGlobalMatrix::removeDirichletRows(Matrix *seqMat, ParallelDofMapping &dofMap, Mat mpiMat) { FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()"); if (!handleDirichletRows) return; int nComponents = seqMat->getSize(); vector dirichletRows; vector dirichletCols; vector dirichletValues; for (int i = 0; i < nComponents; i++) { bool dirichletRow = false; int dirichletMainCol = -1; for (int j = 0; j < nComponents; j++) { DOFMatrix *dofMat = (*seqMat)[i][j]; if (!dofMat) continue; BoundaryIndexMap& bounds = const_cast(dofMat->getBoundaryManager()->getBoundaryConditionMap()); for (BoundaryIndexMap::iterator bIt = bounds.begin(); bIt != bounds.end(); ++bIt) { if (bIt->second && bIt->second->isDirichlet()) { dirichletRow = true; if ((dynamic_cast(bIt->second))->applyBoundaryCondition()) { TEST_EXIT_DBG(dirichletMainCol == -1)("Should not happen!\n"); dirichletMainCol = j; } } } } if (!dirichletRow) continue; DOFMatrix *dofMat = (*seqMat)[i][dirichletMainCol]; const FiniteElemSpace *feSpace = dofMat->getRowFeSpace(); TEST_EXIT(dofMat->getRowFeSpace() == dofMat->getColFeSpace()) ("I have to think about this scenario! Really possible?\n"); std::set &dRows = dofMat->getDirichletRows(); for (std::set::iterator dofIt = dRows.begin(); dofIt != dRows.end(); ++dofIt) { MultiIndex multiIndex; if (dofMap[feSpace].find(*dofIt, multiIndex) && dofMap[feSpace].isRankDof(*dofIt)) { int rowIndex = dofMap.getMatIndex(i, multiIndex.global); int colIndex = dofMap.getMatIndex(dirichletMainCol, multiIndex.global); dirichletRows.push_back(rowIndex); dirichletCols.push_back(colIndex); dirichletValues.push_back(1.0); } } } MatZeroRows(mpiMat, dirichletRows.size(), &(dirichletRows[0]), 0.0, PETSC_NULL, PETSC_NULL); for (int i = 0; i < static_cast(dirichletRows.size()); i++) MatSetValue(mpiMat, dirichletRows[i], dirichletCols[i], dirichletValues[i], INSERT_VALUES); MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY); } void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec, ParallelDofMapping &dofMap, Vec mpiVec) { FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()"); if (!handleDirichletRows) return; int nComponents = seqVec->getSize(); for (int i = 0; i < nComponents; i++) { DOFVector *dofVec = seqVec->getDOFVector(i); const FiniteElemSpace *feSpace = dofVec->getFeSpace(); map& dValues = dofVec->getDirichletValues(); for (map::iterator dofIt = dValues.begin(); dofIt != dValues.end(); ++dofIt) { MultiIndex multiIndex; if (dofMap[feSpace].find(dofIt->first, multiIndex) && dofMap[feSpace].isRankDof(dofIt->first)) VecSetValue(mpiVec, dofMap.getMatIndex(i, multiIndex.global), dofIt->second, INSERT_VALUES); } } VecAssemblyBegin(mpiVec); VecAssemblyEnd(mpiVec); } void PetscSolverGlobalMatrix::createFieldSplit(PC pc) { FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()"); vector isNames; Parameters::get("parallel->solver->is blocks", isNames); int nBlocks = isNames.size(); if (nBlocks == 0) return; for (int i = 0; i < nBlocks; i++) { MSG("Create for block %s\n", isNames[i].c_str()); vector blockComponents; Parameters::get("parallel->solver->is block " + lexical_cast(i), blockComponents); int nComponents = static_cast(blockComponents.size()); TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i); // Check if blocks are continous for (int j = 0; j < nComponents; j++) { TEST_EXIT(blockComponents[j] == blockComponents[0] + j) ("Does not yet support not continous IS blocks! Block %s\n", isNames[i].c_str()); } IS is; interiorMap->createIndexSet(is, blockComponents[0], nComponents); PCFieldSplitSetIS(pc, isNames[i].c_str(), is); ISDestroy(&is); } } void PetscSolverGlobalMatrix::initPreconditioner(PC pc) { FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()"); PCSetFromOptions(pc); createFieldSplit(pc); } void PetscSolverGlobalMatrix::exitPreconditioner(PC pc) { FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()"); } void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat, int nRowMat, int nColMat) { FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()"); TEST_EXIT(seqMat)("No DOFMatrix!\n"); using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; traits::col::type col(seqMat->getBaseMatrix()); traits::const_value::type value(seqMat->getBaseMatrix()); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; vector cols; vector values; cols.reserve(300); values.reserve(300); vector globalCols; // Get periodic mapping object PeriodicMap &perMap = meshDistributor->getPeriodicMap(); const FiniteElemSpace *rowFe = seqMat->getRowFeSpace(); const FiniteElemSpace *colFe = seqMat->getColFeSpace(); // === Traverse all rows of the DOF matrix and insert row wise the values === // === to the PETSc matrix. === for (cursor_type cursor = begin(seqMat->getBaseMatrix()), cend = end(seqMat->getBaseMatrix()); cursor != cend; ++cursor) { // Global index of the current row DOF. MultiIndex rowMultiIndex; if ((*interiorMap)[nRowMat].find(*cursor, rowMultiIndex) == false) continue; int globalRowDof = rowMultiIndex.global; // Test if the current row DOF is a periodic DOF. bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof); if (!periodicRow) { // === Row DOF index is not periodic. === // Get PETSc's mat row index. int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof); cols.clear(); values.clear(); for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // Global index of the current column index. MultiIndex colMultiIndex; if ((*interiorMap)[nColMat].find(col(*icursor), colMultiIndex) == false) continue; int globalColDof = colMultiIndex.global; // Test if the current col dof is a periodic dof. bool periodicCol = perMap.isPeriodic(colFe, globalColDof); // Get PETSc's mat col index. int colIndex = interiorMap->getMatIndex(nColMat, globalColDof); // Ignore all zero entries, expect it is a diagonal entry. if (value(*icursor) == 0.0 && rowIndex != colIndex) continue; if (!periodicCol) { // Calculate the exact position of the column index in the PETSc matrix. cols.push_back(colIndex); values.push_back(value(*icursor)); } else { // === Row index is not periodic, but column index is. === // Create set of all periodic associations of the column index. std::set perAsc; perMap.fillAssociations(colFe, globalColDof, meshDistributor->getElementObjectDb(), perAsc); // Scale value to the number of periodic associations of the column index. double scaledValue = value(*icursor) * pow(0.5, static_cast(perAsc.size())); // === Create set of all matrix column indices due to the periodic === // === associations of the column DOF index. === vector newCols; perMap.mapDof(colFe, globalColDof, perAsc, newCols); for (unsigned int i = 0; i < newCols.size(); i++) { cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i])); values.push_back(scaledValue); } } } MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); } else { // === Row DOF index is periodic. === // Because this row is periodic, we will have to add the entries of this // matrix row to multiple rows. The following maps store to each row an // array of column indices and values of the entries that must be added to // the PETSc matrix. map > colsMap; map > valsMap; // Traverse all column entries. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // Global index of the current column index. int globalColDof = (*interiorMap)[nColMat][col(*icursor)].global; // Ignore all zero entries, expect it is a diagonal entry. if (value(*icursor) == 0.0 && globalRowDof != globalColDof) continue; // === Add all periodic associations of both, the row and the column === // === indices to the set perAsc. === std::set perAsc; perMap.fillAssociations(colFe, globalColDof, meshDistributor->getElementObjectDb(), perAsc); perMap.fillAssociations(rowFe, globalRowDof, meshDistributor->getElementObjectDb(), perAsc); // Scale the value with respect to the number of periodic associations. double scaledValue = value(*icursor) * pow(0.5, static_cast(perAsc.size())); // === Create all matrix entries with respect to the periodic === // === associations of the row and column indices. === vector > entry; perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof), perAsc, entry); // === Translate the matrix entries to PETSc's matrix. for (unsigned int i = 0; i < entry.size(); i++) { int rowIdx = interiorMap->getMatIndex(nRowMat, entry[i].first); int colIdx = interiorMap->getMatIndex(nColMat, entry[i].second); colsMap[rowIdx].push_back(colIdx); valsMap[rowIdx].push_back(scaledValue); } } // === Finally, add all periodic rows to the PETSc matrix. === for (map >::iterator rowIt = colsMap.begin(); rowIt != colsMap.end(); ++rowIt) { TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size()) ("Should not happen!\n"); int rowIndex = rowIt->first; MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(), &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES); } } } } void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, Vec vecCoarse, DOFVector* vec, int nRowVec, bool rankOnly) { FUNCNAME("PetscSolverGlobalMatrix::setDofVector()"); const FiniteElemSpace *feSpace = vec->getFeSpace(); PeriodicMap &perMap = meshDistributor->getPeriodicMap(); ParallelDofMapping *rowCoarseSpace = (coarseSpaceMap.size() ? coarseSpaceMap[nRowVec] : NULL); // Traverse all used DOFs in the dof vector. DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { if (rankOnly && !(*interiorMap)[nRowVec].isRankDof(dofIt.getDOFIndex())) continue; if (isCoarseSpace(nRowVec, dofIt.getDOFIndex())) { TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n"); int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex()); VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES); } else { if ((*interiorMap)[nRowVec].isSet(dofIt.getDOFIndex()) == false) continue; // Calculate global row index of the DOF. DegreeOfFreedom globalRowDof = (*interiorMap)[nRowVec][dofIt.getDOFIndex()].global; // Get PETSc's mat index of the row DOF. int index = 0; if (interiorMap->isMatIndexFromGlobal()) index = interiorMap->getMatIndex(nRowVec, globalRowDof) + rStartInterior; else index = interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior; if (perMap.isPeriodic(feSpace, globalRowDof)) { std::set& perAsc = perMap.getAssociations(feSpace, globalRowDof); double value = *dofIt / (perAsc.size() + 1.0); VecSetValue(vecInterior, index, value, ADD_VALUES); for (std::set::iterator perIt = perAsc.begin(); perIt != perAsc.end(); ++perIt) { int mappedDof = perMap.map(feSpace, *perIt, globalRowDof); int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof); VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES); } } else { // The DOF index is not periodic. VecSetValue(vecInterior, index, *dofIt, ADD_VALUES); } } } } }