diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp index 0739f70371c71badd506efcccc968b71c047af6a..fa3848ed31d4e92deb674d415141005f2c8fda46 100644 --- a/Solver/eigenSolver.cpp +++ b/Solver/eigenSolver.cpp @@ -15,7 +15,7 @@ #define EPSDestroy(e) EPSDestroy(*(e)) #endif -void eigenSolver::_try(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); } +void eigenSolver::_check(int ierr) const { CHKERRABORT(PETSC_COMM_WORLD, ierr); } eigenSolver::eigenSolver(dofManager<double> *manager, std::string A, std::string B, bool hermitian) @@ -45,52 +45,52 @@ bool eigenSolver::solve(int numEigenValues, std::string which, Mat B = _B ? _B->getMatrix() : PETSC_NULL; PetscInt N, M; - _try(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); - _try(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); - _try(MatGetSize(A, &N, &M)); + _check(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); + _check(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); + _check(MatGetSize(A, &N, &M)); PetscInt N2, M2; if(_B) { - _try(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); - _try(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); - _try(MatGetSize(B, &N2, &M2)); + _check(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); + _check(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); + _check(MatGetSize(B, &N2, &M2)); } // generalized eigenvalue problem A x - \lambda B x = 0 EPS eps; - _try(EPSCreate(PETSC_COMM_WORLD, &eps)); - _try(EPSSetOperators(eps, A, B)); + _check(EPSCreate(PETSC_COMM_WORLD, &eps)); + _check(EPSSetOperators(eps, A, B)); if(_hermitian) - _try(EPSSetProblemType(eps, _B ? EPS_GHEP : EPS_HEP)); + _check(EPSSetProblemType(eps, _B ? EPS_GHEP : EPS_HEP)); else - _try(EPSSetProblemType(eps, _B ? EPS_GNHEP : EPS_NHEP)); + _check(EPSSetProblemType(eps, _B ? EPS_GNHEP : EPS_NHEP)); // set some default options - _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); - _try(EPSSetTolerances(eps, tolVal, iterMax)); + _check(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); + _check(EPSSetTolerances(eps, tolVal, iterMax)); if(method == "krylovschur") - _try(EPSSetType(eps, EPSKRYLOVSCHUR)); + _check(EPSSetType(eps, EPSKRYLOVSCHUR)); else if(method == "arnoldi") - _try(EPSSetType(eps, EPSARNOLDI)); + _check(EPSSetType(eps, EPSARNOLDI)); else if(method == "arpack") - _try(EPSSetType(eps, EPSARPACK)); + _check(EPSSetType(eps, EPSARPACK)); else if(method == "power") - _try(EPSSetType(eps, EPSPOWER)); + _check(EPSSetType(eps, EPSPOWER)); else Msg::Fatal("eigenSolver: method '%s' not available", method.c_str()); // override these options at runtime, petsc-style - _try(EPSSetFromOptions(eps)); + _check(EPSSetFromOptions(eps)); // force options specified directly as arguments if(numEigenValues) - _try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); + _check(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE)); if(which == "smallest") - _try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE)); + _check(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE)); else if(which == "smallestReal") - _try(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL)); + _check(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL)); else if(which == "largest") - _try(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE)); + _check(EPSSetWhichEigenpairs(eps, EPS_LARGEST_MAGNITUDE)); // print info #if(SLEPC_VERSION_RELEASE == 0 || \ @@ -100,27 +100,27 @@ bool eigenSolver::solve(int numEigenValues, std::string which, #else const EPSType type; #endif - _try(EPSGetType(eps, &type)); + _check(EPSGetType(eps, &type)); Msg::Debug("SLEPc solution method: %s", type); PetscInt nev; - _try(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL)); + _check(EPSGetDimensions(eps, &nev, PETSC_NULL, PETSC_NULL)); Msg::Debug("SLEPc number of requested eigenvalues: %d", nev); PetscReal tol; PetscInt maxit; - _try(EPSGetTolerances(eps, &tol, &maxit)); + _check(EPSGetTolerances(eps, &tol, &maxit)); Msg::Debug("SLEPc stopping condition: tol=%g, maxit=%d", tol, maxit); // solve Msg::Info("SLEPc solving..."); double t1 = Cpu(); - _try(EPSSolve(eps)); + _check(EPSSolve(eps)); // check convergence int its; - _try(EPSGetIterationNumber(eps, &its)); + _check(EPSGetIterationNumber(eps, &its)); EPSConvergedReason reason; - _try(EPSGetConvergedReason(eps, &reason)); + _check(EPSGetConvergedReason(eps, &reason)); if(reason == EPS_CONVERGED_TOL) { double t2 = Cpu(); Msg::Debug("SLEPc converged in %d iterations (%g s)", its, t2 - t1); @@ -137,7 +137,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which, // get number of converged approximate eigenpairs PetscInt nconv; - _try(EPSGetConverged(eps, &nconv)); + _check(EPSGetConverged(eps, &nconv)); Msg::Debug("SLEPc number of converged eigenpairs: %d", nconv); // ignore additional eigenvalues if we get more than what we asked @@ -146,22 +146,22 @@ bool eigenSolver::solve(int numEigenValues, std::string which, if(nconv > 0) { Vec xr, xi; #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) - _try(MatGetVecs(A, PETSC_NULL, &xr)); - _try(MatGetVecs(A, PETSC_NULL, &xi)); + _check(MatGetVecs(A, PETSC_NULL, &xr)); + _check(MatGetVecs(A, PETSC_NULL, &xi)); #else - _try(MatCreateVecs(A, PETSC_NULL, &xr)); - _try(MatCreateVecs(A, PETSC_NULL, &xi)); + _check(MatCreateVecs(A, PETSC_NULL, &xr)); + _check(MatCreateVecs(A, PETSC_NULL, &xi)); #endif Msg::Debug(" Re[EigenValue] Im[EigenValue]" " Relative error"); for(int i = 0; i < nconv; i++) { PetscScalar kr, ki; - _try(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi)); + _check(EPSGetEigenpair(eps, i, &kr, &ki, xr, xi)); PetscReal error; #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 6) - _try(EPSComputeRelativeError(eps, i, &error)); + _check(EPSComputeRelativeError(eps, i, &error)); #else - _try(EPSComputeError(eps, i, EPS_ERROR_RELATIVE, &error)); + _check(EPSComputeError(eps, i, EPS_ERROR_RELATIVE, &error)); #endif #if defined(PETSC_USE_COMPLEX) PetscReal re = PetscRealPart(kr); @@ -176,8 +176,8 @@ bool eigenSolver::solve(int numEigenValues, std::string which, // store eigenvalues and eigenvectors _eigenValues.push_back(std::complex<double>(re, im)); PetscScalar *tmpr, *tmpi; - _try(VecGetArray(xr, &tmpr)); - _try(VecGetArray(xi, &tmpi)); + _check(VecGetArray(xr, &tmpr)); + _check(VecGetArray(xi, &tmpi)); std::vector<std::complex<double> > ev(N); for(int i = 0; i < N; i++) { #if defined(PETSC_USE_COMPLEX) @@ -188,11 +188,11 @@ bool eigenSolver::solve(int numEigenValues, std::string which, } _eigenVectors.push_back(ev); } - _try(VecDestroy(&xr)); - _try(VecDestroy(&xi)); + _check(VecDestroy(&xr)); + _check(VecDestroy(&xi)); } - _try(EPSDestroy(&eps)); + _check(EPSDestroy(&eps)); if(reason == EPS_CONVERGED_TOL) { Msg::Debug("SLEPc done"); diff --git a/Solver/eigenSolver.h b/Solver/eigenSolver.h index 9a74c79f2c4239fa87185589da114ad56342377b..c722a31ec3baa1df1a89ea8929dc2bbc65b1b5b1 100644 --- a/Solver/eigenSolver.h +++ b/Solver/eigenSolver.h @@ -22,7 +22,7 @@ private: bool _hermitian; std::vector<std::complex<double> > _eigenValues; std::vector<std::vector<std::complex<double> > > _eigenVectors; - void _try(int ierr) const; + void _check(int ierr) const; public: eigenSolver(dofManager<double> *manager, std::string A, std::string B = "", diff --git a/Solver/linearSystemPETSc.cpp b/Solver/linearSystemPETSc.cpp index a9891f2bf8fa0830d656253d1a157322c1569d8b..197d9d17945f3c2af7c77b27fcdf2911458a913a 100644 --- a/Solver/linearSystemPETSc.cpp +++ b/Solver/linearSystemPETSc.cpp @@ -25,9 +25,9 @@ template <> void linearSystemPETSc<double>::getFromRightHandSide(int row, double &val) const { PetscScalar *tmp; - _try(VecGetArray(_b, &tmp)); + _check(VecGetArray(_b, &tmp)); PetscScalar s = tmp[row]; - _try(VecRestoreArray(_b, &tmp)); + _check(VecRestoreArray(_b, &tmp)); val = s.real(); } @@ -37,9 +37,9 @@ template <> void linearSystemPETSc<double>::getFromSolution(int row, double &val) const { PetscScalar *tmp; - _try(VecGetArray(_x, &tmp)); + _check(VecGetArray(_x, &tmp)); PetscScalar s = tmp[row]; - _try(VecRestoreArray(_x, &tmp)); + _check(VecRestoreArray(_x, &tmp)); val = s.real(); } #endif @@ -81,7 +81,7 @@ void linearSystemPETSc<fullMatrix<double> >::addToRightHandSide( { if(!_entriesPreAllocated) preAllocateEntries(); int blockSize; - _try(MatGetBlockSize(_a, &blockSize)); + _check(MatGetBlockSize(_a, &blockSize)); for(int ii = 0; ii < blockSize; ii++) { PetscInt i = row * blockSize + ii; PetscScalar v = val(ii, 0); @@ -94,7 +94,7 @@ void linearSystemPETSc<fullMatrix<double> >::getFromRightHandSide( int row, fullMatrix<double> &val) const { int blockSize; - _try(MatGetBlockSize(_a, &blockSize)); + _check(MatGetBlockSize(_a, &blockSize)); for(int i = 0; i < blockSize; i++) { int ii = row * blockSize + i; #ifdef PETSC_USE_COMPLEX @@ -112,7 +112,7 @@ void linearSystemPETSc<fullMatrix<double> >::addToSolution( int row, const fullMatrix<double> &val) { int blockSize; - _try(MatGetBlockSize(_a, &blockSize)); + _check(MatGetBlockSize(_a, &blockSize)); for(int ii = 0; ii < blockSize; ii++) { PetscInt i = row * blockSize + ii; PetscScalar v = val(ii, 0); @@ -125,7 +125,7 @@ void linearSystemPETSc<fullMatrix<double> >::getFromSolution( int row, fullMatrix<double> &val) const { int blockSize; - _try(MatGetBlockSize(_a, &blockSize)); + _check(MatGetBlockSize(_a, &blockSize)); for(int i = 0; i < blockSize; i++) { int ii = row * blockSize + i; #ifdef PETSC_USE_COMPLEX diff --git a/Solver/linearSystemPETSc.hpp b/Solver/linearSystemPETSc.hpp index 818e730ba56ebf3e897028d4a489458931f1e703..c9754ccbc8da0869d40125012b03da283fc2ab02 100644 --- a/Solver/linearSystemPETSc.hpp +++ b/Solver/linearSystemPETSc.hpp @@ -1,18 +1,21 @@ +// Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// issues on https://gitlab.onelab.info/gmsh/gmsh/issues + #if defined(HAVE_PETSC) #include <petsc.h> #include <petscksp.h> #include "linearSystemPETSc.h" -#if ((PETSC_VERSION_RELEASE == 0) || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 7))) +#if((PETSC_VERSION_RELEASE == 0) || \ + ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 7))) #define PetscOptionsGetInt(A, B, C, D) PetscOptionsGetInt(NULL, A, B, C, D) #define PetscOptionsInsertString(A) PetscOptionsInsertString(NULL, A) #endif -static void _try(int ierr) -{ - CHKERRABORT(PETSC_COMM_WORLD, ierr); -} +static void _check(int ierr) { CHKERRABORT(PETSC_COMM_WORLD, ierr); } template <class scalar> int linearSystemPETSc<scalar>::_getBlockSizeFromParameters() const @@ -20,26 +23,30 @@ int linearSystemPETSc<scalar>::_getBlockSizeFromParameters() const return 1; } -template <class scalar> -void linearSystemPETSc<scalar>::_kspCreate() +template <class scalar> void linearSystemPETSc<scalar>::_kspCreate() { - // Set option given by the user in its (python script) without using argc,argv or .petscrc + // Set option given by the user in its (python script) without using argc,argv + // or .petscrc if(this->_parameters.count("petsc_solver_options")) - _try(PetscOptionsInsertString(this->_parameters["petsc_solver_options"].c_str())); - _try(KSPCreate(_comm, &_ksp)); + _check(PetscOptionsInsertString( + this->_parameters["petsc_solver_options"].c_str())); + _check(KSPCreate(_comm, &_ksp)); PC pc; - _try(KSPGetPC(_ksp, &pc)); + _check(KSPGetPC(_ksp, &pc)); // set some default options - //_try(PCSetType(pc, PCLU));//LU for direct solver and PCILU for indirect solver - /* _try(PCFactorSetMatOrderingType(pc, MATORDERING_RCM)); - _try(PCFactorSetLevels(pc, 1));*/ - _try(KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); + //_check(PCSetType(pc, PCLU));//LU for direct solver and PCILU for indirect + //solver + /* _check(PCFactorSetMatOrderingType(pc, MATORDERING_RCM)); + _check(PCFactorSetLevels(pc, 1));*/ + _check( + KSPSetTolerances(_ksp, 1.e-8, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); // override the default options with the ones from the option // database (if any) - if (this->_parameters.count("petscPrefix")) - _try(KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str())); - _try(KSPSetFromOptions(_ksp)); - _try(PCSetFromOptions(pc)); + if(this->_parameters.count("petscPrefix")) + _check( + KSPAppendOptionsPrefix(_ksp, this->_parameters["petscPrefix"].c_str())); + _check(KSPSetFromOptions(_ksp)); + _check(PCSetFromOptions(pc)); _kspAllocated = true; } @@ -47,7 +54,7 @@ template <class scalar> int linearSystemPETSc<scalar>::getNumKspIteration() const { int n; - _try(KSPGetIterationNumber(_ksp, &n)); + _check(KSPGetIterationNumber(_ksp, &n)); return n; } @@ -62,8 +69,7 @@ linearSystemPETSc<scalar>::linearSystemPETSc(MPI_Comm com) _entriesPreAllocated = false; } -template <class scalar> -linearSystemPETSc<scalar>::linearSystemPETSc() +template <class scalar> linearSystemPETSc<scalar>::linearSystemPETSc() { _comm = PETSC_COMM_WORLD; _isAllocated = false; @@ -73,132 +79,137 @@ linearSystemPETSc<scalar>::linearSystemPETSc() _entriesPreAllocated = false; } -template <class scalar> -linearSystemPETSc<scalar>::~linearSystemPETSc() +template <class scalar> linearSystemPETSc<scalar>::~linearSystemPETSc() { clear(); - if(_kspAllocated) - _try(KSPDestroy(&_ksp)); + if(_kspAllocated) _check(KSPDestroy(&_ksp)); } template <class scalar> void linearSystemPETSc<scalar>::_assembleMatrixIfNeeded() { #if defined(HAVE_MPI) - if (_comm == PETSC_COMM_WORLD){ - if (Msg::GetCommSize() > 1){ + if(_comm == PETSC_COMM_WORLD) { + if(Msg::GetCommSize() > 1) { int value = _valuesNotAssembled ? 1 : 0; int sumValue = 0; - MPI_Allreduce((void*)&value, (void*)&sumValue, 1, MPI_INT, MPI_SUM, _comm); - if ((sumValue > 0) &&(sumValue < Msg::GetCommSize())){ + MPI_Allreduce((void *)&value, (void *)&sumValue, 1, MPI_INT, MPI_SUM, + _comm); + if((sumValue > 0) && (sumValue < Msg::GetCommSize())) { _valuesNotAssembled = 1; } } } #endif - if (_valuesNotAssembled) { - _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY)); - _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY)); + if(_valuesNotAssembled) { + _check(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY)); + _check(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY)); _matrixChangedSinceLastSolve = true; _valuesNotAssembled = false; } } template <class scalar> -void linearSystemPETSc<scalar>::insertInSparsityPattern (int i, int j) +void linearSystemPETSc<scalar>::insertInSparsityPattern(int i, int j) { i -= _localRowStart; - if (i<0 || i>= _localSize) return; - _sparsity.insertEntry (i,j); + if(i < 0 || i >= _localSize) return; + _sparsity.insertEntry(i, j); } -template <class scalar> -void linearSystemPETSc<scalar>::preAllocateEntries() +template <class scalar> void linearSystemPETSc<scalar>::preAllocateEntries() { - if (_entriesPreAllocated) return; - if (!_isAllocated) Msg::Fatal("system must be allocated first"); + if(_entriesPreAllocated) return; + if(!_isAllocated) Msg::Fatal("system must be allocated first"); int blockSize = _getBlockSizeFromParameters(); - std::vector<int> nByRowDiag (_localSize), nByRowOffDiag (_localSize); - if (_sparsity.getNbRows() == 0) { + std::vector<int> nByRowDiag(_localSize), nByRowOffDiag(_localSize); + if(_sparsity.getNbRows() == 0) { PetscInt prealloc = 500; PetscBool set; PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set); prealloc = std::min(prealloc, _localSize); nByRowDiag.resize(0); nByRowDiag.resize(_localSize, prealloc); - } else { - for (int i = 0; i < _localSize; i++) { + } + else { + for(int i = 0; i < _localSize; i++) { int n; const int *r = _sparsity.getRow(i, n); - for (int j = 0; j < n; j++) { - if (r[j] >= _localRowStart && r[j] < _localRowEnd) - nByRowDiag[i] ++; + for(int j = 0; j < n; j++) { + if(r[j] >= _localRowStart && r[j] < _localRowEnd) + nByRowDiag[i]++; else - nByRowOffDiag[i] ++; + nByRowOffDiag[i]++; } } _sparsity.clear(); } - //MatXAIJSetPreallocation is not available in petsc < 3.3 + // MatXAIJSetPreallocation is not available in petsc < 3.3 int commSize = 1; MPI_Comm_size(_comm, &commSize); - if (commSize == 1){ - if (blockSize == 1) - _try(MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0])); + if(commSize == 1) { + if(blockSize == 1) + _check(MatSeqAIJSetPreallocation(_a, 0, &nByRowDiag[0])); else - _try(MatSeqBAIJSetPreallocation(_a, blockSize, 0, &nByRowDiag[0])); + _check(MatSeqBAIJSetPreallocation(_a, blockSize, 0, &nByRowDiag[0])); } else { - if (blockSize == 1) - _try(MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0])); + if(blockSize == 1) + _check( + MatMPIAIJSetPreallocation(_a, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0])); else - _try(MatMPIBAIJSetPreallocation(_a, blockSize, 0, &nByRowDiag[0], 0, &nByRowOffDiag[0])); + _check(MatMPIBAIJSetPreallocation(_a, blockSize, 0, &nByRowDiag[0], 0, + &nByRowOffDiag[0])); } - if (blockSize > 1) - _try(MatSetOption(_a, MAT_ROW_ORIENTED, PETSC_FALSE)); + if(blockSize > 1) _check(MatSetOption(_a, MAT_ROW_ORIENTED, PETSC_FALSE)); _entriesPreAllocated = true; } -template <class scalar> -void linearSystemPETSc<scalar>::allocate(int nbRows) +template <class scalar> void linearSystemPETSc<scalar>::allocate(int nbRows) { int commSize; MPI_Comm_size(_comm, &commSize); int blockSize = _getBlockSizeFromParameters(); clear(); - _try(MatCreate(_comm, &_a)); - _try(MatSetSizes(_a, blockSize * nbRows, blockSize * nbRows, PETSC_DETERMINE, PETSC_DETERMINE)); - if (blockSize > 1) { - if (commSize > 1) { - _try(MatSetType(_a, MATMPIBAIJ)); + _check(MatCreate(_comm, &_a)); + _check(MatSetSizes(_a, blockSize * nbRows, blockSize * nbRows, + PETSC_DETERMINE, PETSC_DETERMINE)); + if(blockSize > 1) { + if(commSize > 1) { + _check(MatSetType(_a, MATMPIBAIJ)); } else { - _try(MatSetType(_a, MATSEQBAIJ)); + _check(MatSetType(_a, MATSEQBAIJ)); } } // override the default options with the ones from the option // database (if any) - if (this->_parameters.count("petscOptions")) - _try(PetscOptionsInsertString(this->_parameters["petscOptions"].c_str())); - if (this->_parameters.count("petscPrefix")) - _try(MatAppendOptionsPrefix(_a, this->_parameters["petscPrefix"].c_str())); - _try(MatSetFromOptions(_a)); - //since PETSc 3.3 GetOwnershipRange and MatGetSize cannot be called before MatXXXSetPreallocation + if(this->_parameters.count("petscOptions")) + _check(PetscOptionsInsertString(this->_parameters["petscOptions"].c_str())); + if(this->_parameters.count("petscPrefix")) + _check( + MatAppendOptionsPrefix(_a, this->_parameters["petscPrefix"].c_str())); + _check(MatSetFromOptions(_a)); + // since PETSc 3.3 GetOwnershipRange and MatGetSize cannot be called before + // MatXXXSetPreallocation _localSize = nbRows; #if defined(HAVE_MPI) - if (commSize > 1){ + if(commSize > 1) { _localRowStart = 0; - if (Msg::GetCommRank() != 0) { + if(Msg::GetCommRank() != 0) { MPI_Status status; - MPI_Recv((void*)&_localRowStart, 1, MPI_INT, Msg::GetCommRank() - 1, 1, MPI_COMM_WORLD, &status); + MPI_Recv((void *)&_localRowStart, 1, MPI_INT, Msg::GetCommRank() - 1, 1, + MPI_COMM_WORLD, &status); } _localRowEnd = _localRowStart + nbRows; - if (Msg::GetCommRank() != Msg::GetCommSize() - 1) { - MPI_Send((void*)&_localRowEnd, 1, MPI_INT, Msg::GetCommRank() + 1, 1, MPI_COMM_WORLD); + if(Msg::GetCommRank() != Msg::GetCommSize() - 1) { + MPI_Send((void *)&_localRowEnd, 1, MPI_INT, Msg::GetCommRank() + 1, 1, + MPI_COMM_WORLD); } - MPI_Allreduce((void*)&_localSize, (void*)&_globalSize, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce((void *)&_localSize, (void *)&_globalSize, 1, MPI_INT, + MPI_SUM, MPI_COMM_WORLD); } - else{ + else { _localRowStart = 0; _localRowEnd = nbRows; _globalSize = _localSize; @@ -209,78 +220,78 @@ void linearSystemPETSc<scalar>::allocate(int nbRows) _globalSize = _localSize; #endif // preallocation option must be set after other options - _try(VecCreate(_comm, &_x)); - _try(VecSetSizes(_x, blockSize * nbRows, PETSC_DETERMINE)); + _check(VecCreate(_comm, &_x)); + _check(VecSetSizes(_x, blockSize * nbRows, PETSC_DETERMINE)); // override the default options with the ones from the option // database (if any) - if (this->_parameters.count("petscPrefix")) - _try(VecAppendOptionsPrefix(_x, this->_parameters["petscPrefix"].c_str())); - _try(VecSetFromOptions(_x)); - _try(VecDuplicate(_x, &_b)); + if(this->_parameters.count("petscPrefix")) + _check( + VecAppendOptionsPrefix(_x, this->_parameters["petscPrefix"].c_str())); + _check(VecSetFromOptions(_x)); + _check(VecDuplicate(_x, &_b)); _isAllocated = true; _entriesPreAllocated = false; } -template <class scalar> -void linearSystemPETSc<scalar>::print() +template <class scalar> void linearSystemPETSc<scalar>::print() { _assembleMatrixIfNeeded(); - _try(VecAssemblyBegin(_b)); - _try(VecAssemblyEnd(_b)); + _check(VecAssemblyBegin(_b)); + _check(VecAssemblyEnd(_b)); /* PetscViewer fd; - _try(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "mat.m", &fd)); - _try(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB)); - _try(PetscObjectSetName((PetscObject)_a, "A")); - _try(MatView(_a, fd)); + _check(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "mat.m", &fd)); + _check(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB)); + _check(PetscObjectSetName((PetscObject)_a, "A")); + _check(MatView(_a, fd)); */ - if(Msg::GetCommRank()==0) - printf("a :\n"); + if(Msg::GetCommRank() == 0) printf("a :\n"); MatView(_a, PETSC_VIEWER_STDOUT_WORLD); - if(Msg::GetCommRank()==0) - printf("b :\n"); + if(Msg::GetCommRank() == 0) printf("b :\n"); VecView(_b, PETSC_VIEWER_STDOUT_WORLD); - if(Msg::GetCommRank()==0) - printf("x :\n"); + if(Msg::GetCommRank() == 0) printf("x :\n"); VecView(_x, PETSC_VIEWER_STDOUT_WORLD); } -template <class scalar> -void linearSystemPETSc<scalar>::clear() +template <class scalar> void linearSystemPETSc<scalar>::clear() { - if(_isAllocated){ - _try(MatDestroy(&_a)); - _try(VecDestroy(&_x)); - _try(VecDestroy(&_b)); + if(_isAllocated) { + _check(MatDestroy(&_a)); + _check(VecDestroy(&_x)); + _check(VecDestroy(&_b)); } _isAllocated = false; } template <class scalar> -void linearSystemPETSc<scalar>::getFromMatrix(int row, int col, scalar &val) const +void linearSystemPETSc<scalar>::getFromMatrix(int row, int col, + scalar &val) const { Msg::Error("getFromMatrix not implemented for PETSc"); } template <class scalar> -void linearSystemPETSc<scalar>::addToRightHandSide(int row, const scalar &val, int ith) +void linearSystemPETSc<scalar>::addToRightHandSide(int row, const scalar &val, + int ith) { PetscInt i = row; PetscScalar s = val; - _try(VecSetValues(_b, 1, &i, &s, ADD_VALUES)); + _check(VecSetValues(_b, 1, &i, &s, ADD_VALUES)); } #if defined(PETSC_USE_COMPLEX) -// this specialization will cast to a double (take the real part) if "val" is a double wheras Petsc is built in complex mode. +// this specialization will cast to a double (take the real part) if "val" is a +// double wheras Petsc is built in complex mode. template <> -void linearSystemPETSc<double>::getFromRightHandSide(int row, double &val) const; +void linearSystemPETSc<double>::getFromRightHandSide(int row, + double &val) const; #endif template <class scalar> void linearSystemPETSc<scalar>::getFromRightHandSide(int row, scalar &val) const { - _try(VecGetValues(_b, 1, &row, &val)); + _check(VecGetValues(_b, 1, &row, &val)); } template <class scalar> @@ -289,18 +300,17 @@ double linearSystemPETSc<scalar>::normInfRightHandSide() const PetscReal nor; VecAssemblyBegin(_b); VecAssemblyEnd(_b); - _try(VecNorm(_b, NORM_INFINITY, &nor)); + _check(VecNorm(_b, NORM_INFINITY, &nor)); return nor; } template <class scalar> void linearSystemPETSc<scalar>::addToMatrix(int row, int col, const scalar &val) { - if (!_entriesPreAllocated) - preAllocateEntries(); + if(!_entriesPreAllocated) preAllocateEntries(); PetscInt i = row, j = col; PetscScalar s = val; - _try(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES)); + _check(MatSetValues(_a, 1, &i, 1, &j, &s, ADD_VALUES)); _valuesNotAssembled = true; } @@ -309,10 +319,11 @@ void linearSystemPETSc<scalar>::addToSolution(int row, const scalar &val) { PetscInt i = row; PetscScalar s = val; - _try(VecSetValues(_x, 1, &i, &s, ADD_VALUES)); + _check(VecSetValues(_x, 1, &i, &s, ADD_VALUES)); } #if defined(PETSC_USE_COMPLEX) -// this specialization will cast to a double (take the real part) if "val" is a double wheras Petsc is built in complex mode. +// this specialization will cast to a double (take the real part) if "val" is a +// double wheras Petsc is built in complex mode. template <> void linearSystemPETSc<double>::getFromSolution(int row, double &val) const; #endif @@ -320,100 +331,96 @@ void linearSystemPETSc<double>::getFromSolution(int row, double &val) const; template <class scalar> void linearSystemPETSc<scalar>::getFromSolution(int row, scalar &val) const { - _try(VecGetValues(_x, 1, &row, &val)); + _check(VecGetValues(_x, 1, &row, &val)); } -template <class scalar> -void linearSystemPETSc<scalar>::zeroMatrix() +template <class scalar> void linearSystemPETSc<scalar>::zeroMatrix() { #if defined(HAVE_MPI) - if (_comm == PETSC_COMM_WORLD){ - if (Msg::GetCommSize() > 1){ + if(_comm == PETSC_COMM_WORLD) { + if(Msg::GetCommSize() > 1) { int value = _entriesPreAllocated ? 1 : 0; int sumValue = 0; - MPI_Allreduce((void*)&value, (void*)&sumValue, 1, MPI_INT, MPI_SUM, _comm); - if ((sumValue >= 0) &&(sumValue < Msg::GetCommSize()) && !_entriesPreAllocated){ + MPI_Allreduce((void *)&value, (void *)&sumValue, 1, MPI_INT, MPI_SUM, + _comm); + if((sumValue >= 0) && (sumValue < Msg::GetCommSize()) && + !_entriesPreAllocated) { preAllocateEntries(); } } } #endif - if (_isAllocated && _entriesPreAllocated) { + if(_isAllocated && _entriesPreAllocated) { _assembleMatrixIfNeeded(); - _try(MatZeroEntries(_a)); + _check(MatZeroEntries(_a)); } } -template <class scalar> -void linearSystemPETSc<scalar>::zeroRightHandSide() +template <class scalar> void linearSystemPETSc<scalar>::zeroRightHandSide() { - if (_isAllocated) { - _try(VecAssemblyBegin(_b)); - _try(VecAssemblyEnd(_b)); - _try(VecZeroEntries(_b)); + if(_isAllocated) { + _check(VecAssemblyBegin(_b)); + _check(VecAssemblyEnd(_b)); + _check(VecZeroEntries(_b)); } } -template <class scalar> -void linearSystemPETSc<scalar>::zeroSolution() +template <class scalar> void linearSystemPETSc<scalar>::zeroSolution() { - if (_isAllocated) { - _try(VecAssemblyBegin(_x)); - _try(VecAssemblyEnd(_x)); - _try(VecZeroEntries(_x)); + if(_isAllocated) { + _check(VecAssemblyBegin(_x)); + _check(VecAssemblyEnd(_x)); + _check(VecZeroEntries(_x)); } } - -template <class scalar> -int linearSystemPETSc<scalar>::systemSolve() +template <class scalar> int linearSystemPETSc<scalar>::systemSolve() { - if (!_kspAllocated) - _kspCreate(); + if(!_kspAllocated) _kspCreate(); _assembleMatrixIfNeeded(); - if (!_matrixChangedSinceLastSolve || linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix") - _try(KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER)); - else if (linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity") - _try(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN)); + if(!_matrixChangedSinceLastSolve || + linearSystem<scalar>::_parameters["matrix_reuse"] == "same_matrix") + _check(KSPSetOperators(_ksp, _a, _a, SAME_PRECONDITIONER)); + else if(linearSystem<scalar>::_parameters["matrix_reuse"] == "same_sparsity") + _check(KSPSetOperators(_ksp, _a, _a, SAME_NONZERO_PATTERN)); else - _try(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN)); + _check(KSPSetOperators(_ksp, _a, _a, DIFFERENT_NONZERO_PATTERN)); _matrixChangedSinceLastSolve = false; /*MatInfo info; MatGetInfo(_a, MAT_LOCAL, &info); - printf("mallocs %.0f nz_allocated %.0f nz_used %.0f nz_unneeded %.0f\n", info.mallocs, info.nz_allocated, info.nz_used, info.nz_unneeded);*/ - _try(VecAssemblyBegin(_b)); - _try(VecAssemblyEnd(_b)); - _try(KSPSolve(_ksp, _b, _x)); - //_try(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF)); - //PetscInt its; - //_try(KSPGetIterationNumber(ksp, &its)); - //Msg::Info("%d iterations", its); + printf("mallocs %.0f nz_allocated %.0f nz_used %.0f nz_unneeded + %.0f\n", info.mallocs, info.nz_allocated, info.nz_used, info.nz_unneeded);*/ + _check(VecAssemblyBegin(_b)); + _check(VecAssemblyEnd(_b)); + _check(KSPSolve(_ksp, _b, _x)); + //_check(KSPView(ksp, PETSC_VIEWER_STDOUT_SELF)); + // PetscInt its; + //_check(KSPGetIterationNumber(ksp, &its)); + // Msg::Info("%d iterations", its); return 1; } -template <class scalar> -int linearSystemPETSc<scalar>::matMult() +template <class scalar> int linearSystemPETSc<scalar>::matMult() { - _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY)); - _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY)); - _try(VecAssemblyBegin(_b)); - _try(VecAssemblyEnd(_b)); - _try(MatMult(_a,_b,_x)); + _check(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY)); + _check(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY)); + _check(VecAssemblyBegin(_b)); + _check(VecAssemblyEnd(_b)); + _check(MatMult(_a, _b, _x)); return 1; } - template <class scalar> void linearSystemPETSc<scalar>::printMatlab(const char *filename) const { - _try(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY)); - _try(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY)); - _try(VecAssemblyBegin(_b)); - _try(VecAssemblyEnd(_b)); + _check(MatAssemblyBegin(_a, MAT_FINAL_ASSEMBLY)); + _check(MatAssemblyEnd(_a, MAT_FINAL_ASSEMBLY)); + _check(VecAssemblyBegin(_b)); + _check(VecAssemblyEnd(_b)); PetscViewer viewer; PetscViewerASCIIOpen(PETSC_COMM_WORLD, filename, &viewer); - PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); + PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB); MatView(_a, viewer); PetscViewerDestroy(&viewer); return; @@ -424,9 +431,9 @@ std::vector<scalar> linearSystemPETSc<scalar>::getData() { _assembleMatrixIfNeeded(); PetscScalar *v; - _try(MatGetArray(_a,&v)); + _check(MatGetArray(_a,&v)); MatInfo info; - _try(MatGetInfo(_a,MAT_LOCAL,&info)); + _check(MatGetInfo(_a,MAT_LOCAL,&info)); std::vector<scalar> data; // Maybe I should reserve or resize (SAM) #if defined(PETSC_USE_COMPLEX) @@ -436,7 +443,7 @@ std::vector<scalar> linearSystemPETSc<scalar>::getData() for (int i = 0; i < info.nz_allocated; i++) data.push_back(v[i]); #endif - _try(MatRestoreArray(_a,&v)); + _check(MatRestoreArray(_a,&v)); return data; }*/ /* @@ -448,11 +455,11 @@ std::vector<int> linearSystemPETSc<scalar>::getRowPointers() const PetscInt *columns; PetscInt n; PetscBool done; - _try(MatGetRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); //case done == PETSC_FALSE should be handled - std::vector<int> rowPointers; // Maybe I should reserve or resize (SAM) - for (int i = 0; i <= n; i++) + _check(MatGetRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); +//case done == PETSC_FALSE should be handled std::vector<int> rowPointers; // +Maybe I should reserve or resize (SAM) for (int i = 0; i <= n; i++) rowPointers.push_back(rows[i]); - _try(MatRestoreRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); + _check(MatRestoreRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); return rowPointers; } @@ -464,13 +471,13 @@ std::vector<int> linearSystemPETSc<scalar>::getColumnsIndices() const PetscInt *columns; PetscInt n; PetscBool done; - _try(MatGetRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); //case done == PETSC_FALSE should be handled - MatInfo info; - _try(MatGetInfo(_a,MAT_LOCAL,&info)); + _check(MatGetRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); +//case done == PETSC_FALSE should be handled MatInfo info; + _check(MatGetInfo(_a,MAT_LOCAL,&info)); std::vector<int> columnIndices; // Maybe I should reserve or resize (SAM) for (int i = 0; i < info.nz_allocated; i++) columnIndices.push_back(columns[i]); - _try(MatRestoreRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); + _check(MatRestoreRowIJ(_a,0,PETSC_FALSE,PETSC_FALSE,&n,&rows,&columns,&done)); return columnIndices; } */