// // 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 *mat) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()"); if (coarseSpaceMap != NULL) { fillPetscMatrixWithCoarseSpace(mat); return; } TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n"); TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n"); TEST_EXIT_DBG(mat)("No DOF matrix defined!\n"); double wtime = MPI::Wtime(); // === Check if mesh was changed and, in this case, recompute matrix === // === nnz structure and matrix indices. === int recvAllValues = 0; int sendValue = static_cast(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz); mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM); if (!d_nnz || recvAllValues != 0 || alwaysCreateNnzStructure) { vector feSpaces = getFeSpaces(mat); interiorMap->setComputeMatIndex(true, true); interiorMap->update(feSpaces); if (d_nnz) { delete [] d_nnz; d_nnz = NULL; delete [] o_nnz; o_nnz = NULL; } createPetscNnzStructure(mat); lastMeshNnz = meshDistributor->getLastMeshChangeIndex(); } // === Create PETSc vector (solution and a temporary vector). === int nRankRows = interiorMap->getRankDofs(); int nOverallRows = interiorMap->getOverallDofs(); VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &petscSolVec); int testddd = 1; Parameters::get("block size", testddd); if (testddd > 1) VecSetBlockSize(petscSolVec, testddd); // === Create PETSc matrix with the computed nnz data structure. === MatCreateMPIAIJ(mpiCommGlobal, nRankRows, nRankRows, nOverallRows, nOverallRows, 0, d_nnz, 0, o_nnz, &matIntInt); if (testddd > 1) { MatSetBlockSize(matIntInt, testddd); MSG("MAT SET BLOCK SIZE: %d\n", testddd); } #if (DEBUG != 0) MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif #if (DEBUG != 0) int a, b; MatGetOwnershipRange(matIntInt, &a, &b); TEST_EXIT(a == interiorMap->getStartDofs())("Wrong matrix ownership range!\n"); TEST_EXIT(b == interiorMap->getStartDofs() + nRankRows) ("Wrong matrix ownership range!\n"); #endif // === Transfer values from DOF matrices to the PETSc matrix. === int nComponents = mat->getNumRows(); for (int i = 0; i < nComponents; i++) for (int j = 0; j < nComponents; j++) if ((*mat)[i][j]) setDofMatrix((*mat)[i][j], i, j); #if (DEBUG != 0) MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime); #endif MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY); // === Init PETSc solver. === KSPCreate(mpiCommGlobal, &kspInterior); KSPGetPC(kspInterior, &pcInterior); KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN); KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT); KSPSetType(kspInterior, KSPBCGS); KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str()); KSPSetFromOptions(kspInterior); PCSetFromOptions(pcInterior); // Do not delete the solution vector, use it for the initial guess. if (!zeroStartVector) KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE); MSG("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime); } void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix *mat) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()"); vector feSpaces = getFeSpaces(mat); int nRowsRankInterior = interiorMap->getRankDofs(); int nRowsOverallInterior = interiorMap->getOverallDofs(); if (subdomainLevel == 0) { nGlobalOverallInterior = nRowsOverallInterior; MatCreateSeqAIJ(mpiCommLocal, nRowsRankInterior, nRowsRankInterior, 60, PETSC_NULL, &matIntInt); } else { MatCreateMPIAIJ(mpiCommLocal, nRowsRankInterior, nRowsRankInterior, nRowsOverallInterior, nRowsOverallInterior, 60, PETSC_NULL, 60, PETSC_NULL, &matIntInt); } if (coarseSpaceMap) { int nRowsRankCoarse = coarseSpaceMap->getRankDofs(); int nRowsOverallCoarse = coarseSpaceMap->getOverallDofs(); MatCreateMPIAIJ(mpiCommGlobal, nRowsRankCoarse, nRowsRankCoarse, nRowsOverallCoarse, nRowsOverallCoarse, 60, PETSC_NULL, 60, PETSC_NULL, &matCoarseCoarse); MatCreateMPIAIJ(mpiCommGlobal, nRowsRankCoarse, nRowsRankInterior, nRowsOverallCoarse, nGlobalOverallInterior, 60, PETSC_NULL, 60, PETSC_NULL, &matCoarseInt); MatCreateMPIAIJ(mpiCommGlobal, nRowsRankInterior, nRowsRankCoarse, nGlobalOverallInterior, nRowsOverallCoarse, 60, PETSC_NULL, 60, PETSC_NULL, &matIntCoarse); } // === 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 = mat->getSize(); for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) { if (!(*mat)[i][j]) continue; traits::col::type col((*mat)[i][j]->getBaseMatrix()); traits::const_value::type value((*mat)[i][j]->getBaseMatrix()); // Traverse all rows. for (cursor_type cursor = begin((*mat)[i][j]->getBaseMatrix()), cend = end((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) { bool rowPrimal = isCoarseSpace(feSpaces[i], *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 colPrimal = isCoarseSpace(feSpaces[j], col(*icursor)); if (colPrimal) { if (rowPrimal) { cols.push_back(col(*icursor)); values.push_back(value(*icursor)); } else { colsOther.push_back(col(*icursor)); valuesOther.push_back(value(*icursor)); } } else { if (rowPrimal) { colsOther.push_back(col(*icursor)); valuesOther.push_back(value(*icursor)); } else { cols.push_back(col(*icursor)); values.push_back(value(*icursor)); } } } // for each nnz in row // === Set matrix values. === if (rowPrimal) { int rowIndex = coarseSpaceMap->getMatIndex(i, *cursor); for (unsigned int k = 0; k < cols.size(); k++) cols[k] = coarseSpaceMap->getMatIndex(j, cols[k]); MatSetValues(matCoarseCoarse, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) { if (subdomainLevel == 0) { for (unsigned int k = 0; k < colsOther.size(); k++) colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]); } else { for (unsigned int k = 0; k < colsOther.size(); k++) colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]) + rStartInterior; } MatSetValues(matCoarseInt, 1, &rowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } } else { int localRowIndex = (subdomainLevel == 0 ? interiorMap->getLocalMatIndex(i, *cursor) : interiorMap->getMatIndex(i, *cursor)); for (unsigned int k = 0; k < cols.size(); k++) { if (subdomainLevel == 0) cols[k] = interiorMap->getLocalMatIndex(j, cols[k]); else cols[k] = interiorMap->getMatIndex(j, cols[k]); } MatSetValues(matIntInt, 1, &localRowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); if (colsOther.size()) { int globalRowIndex = interiorMap->getMatIndex(i, *cursor); if (subdomainLevel != 0) globalRowIndex += rStartInterior; for (unsigned int k = 0; k < colsOther.size(); k++) colsOther[k] = coarseSpaceMap->getMatIndex(j, colsOther[k]); MatSetValues(matIntCoarse, 1, &globalRowIndex, colsOther.size(), &(colsOther[0]), &(valuesOther[0]), ADD_VALUES); } } } } } // === Start global assembly procedure. === MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY); if (coarseSpaceMap) { MatAssemblyBegin(matCoarseCoarse, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matCoarseCoarse, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(matIntCoarse, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matIntCoarse, MAT_FINAL_ASSEMBLY); MatAssemblyBegin(matCoarseInt, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(matCoarseInt, MAT_FINAL_ASSEMBLY); } // === Create solver for the non primal (thus local) variables. === KSPCreate(mpiCommLocal, &kspInterior); KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN); KSPSetOptionsPrefix(kspInterior, "interior_"); KSPSetType(kspInterior, KSPPREONLY); PC pcInterior; KSPGetPC(kspInterior, &pcInterior); PCSetType(pcInterior, PCLU); if (subdomainLevel == 0) PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK); else PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS); KSPSetFromOptions(kspInterior); } void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec) { FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()"); if (coarseSpaceMap) { fillPetscRhsWithCoarseSpace(vec); return; } TEST_EXIT_DBG(vec)("No DOF vector defined!\n"); TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n"); int nRankRows = interiorMap->getRankDofs(); int nOverallRows = interiorMap->getOverallDofs(); VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &rhsInterior); int testddd = 1; Parameters::get("block size", testddd); if (testddd > 1) VecSetBlockSize(rhsInterior, testddd); // === Transfer values from DOF vector to the PETSc vector. === for (int i = 0; i < vec->getSize(); i++) setDofVector(rhsInterior, vec->getDOFVector(i), i); VecAssemblyBegin(rhsInterior); VecAssemblyEnd(rhsInterior); if (removeRhsNullSpace) { MSG("Remove constant null space from the RHS!\n"); MatNullSpace sp; MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &sp); MatNullSpaceRemove(sp, rhsInterior, PETSC_NULL); MatNullSpaceDestroy(&sp); } } void PetscSolverGlobalMatrix::fillPetscRhsWithCoarseSpace(SystemVector *vec) { FUNCNAME("SubDomainSolver::fillPetscRhs()"); VecCreateMPI(mpiCommGlobal, interiorMap->getRankDofs(), nGlobalOverallInterior, &rhsInterior); if (coarseSpaceMap) VecCreateMPI(mpiCommGlobal, coarseSpaceMap->getRankDofs(), coarseSpaceMap->getOverallDofs(), &rhsCoarseSpace); for (int i = 0; i < vec->getSize(); i++) { const FiniteElemSpace *feSpace = vec->getDOFVector(i)->getFeSpace(); DOFVector::Iterator dofIt(vec->getDOFVector(i), USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { int index = dofIt.getDOFIndex(); if (isCoarseSpace(feSpace, index)) { index = coarseSpaceMap->getMatIndex(i, index); VecSetValue(rhsCoarseSpace, index, *dofIt, ADD_VALUES); } else { index = interiorMap->getMatIndex(i, index) + rStartInterior; VecSetValue(rhsInterior, index, *dofIt, ADD_VALUES); } } } VecAssemblyBegin(rhsInterior); VecAssemblyEnd(rhsInterior); if (coarseSpaceMap) { VecAssemblyBegin(rhsCoarseSpace); VecAssemblyEnd(rhsCoarseSpace); } } 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) { VecSet(petscSolVec, 0.0); for (int i = 0; i < nComponents; i++) setDofVector(petscSolVec, vec.getDOFVector(i), i, true); VecAssemblyBegin(petscSolVec); VecAssemblyEnd(petscSolVec); } // PETSc. solve(rhsInterior, petscSolVec); // === Transfere values from PETSc's solution vectors to the DOF vectors. === PetscScalar *vecPointer; VecGetArray(petscSolVec, &vecPointer); int c = 0; for (int i = 0; i < nComponents; i++) { DOFVector &dv = *(vec.getDOFVector(i)); DofMap& d = (*interiorMap)[dv.getFeSpace()].getMap(); for (DofMap::iterator it = d.begin(); it != d.end(); ++it) if (it->second.local != -1) dv[it->first] = vecPointer[c++]; } VecRestoreArray(petscSolVec, &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()"); Vec tmp; if (mpiCommLocal.Get_size() == 1) VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp); else VecCreateMPI(mpiCommLocal, interiorMap->getRankDofs(), interiorMap->getOverallDofs(), &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); KSPSolve(kspInterior, tmp, tmp); 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); } void PetscSolverGlobalMatrix::destroyMatrixData() { FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()"); MatDestroy(&matIntInt); KSPDestroy(&kspInterior); VecDestroy(&petscSolVec); } void PetscSolverGlobalMatrix::destroyVectorData() { FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()"); VecDestroy(&rhsInterior); } void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* mat, int nRowMat, int nColMat) { FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()"); TEST_EXIT(mat)("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(mat->getBaseMatrix()); traits::const_value::type value(mat->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(); // === Traverse all rows of the dof matrix and insert row wise the values === // === to the PETSc matrix. === for (cursor_type cursor = begin(mat->getBaseMatrix()), cend = end(mat->getBaseMatrix()); cursor != cend; ++cursor) { const FiniteElemSpace *rowFe = mat->getRowFeSpace(); const FiniteElemSpace *colFe = mat->getColFeSpace(); // Global index of the current row DOF. int globalRowDof = (*interiorMap)[rowFe][*cursor].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. int globalColDof = (*interiorMap)[colFe][col(*icursor)].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; std::set& perColAsc = perMap.getAssociations(colFe, globalColDof); for (std::set::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it) if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it)) perAsc.insert(*it); // 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; // First, add the original matrix index. newCols.push_back(globalColDof); // And add all periodic matrix indices. for (std::set::iterator it = perAsc.begin(); it != perAsc.end(); ++it) { int nCols = static_cast(newCols.size()); for (int i = 0; i < nCols; i++) { TEST_EXIT_DBG(perMap.isPeriodic(colFe, *it, newCols[i])) ("Wrong periodic DOF associations at boundary %d with DOF %d!\n", *it, newCols[i]); newCols.push_back(perMap.map(colFe, *it, newCols[i])); } } for (unsigned int i = 0; i < newCols.size(); i++) { cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i])); values.push_back(scaledValue); } } } MatSetValues(matIntInt, 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)[colFe][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; if (perMap.isPeriodic(colFe, globalColDof)) { std::set& perColAsc = perMap.getAssociations(colFe, globalColDof); for (std::set::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it) if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it)) perAsc.insert(*it); } std::set& perRowAsc = perMap.getAssociations(rowFe, globalRowDof); for (std::set::iterator it = perRowAsc.begin(); it != perRowAsc.end(); ++it) if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it)) perAsc.insert(*it); // 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; // First, add the original entry. entry.push_back(make_pair(globalRowDof, globalColDof)); // Then, traverse the periodic associations of the row and column // indices and create the corresponding entries. for (std::set::iterator it = perAsc.begin(); it != perAsc.end(); ++it) { int nEntry = static_cast(entry.size()); for (int i = 0; i < nEntry; i++) { int perRowDof = 0; if (perMap.isPeriodic(rowFe, *it, entry[i].first)) perRowDof = perMap.map(rowFe, *it, entry[i].first); else perRowDof = entry[i].first; int perColDof; if (perMap.isPeriodic(colFe, *it, entry[i].second)) perColDof = perMap.map(colFe, *it, entry[i].second); else perColDof = entry[i].second; entry.push_back(make_pair(perRowDof, perColDof)); } } // === 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(matIntInt, 1, &rowIndex, rowIt->second.size(), &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES); } } } } void PetscSolverGlobalMatrix::setDofVector(Vec& petscVec, DOFVector* vec, int nRowVec, bool rankOnly) { FUNCNAME("PetscSolverGlobalMatrix::setDofVector()"); const FiniteElemSpace *feSpace = vec->getFeSpace(); PeriodicMap &perMap = meshDistributor->getPeriodicMap(); // Traverse all used DOFs in the dof vector. DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex())) continue; // Calculate global row index of the DOF. DegreeOfFreedom globalRowDof = (*interiorMap)[feSpace][dofIt.getDOFIndex()].global; // Get PETSc's mat index of the row DOF. int index = interiorMap->getMatIndex(nRowVec, globalRowDof); if (perMap.isPeriodic(feSpace, globalRowDof)) { std::set& perAsc = perMap.getAssociations(feSpace, globalRowDof); double value = *dofIt / (perAsc.size() + 1.0); VecSetValues(petscVec, 1, &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); VecSetValues(petscVec, 1, &mappedIndex, &value, ADD_VALUES); } } else { // The DOF index is not periodic. double value = *dofIt; VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); } } } void PetscSolverGlobalMatrix::createPetscNnzStructure(Matrix *mat) { FUNCNAME("PetscSolverGlobalMatrix::createPetscNnzStructure()"); TEST_EXIT_DBG(!d_nnz)("There is something wrong!\n"); TEST_EXIT_DBG(!o_nnz)("There is something wrong!\n"); vector feSpaces = getFeSpaces(mat); int nRankRows = interiorMap->getRankDofs(); int rankStartIndex = interiorMap->getStartDofs(); d_nnz = new int[nRankRows]; o_nnz = new int[nRankRows]; for (int i = 0; i < nRankRows; i++) { d_nnz[i] = 0; o_nnz[i] = 0; } 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 vector > MatrixNnzEntry; typedef map RankToDofContainer; // Stores to each rank a list of nnz entries (i.e. pairs of row and column // index) that this rank will send to. These nnz entries will be assembled // on this rank, but because the row DOFs are not DOFs of this rank they // will be send to the owner of the row DOFs. map sendMatrixEntry; // Maps to each DOF that must be send to another rank the rank number of the // receiving rank. map, int> sendDofToRank; // First, create for all ranks, to which we send data to, MatrixNnzEntry // object with 0 entries. for (unsigned int i = 0; i < feSpaces.size(); i++) { for (DofComm::Iterator it(meshDistributor->getDofComm().getRecvDofs(), feSpaces[i]); !it.end(); it.nextRank()) { sendMatrixEntry[it.getRank()].resize(0); for (; !it.endDofIter(); it.nextDof()) sendDofToRank[make_pair(it.getDofIndex(), i)] = it.getRank(); } } // Create list of ranks from which we receive data from. std::set recvFromRank; for (unsigned int i = 0; i < feSpaces.size(); i++) for (DofComm::Iterator it(meshDistributor->getDofComm().getSendDofs(), feSpaces[i]); !it.end(); it.nextRank()) recvFromRank.insert(it.getRank()); // === Traverse matrices to create nnz data. === int nComponents = mat->getNumRows(); for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) { if (!(*mat)[i][j]) continue; TEST_EXIT_DBG((*mat)[i][j]->getRowFeSpace() == feSpaces[i]) ("Should not happen!\n"); TEST_EXIT_DBG((*mat)[i][j]->getColFeSpace() == feSpaces[j]) ("Should not happen!\n"); Matrix bmat = (*mat)[i][j]->getBaseMatrix(); traits::col::type col(bmat); traits::const_value::type value(bmat); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; for (cursor_type cursor = begin(bmat), cend = end(bmat); cursor != cend; ++cursor) { int globalRowDof = (*interiorMap)[feSpaces[i]][*cursor].global; // The corresponding global matrix row index of the current row DOF. int petscRowIdx = interiorMap->getMatIndex(i, globalRowDof); if ((*interiorMap)[feSpaces[i]].isRankDof(*cursor)) { // === The current row DOF is a rank DOF, so create the === // === corresponding nnz values directly on rank's nnz data. === // This is the local row index of the local PETSc matrix. int localPetscRowIdx = petscRowIdx - rankStartIndex; TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows) ("Should not happen! \n Debug info: localRowIdx = %d globalRowIndx = %d petscRowIdx = %d localPetscRowIdx = %d rStart = %d nCompontens = %d nRankRows = %d\n", *cursor, (*interiorMap)[feSpaces[i]][*cursor].global, petscRowIdx, localPetscRowIdx, rankStartIndex, nComponents, nRankRows); // Traverse all non zero entries in this row. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { int globalColDof = (*interiorMap)[feSpaces[j]][col(*icursor)].global; int petscColIdx = interiorMap->getMatIndex(j, globalColDof); if (value(*icursor) != 0.0 || petscRowIdx == petscColIdx) { // The row DOF is a rank DOF, if also the column is a rank DOF, // increment the d_nnz values for this row, otherwise the // o_nnz value. if (petscColIdx >= rankStartIndex && petscColIdx < rankStartIndex + nRankRows) d_nnz[localPetscRowIdx]++; else o_nnz[localPetscRowIdx]++; } } } else { // === The current row DOF is not a rank DOF, i.e., its values === // === are also created on this rank, but afterthere they will === // === be send to another rank. So we need to send also the === // === corresponding nnz structure of this row to the corres- === // === ponding rank. === // Send all non zero entries to the member of the row DOF. int sendToRank = sendDofToRank[make_pair(*cursor, i)]; for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { if (value(*icursor) != 0.0) { int globalColDof = (*interiorMap)[feSpaces[j]][col(*icursor)].global; int petscColIdx = interiorMap->getMatIndex(j, globalColDof); sendMatrixEntry[sendToRank]. push_back(make_pair(petscRowIdx, petscColIdx)); } } } // if (isRankDof[*cursor]) ... else ... } // for each row in mat[i][j] } } // === Send and recv the nnz row structure to/from other ranks. === StdMpi stdMpi(mpiCommGlobal, true); stdMpi.send(sendMatrixEntry); for (std::set::iterator it = recvFromRank.begin(); it != recvFromRank.end(); ++it) stdMpi.recv(*it); stdMpi.startCommunication(); // === Evaluate the nnz structure this rank got from other ranks and add === // === it to the PETSc nnz data structure. === for (map::iterator it = stdMpi.getRecvData().begin(); it != stdMpi.getRecvData().end(); ++it) { if (it->second.size() > 0) { for (unsigned int i = 0; i < it->second.size(); i++) { int r = it->second[i].first; int c = it->second[i].second; int localRowIdx = r - rankStartIndex; TEST_EXIT_DBG(localRowIdx >= 0 && localRowIdx < nRankRows) ("Got row index %d/%d (nRankRows = %d) from rank %d. Should not happen!\n", r, localRowIdx, nRankRows, it->first); if (c < rankStartIndex || c >= rankStartIndex + nRankRows) o_nnz[localRowIdx]++; else d_nnz[localRowIdx]++; } } } // The above algorithm for calculating the number of nnz per row over- // approximates the value, i.e., the number is always equal or larger to // the real number of nnz values in the global parallel matrix. For small // matrices, the problem may arise, that the result is larger than the // number of elements in a row. This is fixed in the following. if (nRankRows < 100) for (int i = 0; i < nRankRows; i++) d_nnz[i] = std::min(d_nnz[i], nRankRows); } }