is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| LinAlg_PETSC.cpp (getdp-3.4.0-source.tgz) | : | LinAlg_PETSC.cpp (getdp-3.5.0-source.tgz) |
|---|
|
| // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege | | // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege |
| // | | // |
| // See the LICENSE.txt file for license information. Please report all | | // See the LICENSE.txt file for license information. Please report all |
| // issues on https://gitlab.onelab.info/getdp/getdp/issues. | | // issues on https://gitlab.onelab.info/getdp/getdp/issues. |
| // | | // |
| // Contributor(s): | | // Contributor(s): |
| // David Colignon | | // David Colignon |
| // Ruth Sabariego | | // Ruth Sabariego |
| // Jose Geraldo A. Brito Neto | | // Jose Geraldo A. Brito Neto |
| // | | // |
| | |
| | |
| skipping to change at line 28 | | skipping to change at line 28 |
|---|
| #include "GetDPConfig.h" | | #include "GetDPConfig.h" |
| #include "LinAlg.h" | | #include "LinAlg.h" |
| #include "MallocUtils.h" | | #include "MallocUtils.h" |
| #include "Message.h" | | #include "Message.h" |
| #include "OS.h" | | #include "OS.h" |
| | |
| #if defined(HAVE_SLEPC) | | #if defined(HAVE_SLEPC) |
| #include <slepc.h> | | #include <slepc.h> |
| #endif | | #endif |
| | |
|
| // Johan, we curse you for a thousand generations! | | // FIXME: this dependency should be removed |
| #include "ProData.h" | | #include "ProData.h" |
| #include "DofData.h" | | #include "DofData.h" |
|
| extern struct CurrentData Current ; | | extern struct CurrentData Current; |
| | |
| #if defined(HAVE_PETSC) | | #if defined(HAVE_PETSC) |
| | |
| // Options for PETSc can be provided on the command line, or in the file | | // Options for PETSc can be provided on the command line, or in the file |
| // ~/.petscrc. | | // ~/.petscrc. |
| // | | // |
| // By default we use a direct solver (MUMPS, UMFPACK or the PETSc LU). | | // By default we use a direct solver (MUMPS, UMFPACK or the PETSc LU). |
| // | | // |
| // All these options can be changed at runtime. For example you could | | // All these options can be changed at runtime. For example you could |
| // use | | // use |
| | |
| skipping to change at line 59 | | skipping to change at line 59 |
|---|
| // | | // |
| // for GMRES with ILU(0), with a restart of 500 and a stopping | | // for GMRES with ILU(0), with a restart of 500 and a stopping |
| // criterion of 1e-6. | | // criterion of 1e-6. |
| | |
| static MPI_Comm MyComm = MPI_COMM_SELF; | | static MPI_Comm MyComm = MPI_COMM_SELF; |
| static PetscViewer MyPetscViewer; | | static PetscViewer MyPetscViewer; |
| | |
| static void _try(int ierr) | | static void _try(int ierr) |
| { | | { |
| CHKERRCONTINUE(ierr); | | CHKERRCONTINUE(ierr); |
|
| if(PetscUnlikely(ierr)){ | | if(PetscUnlikely(ierr)) { |
| const char *text; | | const char *text; |
| PetscErrorMessage(ierr, &text, 0); | | PetscErrorMessage(ierr, &text, 0); |
| // Do not produce an error in case of a PETSc-crash when we are in | | // Do not produce an error in case of a PETSc-crash when we are in |
| // TimeLoopAdaptive loop | | // TimeLoopAdaptive loop |
|
| if (Message::GetOperatingInTimeLoopAdaptive()) | | if(Message::GetOperatingInTimeLoopAdaptive()) |
| Message::Warning("PETSc error: %s", text); | | Message::Warning("PETSc error: %s", text); |
| else | | else |
| Message::Error("PETSc error: %s", text); | | Message::Error("PETSc error: %s", text); |
| Message::SetLastPETScError(ierr); | | Message::SetLastPETScError(ierr); |
| } | | } |
| } | | } |
| | |
| static int SolverInitialized = 0; | | static int SolverInitialized = 0; |
| | |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 7) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 7) |
| #define PetscTruth PetscBool | | #define PetscTruth PetscBool |
| #define PetscOptionsGetTruth(A, B, C, D) PetscOptionsGetBool(A, NULL, B, C, D) | | #define PetscOptionsGetTruth(A, B, C, D) PetscOptionsGetBool(A, NULL, B, C, D) |
| #define PetscOptionsInsertFile(A, B, C) PetscOptionsInsertFile(A, NULL, B, C) | | #define PetscOptionsInsertFile(A, B, C) PetscOptionsInsertFile(A, NULL, B, C) |
| #define PetscOptionsGetInt(A, B, C, D) PetscOptionsGetInt(NULL, A, B, C, D); | | #define PetscOptionsGetInt(A, B, C, D) PetscOptionsGetInt(NULL, A, B, C, D); |
| #define PetscOptionsSetValue(A, B) PetscOptionsSetValue(NULL, A, B) | | #define PetscOptionsSetValue(A, B) PetscOptionsSetValue(NULL, A, B) |
| #define PetscOptionsInsertString(A) PetscOptionsInsertString(NULL, A) | | #define PetscOptionsInsertString(A) PetscOptionsInsertString(NULL, A) |
| #define PetscViewerSetFormat(A, B) PetscViewerPushFormat(A, B) | | #define PetscViewerSetFormat(A, B) PetscViewerPushFormat(A, B) |
|
| #elif ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)) | | #elif((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2)) |
| #define PetscTruth PetscBool | | #define PetscTruth PetscBool |
| #define PetscOptionsGetTruth PetscOptionsGetBool | | #define PetscOptionsGetTruth PetscOptionsGetBool |
| #endif | | #endif |
| | |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 9) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 9) |
| #define PCFactorSetMatSolverPackage PCFactorSetMatSolverType | | #define PCFactorSetMatSolverPackage PCFactorSetMatSolverType |
| #endif | | #endif |
| | |
|
| void LinAlg_InitializeSolver(int* argc, char*** argv) | | void LinAlg_InitializeSolver(int *argc, char ***argv) |
| { | | { |
| if(SolverInitialized) return; | | if(SolverInitialized) return; |
| SolverInitialized = 1; | | SolverInitialized = 1; |
| | |
| // This function detects if MPI is initialized | | // This function detects if MPI is initialized |
| PetscInitialize(argc, argv, PETSC_NULL, PETSC_NULL); | | PetscInitialize(argc, argv, PETSC_NULL, PETSC_NULL); |
| PetscPopSignalHandler(); | | PetscPopSignalHandler(); |
| MyPetscViewer = PETSC_VIEWER_STDOUT_SELF; | | MyPetscViewer = PETSC_VIEWER_STDOUT_SELF; |
| MyComm = PETSC_COMM_WORLD; | | MyComm = PETSC_COMM_WORLD; |
| #if defined(HAVE_SLEPC) | | #if defined(HAVE_SLEPC) |
| SlepcInitialize(argc, argv, PETSC_NULL, PETSC_NULL); | | SlepcInitialize(argc, argv, PETSC_NULL, PETSC_NULL); |
| #endif | | #endif |
| | |
| // get additional petsc options from specified file (useful e.g. on | | // get additional petsc options from specified file (useful e.g. on |
| // Windows where we don't know where to search for ~/.petscrc) | | // Windows where we don't know where to search for ~/.petscrc) |
|
| for(int i = 0; i < *argc - 1; i++){ | | for(int i = 0; i < *argc - 1; i++) { |
| if (!strcmp((*argv)[i], "-solver")){ | | if(!strcmp((*argv)[i], "-solver")) { |
| #if (PETSC_VERSION_MAJOR == 2) | | #if(PETSC_VERSION_MAJOR == 2) |
| PetscOptionsInsertFile((*argv)[i+1]); | | PetscOptionsInsertFile((*argv)[i + 1]); |
| #else | | #else |
|
| PetscOptionsInsertFile(MyComm, (*argv)[i+1], PETSC_FALSE); | | PetscOptionsInsertFile(MyComm, (*argv)[i + 1], PETSC_FALSE); |
| #endif | | #endif |
| } | | } |
| } | | } |
| } | | } |
| | |
| void LinAlg_FinalizeSolver() | | void LinAlg_FinalizeSolver() |
| { | | { |
|
| if(SolverInitialized){ | | if(SolverInitialized) { |
| #if defined(HAVE_SLEPC) | | #if defined(HAVE_SLEPC) |
| SlepcFinalize(); | | SlepcFinalize(); |
| #endif | | #endif |
| PetscFinalize(); | | PetscFinalize(); |
| SolverInitialized = 0; | | SolverInitialized = 0; |
| } | | } |
| } | | } |
| | |
| void LinAlg_SetCommSelf() | | void LinAlg_SetCommSelf() |
| { | | { |
| | |
| skipping to change at line 145 | | skipping to change at line 145 |
|---|
| | |
| void LinAlg_SetCommWorld() | | void LinAlg_SetCommWorld() |
| { | | { |
| Message::Info("Set communicator to WORLD"); | | Message::Info("Set communicator to WORLD"); |
| MyComm = PETSC_COMM_WORLD; | | MyComm = PETSC_COMM_WORLD; |
| Message::SetIsCommWorld(1); | | Message::SetIsCommWorld(1); |
| } | | } |
| | |
| void LinAlg_CreateSolver(gSolver *Solver, const char *SolverDataFileName) | | void LinAlg_CreateSolver(gSolver *Solver, const char *SolverDataFileName) |
| { | | { |
|
| for(int i = 0; i < 10; i++){ | | for(int i = 0; i < 10; i++) { |
| Solver->ksp[i] = NULL; | | Solver->ksp[i] = NULL; |
| Solver->snes[i] = NULL; | | Solver->snes[i] = NULL; |
| } | | } |
| } | | } |
| | |
| void LinAlg_CreateVector(gVector *V, gSolver *Solver, int n) | | void LinAlg_CreateVector(gVector *V, gSolver *Solver, int n) |
| { | | { |
| _try(VecCreate(MyComm, &V->V)); | | _try(VecCreate(MyComm, &V->V)); |
| _try(VecSetSizes(V->V, PETSC_DECIDE, n)); | | _try(VecSetSizes(V->V, PETSC_DECIDE, n)); |
| | |
| // override the default options with the ones from the option | | // override the default options with the ones from the option |
| // database (if any) | | // database (if any) |
| _try(VecSetFromOptions(V->V)); | | _try(VecSetFromOptions(V->V)); |
| | |
| // create sequential vector that will contain all the values on all | | // create sequential vector that will contain all the values on all |
| // the procs | | // the procs |
|
| if(Message::GetCommSize() > 1 && MyComm != PETSC_COMM_SELF){ | | if(Message::GetCommSize() > 1 && MyComm != PETSC_COMM_SELF) { |
| _try(VecCreateSeq(PETSC_COMM_SELF, n, &V->Vseq)); | | _try(VecCreateSeq(PETSC_COMM_SELF, n, &V->Vseq)); |
| V->haveSeq = 1; | | V->haveSeq = 1; |
| } | | } |
|
| else{ | | else { |
| V->haveSeq = 0; | | V->haveSeq = 0; |
| } | | } |
| } | | } |
| | |
| void _fillseq(Vec &V, Vec &Vseq) | | void _fillseq(Vec &V, Vec &Vseq) |
| { | | { |
| // collect all the values from the parallel petsc vector into a sequential | | // collect all the values from the parallel petsc vector into a sequential |
| // vector on each processor | | // vector on each processor |
| VecScatter ctx; | | VecScatter ctx; |
| VecScatterCreateToAll(V, &ctx, NULL); | | VecScatterCreateToAll(V, &ctx, NULL); |
|
| #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR == 3) && (PETSC_VERSION_S | | #if(PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR == 3) && \ |
| UBMINOR < 3) | | (PETSC_VERSION_SUBMINOR < 3) |
| VecScatterBegin(V, Vseq, INSERT_VALUES, SCATTER_FORWARD, ctx); | | VecScatterBegin(V, Vseq, INSERT_VALUES, SCATTER_FORWARD, ctx); |
| VecScatterEnd(V, Vseq, INSERT_VALUES, SCATTER_FORWARD, ctx); | | VecScatterEnd(V, Vseq, INSERT_VALUES, SCATTER_FORWARD, ctx); |
| #else | | #else |
| VecScatterBegin(ctx, V, Vseq, INSERT_VALUES, SCATTER_FORWARD); | | VecScatterBegin(ctx, V, Vseq, INSERT_VALUES, SCATTER_FORWARD); |
| VecScatterEnd(ctx, V, Vseq, INSERT_VALUES, SCATTER_FORWARD); | | VecScatterEnd(ctx, V, Vseq, INSERT_VALUES, SCATTER_FORWARD); |
| #endif | | #endif |
| | |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) |
| VecScatterDestroy(&ctx); | | VecScatterDestroy(&ctx); |
| #else | | #else |
| VecScatterDestroy(ctx); | | VecScatterDestroy(ctx); |
| #endif | | #endif |
| } | | } |
| | |
| static void _fillseq(gVector *V) | | static void _fillseq(gVector *V) |
| { | | { |
| if(V->haveSeq) _fillseq(V->V, V->Vseq); | | if(V->haveSeq) _fillseq(V->V, V->Vseq); |
| } | | } |
| | |
| void LinAlg_CreateMatrix(gMatrix *M, gSolver *Solver, int n, int m) | | void LinAlg_CreateMatrix(gMatrix *M, gSolver *Solver, int n, int m) |
| { | | { |
|
| | | PetscInt prealloc = 100.; |
| PetscInt prealloc = 100. ; | | PetscInt prealloc_full = n; |
| PetscInt prealloc_full = n ; | | |
| int nonloc = Current.DofData->NonLocalEquations.size(); | | int nonloc = Current.DofData->NonLocalEquations.size(); |
| // heuristic for preallocation of global rows: don't prelloc more than 500 Mb | | // heuristic for preallocation of global rows: don't prelloc more than 500 Mb |
| double limit = 500. * 1024. * 1024. / (gSCALAR_SIZE * sizeof(double)); | | double limit = 500. * 1024. * 1024. / (gSCALAR_SIZE * sizeof(double)); |
| double estim = (double)nonloc * (double)n; | | double estim = (double)nonloc * (double)n; |
|
| if(estim > limit){ | | if(estim > limit) { |
| prealloc_full = (int)(limit / nonloc); | | prealloc_full = (int)(limit / nonloc); |
|
| Message::Debug("Heuristic -petsc_prealloc_full changed to %d", prealloc_full | | Message::Debug("Heuristic -petsc_prealloc_full changed to %d", |
| ); | | prealloc_full); |
| } | | } |
| | |
| PetscTruth set; | | PetscTruth set; |
| PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set); | | PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc", &prealloc, &set); |
| PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc_full", &prealloc_full, &set); | | PetscOptionsGetInt(PETSC_NULL, "-petsc_prealloc_full", &prealloc_full, &set); |
| | |
| // prealloc cannot be bigger than the number of rows! | | // prealloc cannot be bigger than the number of rows! |
| prealloc = (n < prealloc) ? n : prealloc; | | prealloc = (n < prealloc) ? n : prealloc; |
| prealloc_full = (n < prealloc_full) ? n : prealloc_full; | | prealloc_full = (n < prealloc_full) ? n : prealloc_full; |
| | |
| std::vector<PetscInt> nnz(n, prealloc); | | std::vector<PetscInt> nnz(n, prealloc); |
| | |
| // preallocate non local equations as full lines (this is not | | // preallocate non local equations as full lines (this is not |
| // optimal, but preallocating too few elements leads to horrible | | // optimal, but preallocating too few elements leads to horrible |
| // assembly performance: petsc really sucks at dynamic reallocation | | // assembly performance: petsc really sucks at dynamic reallocation |
| // in the AIJ matrix format) | | // in the AIJ matrix format) |
| for(int i = 0; i < nonloc; i++) | | for(int i = 0; i < nonloc; i++) |
| nnz[Current.DofData->NonLocalEquations[i] - 1] = prealloc_full; | | nnz[Current.DofData->NonLocalEquations[i] - 1] = prealloc_full; |
| | |
|
| if(Message::GetCommSize() > 1){ // FIXME: alloc full lines... | | if(Message::GetCommSize() > 1) { // FIXME: alloc full lines... |
| #if ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 3)) | | #if((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 3)) |
| _try(MatCreateAIJ(MyComm, PETSC_DECIDE, PETSC_DECIDE, n, m, | | _try(MatCreateAIJ(MyComm, PETSC_DECIDE, PETSC_DECIDE, n, m, prealloc, |
| prealloc, PETSC_NULL, prealloc, PETSC_NULL, &M->M)); | | PETSC_NULL, prealloc, PETSC_NULL, &M->M)); |
| #else | | #else |
|
| _try(MatCreateMPIAIJ(MyComm, PETSC_DECIDE, PETSC_DECIDE, n, m, | | _try(MatCreateMPIAIJ(MyComm, PETSC_DECIDE, PETSC_DECIDE, n, m, prealloc, |
| prealloc, PETSC_NULL, prealloc, PETSC_NULL, &M->M)); | | PETSC_NULL, prealloc, PETSC_NULL, &M->M)); |
| #endif | | #endif |
| } | | } |
|
| else{ | | else { |
| _try(MatCreateSeqAIJ(PETSC_COMM_SELF, n, m, 0, &nnz[0], &M->M)); | | _try(MatCreateSeqAIJ(PETSC_COMM_SELF, n, m, 0, &nnz[0], &M->M)); |
| // PETSc (I)LU does not like matrices with empty (non assembled) diagonals | | // PETSc (I)LU does not like matrices with empty (non assembled) diagonals |
| | |
|
| for(int i = 0; i < n; i++){ | | for(int i = 0; i < n; i++) { |
| PetscScalar d = 0.; | | PetscScalar d = 0.; |
| _try(MatSetValues(M->M, 1, &i, 1, &i, &d, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &i, 1, &i, &d, INSERT_VALUES)); |
| } | | } |
| _try(MatAssemblyBegin(M->M, MAT_FLUSH_ASSEMBLY)); | | _try(MatAssemblyBegin(M->M, MAT_FLUSH_ASSEMBLY)); |
| _try(MatAssemblyEnd(M->M, MAT_FLUSH_ASSEMBLY)); | | _try(MatAssemblyEnd(M->M, MAT_FLUSH_ASSEMBLY)); |
|
| | | |
| } | | } |
| | |
|
| //MatSetOption(M->M, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); | | // MatSetOption(M->M, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); |
| | |
|
| #if ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 3)) | | #if((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 3)) |
| // Preallocation routines automatically set now MAT_NEW_NONZERO_ALLOCATION_ERR | | // Preallocation routines automatically set now |
| , | | // MAT_NEW_NONZERO_ALLOCATION_ERR, what causes a problem when the mask of the |
| // what causes a problem when the mask of the matrix changes (e.g. moving band | | // matrix changes (e.g. moving band) We must disable the error generation and |
| ) | | // allow new allocation (if needed) |
| // We must disable the error generation and allow new allocation (if needed) | | _try(MatSetOption(M->M, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); |
| _try(MatSetOption(M->M,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); | | |
| #endif | | #endif |
| | |
| // override the default options with the ones from the option | | // override the default options with the ones from the option |
| // database (if any) | | // database (if any) |
| _try(MatSetFromOptions(M->M)); | | _try(MatSetFromOptions(M->M)); |
|
| | | |
| } | | } |
| | |
| void LinAlg_DestroySolver(gSolver *Solver) | | void LinAlg_DestroySolver(gSolver *Solver) |
| { | | { |
|
| for(int i = 0; i < 10; i++){ | | for(int i = 0; i < 10; i++) { |
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) |
| if(Solver->ksp[i]){ | | if(Solver->ksp[i]) { |
| _try(KSPDestroy(&Solver->ksp[i])); | | _try(KSPDestroy(&Solver->ksp[i])); |
| Solver->ksp[i] = NULL; | | Solver->ksp[i] = NULL; |
| } | | } |
|
| if(Solver->snes[i]){ | | if(Solver->snes[i]) { |
| _try(SNESDestroy(&Solver->snes[i])); | | _try(SNESDestroy(&Solver->snes[i])); |
| Solver->snes[i] = NULL; | | Solver->snes[i] = NULL; |
| } | | } |
| #else | | #else |
|
| if(Solver->ksp[i]){ | | if(Solver->ksp[i]) { |
| _try(KSPDestroy(Solver->ksp[i])); | | _try(KSPDestroy(Solver->ksp[i])); |
| Solver->ksp[i] = NULL; | | Solver->ksp[i] = NULL; |
| } | | } |
|
| if(Solver->snes[i]){ | | if(Solver->snes[i]) { |
| _try(SNESDestroy(Solver->snes[i])); | | _try(SNESDestroy(Solver->snes[i])); |
| Solver->snes[i] = NULL; | | Solver->snes[i] = NULL; |
| } | | } |
| #endif | | #endif |
| } | | } |
| } | | } |
| | |
| void LinAlg_DestroyVector(gVector *V) | | void LinAlg_DestroyVector(gVector *V) |
| { | | { |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) |
| _try(VecDestroy(&V->V)); | | _try(VecDestroy(&V->V)); |
| if(V->haveSeq) _try(VecDestroy(&V->Vseq)); | | if(V->haveSeq) _try(VecDestroy(&V->Vseq)); |
| #else | | #else |
| _try(VecDestroy(V->V)); | | _try(VecDestroy(V->V)); |
| if(V->haveSeq) _try(VecDestroy(V->Vseq)); | | if(V->haveSeq) _try(VecDestroy(V->Vseq)); |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_DestroyMatrix(gMatrix *M) | | void LinAlg_DestroyMatrix(gMatrix *M) |
| { | | { |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) |
| _try(MatDestroy(&M->M)); | | _try(MatDestroy(&M->M)); |
| #else | | #else |
| _try(MatDestroy(M->M)); | | _try(MatDestroy(M->M)); |
| #endif | | #endif |
| } | | } |
| | |
|
| void LinAlg_CopyScalar(gScalar *S1, gScalar *S2) | | void LinAlg_CopyScalar(gScalar *S1, gScalar *S2) { S1->s = S2->s; } |
| { | | |
| S1->s = S2->s; | | |
| } | | |
| | |
| void LinAlg_CopyVector(gVector *V1, gVector *V2) | | void LinAlg_CopyVector(gVector *V1, gVector *V2) |
| { | | { |
| _try(VecCopy(V1->V, V2->V)); | | _try(VecCopy(V1->V, V2->V)); |
| if(V1->haveSeq) _try(VecCopy(V1->Vseq, V2->Vseq)); | | if(V1->haveSeq) _try(VecCopy(V1->Vseq, V2->Vseq)); |
| } | | } |
| | |
| void LinAlg_SwapVector(gVector *V1, gVector *V2) | | void LinAlg_SwapVector(gVector *V1, gVector *V2) |
| { | | { |
| _try(VecSwap(V1->V, V2->V)); | | _try(VecSwap(V1->V, V2->V)); |
| if(V1->haveSeq) _try(VecSwap(V1->Vseq, V2->Vseq)); | | if(V1->haveSeq) _try(VecSwap(V1->Vseq, V2->Vseq)); |
| } | | } |
| | |
| void LinAlg_CopyMatrix(gMatrix *M1, gMatrix *M2) | | void LinAlg_CopyMatrix(gMatrix *M1, gMatrix *M2) |
| { | | { |
| _try(MatCopy(M1->M, M2->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatCopy(M1->M, M2->M, DIFFERENT_NONZERO_PATTERN)); |
| } | | } |
| | |
|
| void LinAlg_ZeroScalar(gScalar *S) | | void LinAlg_ZeroScalar(gScalar *S) { S->s = 0.; } |
| { | | |
| S->s = 0.; | | |
| } | | |
| | |
| void LinAlg_ZeroVector(gVector *V) | | void LinAlg_ZeroVector(gVector *V) |
| { | | { |
| PetscScalar zero = 0.0; | | PetscScalar zero = 0.0; |
| | |
| _try(VecSet(V->V, zero)); | | _try(VecSet(V->V, zero)); |
| if(V->haveSeq) _try(VecSet(V->Vseq, zero)); | | if(V->haveSeq) _try(VecSet(V->Vseq, zero)); |
| } | | } |
| | |
|
| void LinAlg_ZeroMatrix(gMatrix *M) | | void LinAlg_ZeroMatrix(gMatrix *M) { _try(MatZeroEntries(M->M)); } |
| { | | |
| _try(MatZeroEntries(M->M)); | | |
| } | | |
| | |
| void LinAlg_ScanScalar(FILE *file, gScalar *S) | | void LinAlg_ScanScalar(FILE *file, gScalar *S) |
| { | | { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| double a, b; | | double a, b; |
| if(fscanf(file, "%lf %lf", &a, &b) != 2) | | if(fscanf(file, "%lf %lf", &a, &b) != 2) |
| Message::Error("Could not scan scalar"); | | Message::Error("Could not scan scalar"); |
| S->s = a + PETSC_i * b; | | S->s = a + PETSC_i * b; |
| #else | | #else |
|
| if(fscanf(file, "%lf", &S->s) != 1) | | if(fscanf(file, "%lf", &S->s) != 1) Message::Error("Could not scan scalar"); |
| Message::Error("Could not scan scalar"); | | |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_ScanVector(FILE *file, gVector *V) | | void LinAlg_ScanVector(FILE *file, gVector *V) |
| { | | { |
| PetscInt n; | | PetscInt n; |
| _try(VecGetSize(V->V, &n)); | | _try(VecGetSize(V->V, &n)); |
|
| for(PetscInt i = 0; i < n; i++){ | | for(PetscInt i = 0; i < n; i++) { |
| PetscScalar tmp; | | PetscScalar tmp; |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| double a, b; | | double a, b; |
| if(fscanf(file, "%lf %lf", &a, &b) != 2) | | if(fscanf(file, "%lf %lf", &a, &b) != 2) |
| Message::Error("Could not read data in vector"); | | Message::Error("Could not read data in vector"); |
| tmp = a + PETSC_i * b; | | tmp = a + PETSC_i * b; |
| #else | | #else |
| double a; | | double a; |
| if(fscanf(file, "%lf", &a) != 1) | | if(fscanf(file, "%lf", &a) != 1) |
| Message::Error("Could not read data in vector"); | | Message::Error("Could not read data in vector"); |
| | |
| skipping to change at line 399 | | skipping to change at line 389 |
|---|
| | |
| void LinAlg_ReadScalar(FILE *file, gScalar *S) | | void LinAlg_ReadScalar(FILE *file, gScalar *S) |
| { | | { |
| Message::Error("ReadScalar not yet implemented"); | | Message::Error("ReadScalar not yet implemented"); |
| } | | } |
| | |
| void LinAlg_ReadVector(FILE *file, gVector *V) | | void LinAlg_ReadVector(FILE *file, gVector *V) |
| { | | { |
| PetscInt n; | | PetscInt n; |
| _try(VecGetSize(V->V, &n)); | | _try(VecGetSize(V->V, &n)); |
|
| PetscScalar *tmp = (PetscScalar*)Malloc(n*sizeof(PetscScalar)); | | PetscScalar *tmp = (PetscScalar *)Malloc(n * sizeof(PetscScalar)); |
| if(!fread(tmp, sizeof(PetscScalar), n, file)){ | | if(!fread(tmp, sizeof(PetscScalar), n, file)) { |
| Message::Error("Could not read vector"); | | Message::Error("Could not read vector"); |
| return; | | return; |
| } | | } |
|
| for(PetscInt i = 0; i < n; i++){ | | for(PetscInt i = 0; i < n; i++) { |
| _try(VecSetValues(V->V, 1, &i, &tmp[i], INSERT_VALUES)); | | _try(VecSetValues(V->V, 1, &i, &tmp[i], INSERT_VALUES)); |
| } | | } |
| LinAlg_AssembleVector(V); | | LinAlg_AssembleVector(V); |
| Free(tmp); | | Free(tmp); |
| } | | } |
| | |
| void LinAlg_ReadMatrix(FILE *file, gMatrix *M) | | void LinAlg_ReadMatrix(FILE *file, gMatrix *M) |
| { | | { |
| Message::Error("ReadMatrix not yet implemented"); | | Message::Error("ReadMatrix not yet implemented"); |
| } | | } |
| | |
| skipping to change at line 426 | | skipping to change at line 416 |
|---|
| void LinAlg_PrintScalar(FILE *file, gScalar *S) | | void LinAlg_PrintScalar(FILE *file, gScalar *S) |
| { | | { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| fprintf(file, "%.16g %.16g", real(S->s), imag(S->s)); | | fprintf(file, "%.16g %.16g", real(S->s), imag(S->s)); |
| #else | | #else |
| fprintf(file, "%.16g", S->s); | | fprintf(file, "%.16g", S->s); |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_PrintVector(FILE *file, gVector *V, bool matlab, | | void LinAlg_PrintVector(FILE *file, gVector *V, bool matlab, |
|
| const char* fileName, const char* varName) | | const char *fileName, const char *varName) |
| { | | { |
|
| if(!matlab){ | | if(!matlab) { |
| PetscInt n; | | PetscInt n; |
| _try(VecGetSize(V->V, &n)); | | _try(VecGetSize(V->V, &n)); |
| Vec VV = V->haveSeq ? V->Vseq : V->V; | | Vec VV = V->haveSeq ? V->Vseq : V->V; |
| PetscScalar *tmp; | | PetscScalar *tmp; |
| _try(VecGetArray(VV, &tmp)); | | _try(VecGetArray(VV, &tmp)); |
|
| for (int i = 0; i < n; i++){ | | for(int i = 0; i < n; i++) { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| fprintf(file, "%.16g %.16g\n", real(tmp[i]), imag(tmp[i])); | | fprintf(file, "%.16g %.16g\n", real(tmp[i]), imag(tmp[i])); |
| #else | | #else |
| fprintf(file, "%.16g\n", tmp[i]); | | fprintf(file, "%.16g\n", tmp[i]); |
| #endif | | #endif |
| } | | } |
| fflush(file); | | fflush(file); |
| _try(VecRestoreArray(VV, &tmp)); | | _try(VecRestoreArray(VV, &tmp)); |
| } | | } |
|
| else{ | | else { |
| PetscViewer fd; | | PetscViewer fd; |
| _try(PetscViewerASCIIOpen(MyComm, fileName, &fd)); | | _try(PetscViewerASCIIOpen(MyComm, fileName, &fd)); |
| _try(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB)); | | _try(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB)); |
| _try(PetscObjectSetName((PetscObject)V->V, varName)); | | _try(PetscObjectSetName((PetscObject)V->V, varName)); |
| _try(VecView(V->V, fd)); | | _try(VecView(V->V, fd)); |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) |
| _try(PetscViewerDestroy(&fd)); | | _try(PetscViewerDestroy(&fd)); |
| #else | | #else |
| _try(PetscViewerDestroy(fd)); | | _try(PetscViewerDestroy(fd)); |
| #endif | | #endif |
| } | | } |
| } | | } |
| | |
| void LinAlg_PrintMatrix(FILE *file, gMatrix *M, bool matlab, | | void LinAlg_PrintMatrix(FILE *file, gMatrix *M, bool matlab, |
|
| const char* fileName, const char* varName) | | const char *fileName, const char *varName) |
| { | | { |
|
| if(!matlab){ | | if(!matlab) { |
| PetscInt n, m; | | PetscInt n, m; |
| _try(MatGetSize(M->M, &n, &m)); | | _try(MatGetSize(M->M, &n, &m)); |
|
| for(int i = 0; i < n; i++){ | | for(int i = 0; i < n; i++) { |
| PetscInt ncols; | | PetscInt ncols; |
| const PetscInt *cols; | | const PetscInt *cols; |
| const PetscScalar *vals; | | const PetscScalar *vals; |
| _try(MatGetRow(M->M, i, &ncols, &cols, &vals)); | | _try(MatGetRow(M->M, i, &ncols, &cols, &vals)); |
|
| for(int j = 0; j < m; j++){ | | for(int j = 0; j < m; j++) { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
|
| fprintf(file, "[%d, %d] %.16g %.16g\n", i, j, real(vals[j]), imag(vals[j | | fprintf(file, "[%d, %d] %.16g %.16g\n", i, j, real(vals[j]), |
| ]) ); | | imag(vals[j])); |
| #else | | #else |
| fprintf(file, "[%d, %d] %.16g\n", i, j, vals[j]); | | fprintf(file, "[%d, %d] %.16g\n", i, j, vals[j]); |
| #endif | | #endif |
| } | | } |
| _try(MatRestoreRow(M->M, i, &ncols, &cols, &vals)); | | _try(MatRestoreRow(M->M, i, &ncols, &cols, &vals)); |
| } | | } |
| } | | } |
|
| else{ | | else { |
| // ASCII | | // ASCII |
| PetscViewer fd; | | PetscViewer fd; |
| _try(PetscViewerASCIIOpen(MyComm, fileName, &fd)); | | _try(PetscViewerASCIIOpen(MyComm, fileName, &fd)); |
| _try(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB)); | | _try(PetscViewerSetFormat(fd, PETSC_VIEWER_ASCII_MATLAB)); |
| _try(PetscObjectSetName((PetscObject)M->M, varName)); | | _try(PetscObjectSetName((PetscObject)M->M, varName)); |
| _try(MatView(M->M, fd)); | | _try(MatView(M->M, fd)); |
| | |
| // Binary | | // Binary |
| PetscViewer fd2; | | PetscViewer fd2; |
| std::string tmp(fileName); | | std::string tmp(fileName); |
|
| _try(PetscViewerBinaryOpen(MyComm, (tmp + ".bin").c_str(), FILE_MODE_WRITE, | | _try(PetscViewerBinaryOpen(MyComm, (tmp + ".bin").c_str(), FILE_MODE_WRITE, |
| &fd2)); | | &fd2)); |
| _try(PetscViewerSetFormat(fd2, PETSC_VIEWER_DEFAULT)); | | _try(PetscViewerSetFormat(fd2, PETSC_VIEWER_DEFAULT)); |
| _try(MatView(M->M, fd2)); | | _try(MatView(M->M, fd2)); |
| | |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 2) |
| _try(PetscViewerDestroy(&fd)); | | _try(PetscViewerDestroy(&fd)); |
| _try(PetscViewerDestroy(&fd2)); | | _try(PetscViewerDestroy(&fd2)); |
| #else | | #else |
| _try(PetscViewerDestroy(fd)); | | _try(PetscViewerDestroy(fd)); |
| _try(PetscViewerDestroy(fd2)); | | _try(PetscViewerDestroy(fd2)); |
| #endif | | #endif |
| } | | } |
| } | | } |
| | |
| void LinAlg_WriteScalar(FILE *file, gScalar *S) | | void LinAlg_WriteScalar(FILE *file, gScalar *S) |
| | |
| skipping to change at line 657 | | skipping to change at line 649 |
|---|
| { | | { |
| if(!_isInLocalRange(M, i)) return; | | if(!_isInLocalRange(M, i)) return; |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
| _try(MatGetValues(M->M, 1, &ti, 1, &tj, &S->s)); | | _try(MatGetValues(M->M, 1, &ti, 1, &tj, &S->s)); |
| } | | } |
| | |
| void LinAlg_GetDoubleInMatrix(double *d, gMatrix *M, int i, int j) | | void LinAlg_GetDoubleInMatrix(double *d, gMatrix *M, int i, int j) |
| { | | { |
| if(!_isInLocalRange(M, i)) return; | | if(!_isInLocalRange(M, i)) return; |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| _try(MatGetValues(M->M, 1, &ti, 1, &tj, (PetscScalar*)d)); | | _try(MatGetValues(M->M, 1, &ti, 1, &tj, (PetscScalar *)d)); |
| } | | } |
| | |
|
| void LinAlg_GetComplexInMatrix(double *d1, double *d2, gMatrix *M, | | void LinAlg_GetComplexInMatrix(double *d1, double *d2, gMatrix *M, int i, int j, |
| int i, int j, int k, int l) | | int k, int l) |
| { | | { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| PetscScalar tmp; | | PetscScalar tmp; |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| if(_isInLocalRange(M, i)){ | | if(_isInLocalRange(M, i)) { |
| _try(MatGetValues(M->M, 1, &ti, 1, &tj, &tmp)); | | _try(MatGetValues(M->M, 1, &ti, 1, &tj, &tmp)); |
|
| *d1 = real(tmp) ; | | *d1 = real(tmp); |
| *d2 = imag(tmp) ; | | *d2 = imag(tmp); |
| } | | } |
| #else | | #else |
| PetscInt ti = i, tj = j, tk = k, tl = l; | | PetscInt ti = i, tj = j, tk = k, tl = l; |
| if(_isInLocalRange(M, i)) | | if(_isInLocalRange(M, i)) |
|
| _try(MatGetValues(M->M, 1, &ti, 1, &tj, (PetscScalar*)d1)); | | _try(MatGetValues(M->M, 1, &ti, 1, &tj, (PetscScalar *)d1)); |
| if(_isInLocalRange(M, k)) | | if(_isInLocalRange(M, k)) |
|
| _try(MatGetValues(M->M, 1, &tk, 1, &tl, (PetscScalar*)d2)); | | _try(MatGetValues(M->M, 1, &tk, 1, &tl, (PetscScalar *)d2)); |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_GetColumnInMatrix(gMatrix *M, int col, gVector *V1) | | void LinAlg_GetColumnInMatrix(gMatrix *M, int col, gVector *V1) |
| { | | { |
| Message::Error("GetColumnInMatrix not yet implemented"); | | Message::Error("GetColumnInMatrix not yet implemented"); |
| } | | } |
| | |
| void LinAlg_SetScalar(gScalar *S, double *d) | | void LinAlg_SetScalar(gScalar *S, double *d) |
| { | | { |
| | |
| skipping to change at line 720 | | skipping to change at line 712 |
|---|
| if(!_isInLocalRange(V, i)) return; | | if(!_isInLocalRange(V, i)) return; |
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
| PetscInt ti = i; | | PetscInt ti = i; |
| _try(VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES)); |
| } | | } |
| | |
| void LinAlg_SetComplexInVector(double d1, double d2, gVector *V, int i, int j) | | void LinAlg_SetComplexInVector(double d1, double d2, gVector *V, int i, int j) |
| { | | { |
| PetscScalar tmp; | | PetscScalar tmp; |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
|
| if(_isInLocalRange(V, i)){ | | if(_isInLocalRange(V, i)) { |
| PetscInt ti = i; | | PetscInt ti = i; |
| tmp = d1 + PETSC_i * d2; | | tmp = d1 + PETSC_i * d2; |
| _try(VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES)); |
| } | | } |
| #else | | #else |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| if(_isInLocalRange(V, i)){ | | if(_isInLocalRange(V, i)) { |
| tmp = d1; | | tmp = d1; |
| _try(VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &tmp, INSERT_VALUES)); |
| } | | } |
|
| if(_isInLocalRange(V, j)){ | | if(_isInLocalRange(V, j)) { |
| tmp = d2; | | tmp = d2; |
| _try(VecSetValues(V->V, 1, &tj, &tmp, INSERT_VALUES)); | | _try(VecSetValues(V->V, 1, &tj, &tmp, INSERT_VALUES)); |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_SetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) | | void LinAlg_SetScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) |
| { | | { |
| if(!_isInLocalRange(M, i)) return; | | if(!_isInLocalRange(M, i)) return; |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &S->s, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &S->s, INSERT_VALUES)); |
| } | | } |
| | |
| void LinAlg_SetDoubleInMatrix(double d, gMatrix *M, int i, int j) | | void LinAlg_SetDoubleInMatrix(double d, gMatrix *M, int i, int j) |
| { | | { |
| if(!_isInLocalRange(M, i)) return; | | if(!_isInLocalRange(M, i)) return; |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, (PetscScalar*)&d, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, (PetscScalar *)&d, INSERT_VALUES)); |
| } | | } |
| | |
|
| void LinAlg_SetComplexInMatrix(double d1, double d2, gMatrix *M, | | void LinAlg_SetComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, |
| int i, int j, int k, int l) | | int k, int l) |
| { | | { |
| PetscScalar tmp; | | PetscScalar tmp; |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| if(_isInLocalRange(M, i)){ | | if(_isInLocalRange(M, i)) { |
| tmp = d1 + PETSC_i * d2; | | tmp = d1 + PETSC_i * d2; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, INSERT_VALUES)); |
| } | | } |
| #else | | #else |
| PetscInt ti = i, tj = j, tk = k, tl = l; | | PetscInt ti = i, tj = j, tk = k, tl = l; |
|
| if(d1){ | | if(d1) { |
| tmp = d1; | | tmp = d1; |
| if(_isInLocalRange(M, i)) | | if(_isInLocalRange(M, i)) |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, INSERT_VALUES)); |
| if(_isInLocalRange(M, k)) | | if(_isInLocalRange(M, k)) |
| _try(MatSetValues(M->M, 1, &tk, 1, &tl, &tmp, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &tk, 1, &tl, &tmp, INSERT_VALUES)); |
| } | | } |
|
| if(d2){ | | if(d2) { |
| if(_isInLocalRange(M, i)){ | | if(_isInLocalRange(M, i)) { |
| tmp = -d2; | | tmp = -d2; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tl, &tmp, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tl, &tmp, INSERT_VALUES)); |
| } | | } |
|
| if(_isInLocalRange(M, k)){ | | if(_isInLocalRange(M, k)) { |
| tmp = d2; | | tmp = d2; |
| _try(MatSetValues(M->M, 1, &tk, 1, &tj, &tmp, INSERT_VALUES)); | | _try(MatSetValues(M->M, 1, &tk, 1, &tj, &tmp, INSERT_VALUES)); |
| } | | } |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_AddScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_AddScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
| S3->s = S1->s + S2->s; | | S3->s = S1->s + S2->s; |
| } | | } |
| | |
| void LinAlg_DummyVector(gVector *V) | | void LinAlg_DummyVector(gVector *V) |
| { | | { |
| PetscInt n; | | PetscInt n; |
| PetscScalar zero = 0.0; | | PetscScalar zero = 0.0; |
| | |
|
| if (Current.DofData->DummyDof == NULL) return ; | | if(Current.DofData->DummyDof == NULL) return; |
| | |
| _try(VecGetSize(V->V, &n)); | | _try(VecGetSize(V->V, &n)); |
| for(PetscInt i = 0; i < n; i++) | | for(PetscInt i = 0; i < n; i++) |
|
| if (Current.DofData->DummyDof[i]==1) | | if(Current.DofData->DummyDof[i] == 1) |
| _try(VecSetValues(V->V, 1, &i, &zero, INSERT_VALUES)); | | _try(VecSetValues(V->V, 1, &i, &zero, INSERT_VALUES)); |
| } | | } |
| | |
| void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i) | | void LinAlg_AddScalarInVector(gScalar *S, gVector *V, int i) |
| { | | { |
| if(!_isInLocalRange(V, i)) return; | | if(!_isInLocalRange(V, i)) return; |
| | |
| if(Current.DofData->DummyDof) | | if(Current.DofData->DummyDof) |
|
| if(Current.DofData->DummyDof[i]==1) return ; | | if(Current.DofData->DummyDof[i] == 1) return; |
| | |
| PetscInt ti = i; | | PetscInt ti = i; |
| _try(VecSetValues(V->V, 1, &ti, &S->s, ADD_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &S->s, ADD_VALUES)); |
| } | | } |
| | |
| void LinAlg_AddDoubleInVector(double d, gVector *V, int i) | | void LinAlg_AddDoubleInVector(double d, gVector *V, int i) |
| { | | { |
| if(!_isInLocalRange(V, i)) return; | | if(!_isInLocalRange(V, i)) return; |
| | |
| if(Current.DofData->DummyDof) | | if(Current.DofData->DummyDof) |
|
| if(Current.DofData->DummyDof[i]==1) return ; | | if(Current.DofData->DummyDof[i] == 1) return; |
| | |
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
| PetscInt ti = i; | | PetscInt ti = i; |
| _try(VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES)); |
| } | | } |
| | |
| void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j) | | void LinAlg_AddComplexInVector(double d1, double d2, gVector *V, int i, int j) |
| { | | { |
| PetscScalar tmp; | | PetscScalar tmp; |
|
| int iok=1, jok=1; | | int iok = 1, jok = 1; |
| | |
|
| if(Current.DofData->DummyDof){ | | if(Current.DofData->DummyDof) { |
| if(Current.DofData->DummyDof[i]==1) iok=0; | | if(Current.DofData->DummyDof[i] == 1) iok = 0; |
| if(Current.DofData->DummyDof[j]==1) jok=0; | | if(Current.DofData->DummyDof[j] == 1) jok = 0; |
| } | | } |
| | |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
|
| if(_isInLocalRange(V, i) && iok && jok){ | | if(_isInLocalRange(V, i) && iok && jok) { |
| PetscInt ti = i; | | PetscInt ti = i; |
| tmp = d1 + PETSC_i * d2; | | tmp = d1 + PETSC_i * d2; |
| _try(VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES)); |
| } | | } |
| #else | | #else |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| if(_isInLocalRange(V, i) && iok){ | | if(_isInLocalRange(V, i) && iok) { |
| tmp = d1; | | tmp = d1; |
| _try(VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES)); | | _try(VecSetValues(V->V, 1, &ti, &tmp, ADD_VALUES)); |
| } | | } |
|
| if(_isInLocalRange(V, j) && jok){ | | if(_isInLocalRange(V, j) && jok) { |
| tmp = d2; | | tmp = d2; |
| _try(VecSetValues(V->V, 1, &tj, &tmp, ADD_VALUES)); | | _try(VecSetValues(V->V, 1, &tj, &tmp, ADD_VALUES)); |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_AddScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) | | void LinAlg_AddScalarInMatrix(gScalar *S, gMatrix *M, int i, int j) |
| { | | { |
| if(!_isInLocalRange(M, i)) return; | | if(!_isInLocalRange(M, i)) return; |
| | |
|
| if (Current.DofData->DummyDof) | | if(Current.DofData->DummyDof) |
| if ( (Current.DofData->DummyDof[i]==1 || Current.DofData->DummyDof[j]==1) && | | if((Current.DofData->DummyDof[i] == 1 || |
| (i!=j) ) | | Current.DofData->DummyDof[j] == 1) && |
| | (i != j)) |
| return; | | return; |
| | |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &S->s, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &S->s, ADD_VALUES)); |
| } | | } |
| | |
| void LinAlg_AddDoubleInMatrix(double d, gMatrix *M, int i, int j) | | void LinAlg_AddDoubleInMatrix(double d, gMatrix *M, int i, int j) |
| { | | { |
| if(!_isInLocalRange(M, i)) return; | | if(!_isInLocalRange(M, i)) return; |
| | |
|
| if (Current.DofData->DummyDof) | | if(Current.DofData->DummyDof) |
| if ( (Current.DofData->DummyDof[i]==1 || Current.DofData->DummyDof[j]==1) && | | if((Current.DofData->DummyDof[i] == 1 || |
| (i!=j) ) | | Current.DofData->DummyDof[j] == 1) && |
| | (i != j)) |
| return; | | return; |
| | |
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES)); |
| } | | } |
| | |
|
| void LinAlg_AddComplexInMatrix(double d1, double d2, gMatrix *M, | | void LinAlg_AddComplexInMatrix(double d1, double d2, gMatrix *M, int i, int j, |
| int i, int j, int k, int l) | | int k, int l) |
| { | | { |
| PetscScalar tmp; | | PetscScalar tmp; |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| PetscInt ti = i, tj = j; | | PetscInt ti = i, tj = j; |
|
| if(_isInLocalRange(M, i)){ | | if(_isInLocalRange(M, i)) { |
| tmp = d1 + PETSC_i * d2; | | tmp = d1 + PETSC_i * d2; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES)); |
| } | | } |
| #else | | #else |
| PetscInt ti = i, tj = j, tk = k, tl = l; | | PetscInt ti = i, tj = j, tk = k, tl = l; |
|
| if(d1){ | | if(d1) { |
| tmp = d1; | | tmp = d1; |
| if(_isInLocalRange(M, i)) | | if(_isInLocalRange(M, i)) |
| _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tj, &tmp, ADD_VALUES)); |
| if(_isInLocalRange(M, k)) | | if(_isInLocalRange(M, k)) |
| _try(MatSetValues(M->M, 1, &tk, 1, &tl, &tmp, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &tk, 1, &tl, &tmp, ADD_VALUES)); |
| } | | } |
|
| if(d2){ | | if(d2) { |
| if(_isInLocalRange(M, i)){ | | if(_isInLocalRange(M, i)) { |
| tmp = -d2; | | tmp = -d2; |
| _try(MatSetValues(M->M, 1, &ti, 1, &tl, &tmp, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &ti, 1, &tl, &tmp, ADD_VALUES)); |
| } | | } |
|
| if(_isInLocalRange(M, k)){ | | if(_isInLocalRange(M, k)) { |
| tmp = d2; | | tmp = d2; |
| _try(MatSetValues(M->M, 1, &tk, 1, &tj, &tmp, ADD_VALUES)); | | _try(MatSetValues(M->M, 1, &tk, 1, &tj, &tmp, ADD_VALUES)); |
| } | | } |
| } | | } |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_AddVectorVector(gVector *V1, gVector *V2, gVector *V3) | | void LinAlg_AddVectorVector(gVector *V1, gVector *V2, gVector *V3) |
| { | | { |
| PetscScalar tmp = 1.0; | | PetscScalar tmp = 1.0; |
|
| if(V3 == V1){ | | if(V3 == V1) { |
| _try(VecAXPY(V1->V, tmp, V2->V)); | | _try(VecAXPY(V1->V, tmp, V2->V)); |
| _fillseq(V1); | | _fillseq(V1); |
| } | | } |
|
| else if(V3 == V2){ | | else if(V3 == V2) { |
| _try(VecAXPY(V2->V, tmp, V1->V)); | | _try(VecAXPY(V2->V, tmp, V1->V)); |
| _fillseq(V2); | | _fillseq(V2); |
| } | | } |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddVectorVector'"); | | Message::Error("Wrong arguments in 'LinAlg_AddVectorVector'"); |
| } | | } |
| | |
|
| void LinAlg_AddVectorProdVectorDouble(gVector *V1, gVector *V2, double d, gVecto | | void LinAlg_AddVectorProdVectorDouble(gVector *V1, gVector *V2, double d, |
| r *V3) | | gVector *V3) |
| { | | { |
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
|
| if(V3 == V1){ | | if(V3 == V1) { |
| _try(VecAXPY(V1->V, tmp, V2->V)); | | _try(VecAXPY(V1->V, tmp, V2->V)); |
| _fillseq(V1); | | _fillseq(V1); |
| } | | } |
|
| else if(V3 == V2){ | | else if(V3 == V2) { |
| _try(VecAYPX(V2->V, tmp, V1->V)); | | _try(VecAYPX(V2->V, tmp, V1->V)); |
| _fillseq(V2); | | _fillseq(V2); |
| } | | } |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddVectorProdVectorDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_AddVectorProdVectorDouble'"); |
| } | | } |
| | |
| void LinAlg_AddProdVectorDoubleProdVectorDouble(double alpha, gVector *V1, | | void LinAlg_AddProdVectorDoubleProdVectorDouble(double alpha, gVector *V1, |
|
| double beta, gVector *V2, gVector *V3) | | double beta, gVector *V2, |
| | gVector *V3) |
| { | | { |
| PetscScalar alpha1 = alpha, beta1 = beta; | | PetscScalar alpha1 = alpha, beta1 = beta; |
| PetscScalar gamma1 = 0.0; | | PetscScalar gamma1 = 0.0; |
| _try(VecAXPBYPCZ(V3->V, alpha1, beta1, gamma1, V1->V, V2->V)); | | _try(VecAXPBYPCZ(V3->V, alpha1, beta1, gamma1, V1->V, V2->V)); |
| } | | } |
| | |
| void LinAlg_AddMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) | | void LinAlg_AddMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) |
| { | | { |
| PetscScalar tmp = 1.0; | | PetscScalar tmp = 1.0; |
| if(M3 == M1) | | if(M3 == M1) |
| _try(MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN)); |
| else if(M3 == M2) | | else if(M3 == M2) |
| _try(MatAXPY(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatAXPY(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN)); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddMatrixMatrix'"); | | Message::Error("Wrong arguments in 'LinAlg_AddMatrixMatrix'"); |
| } | | } |
| | |
|
| void LinAlg_AddMatrixProdMatrixDouble(gMatrix *M1, gMatrix *M2, double d, gMatri | | void LinAlg_AddMatrixProdMatrixDouble(gMatrix *M1, gMatrix *M2, double d, |
| x *M3) | | gMatrix *M3) |
| { | | { |
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
| if(M3 == M1) | | if(M3 == M1) |
| _try(MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN)); |
| else if(M3 == M2) | | else if(M3 == M2) |
|
| #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR == 3) && (PETSC_VERSION_S | | #if(PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR == 3) && \ |
| UBMINOR < 2) | | (PETSC_VERSION_SUBMINOR < 2) |
| _try(MatAYPX(M2->M, tmp, M1->M)); | | _try(MatAYPX(M2->M, tmp, M1->M)); |
| #else | | #else |
| _try(MatAYPX(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatAYPX(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN)); |
| #endif | | #endif |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_AddMatrixProdMatrixDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_AddMatrixProdMatrixDouble'"); |
| } | | } |
| | |
| void LinAlg_SubScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_SubScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
| S3->s = S1->s - S2->s; | | S3->s = S1->s - S2->s; |
| } | | } |
| | |
| void LinAlg_SubVectorVector(gVector *V1, gVector *V2, gVector *V3) | | void LinAlg_SubVectorVector(gVector *V1, gVector *V2, gVector *V3) |
| { | | { |
| PetscScalar tmp = -1.0; | | PetscScalar tmp = -1.0; |
|
| if(V3 == V1){ | | if(V3 == V1) { |
| _try(VecAXPY(V1->V, tmp, V2->V)); // V1->V = V1->V - V2->V | | _try(VecAXPY(V1->V, tmp, V2->V)); // V1->V = V1->V - V2->V |
| _fillseq(V1); | | _fillseq(V1); |
| } | | } |
|
| else if(V3 == V2){ | | else if(V3 == V2) { |
| _try(VecAYPX(V2->V, tmp, V1->V)); // V2->V = V1->V - V2->V | | _try(VecAYPX(V2->V, tmp, V1->V)); // V2->V = V1->V - V2->V |
| _fillseq(V2); | | _fillseq(V2); |
| } | | } |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_SubVectorVector'"); | | Message::Error("Wrong arguments in 'LinAlg_SubVectorVector'"); |
| } | | } |
| | |
| void LinAlg_SubMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) | | void LinAlg_SubMatrixMatrix(gMatrix *M1, gMatrix *M2, gMatrix *M3) |
| { | | { |
|
| | | |
| PetscScalar tmp = -1.0; | | PetscScalar tmp = -1.0; |
| if(M3 == M1) // M1->M = M1->M - M2->M | | if(M3 == M1) // M1->M = M1->M - M2->M |
| _try(MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatAXPY(M1->M, tmp, M2->M, DIFFERENT_NONZERO_PATTERN)); |
| else if(M3 == M2) // M2->M = M1->M - M2->M | | else if(M3 == M2) // M2->M = M1->M - M2->M |
| _try(MatAYPX(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN)); | | _try(MatAYPX(M2->M, tmp, M1->M, DIFFERENT_NONZERO_PATTERN)); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_SubMatrixMatrix'"); | | Message::Error("Wrong arguments in 'LinAlg_SubMatrixMatrix'"); |
| } | | } |
| | |
| void LinAlg_ProdScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) | | void LinAlg_ProdScalarScalar(gScalar *S1, gScalar *S2, gScalar *S3) |
| { | | { |
| S3->s = S1->s * S2->s; | | S3->s = S1->s * S2->s; |
| } | | } |
| | |
| void LinAlg_ProdScalarDouble(gScalar *S1, double d, gScalar *S2) | | void LinAlg_ProdScalarDouble(gScalar *S1, double d, gScalar *S2) |
| { | | { |
| S2->s = S1->s * d; | | S2->s = S1->s * d; |
| } | | } |
| | |
|
| void LinAlg_ProdScalarComplex(gScalar *S, double d1, double d2, double *d3, doub | | void LinAlg_ProdScalarComplex(gScalar *S, double d1, double d2, double *d3, |
| le *d4) | | double *d4) |
| { | | { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| PetscScalar tmp; | | PetscScalar tmp; |
| #endif | | #endif |
| | |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| tmp = S->s * (d1 + PETSC_i * d2); | | tmp = S->s * (d1 + PETSC_i * d2); |
| *d3 = real(tmp); | | *d3 = real(tmp); |
| *d4 = imag(tmp); | | *d4 = imag(tmp); |
| #else | | #else |
| *d3 = S->s * d1; | | *d3 = S->s * d1; |
| *d4 = S->s * d2; | | *d4 = S->s * d2; |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_ProdVectorScalar(gVector *V1, gScalar *S, gVector *V2) | | void LinAlg_ProdVectorScalar(gVector *V1, gScalar *S, gVector *V2) |
| { | | { |
|
| if(V2 == V1){ | | if(V2 == V1) { |
| _try(VecScale(V1->V, S->s)); | | _try(VecScale(V1->V, S->s)); |
| _fillseq(V1); | | _fillseq(V1); |
| } | | } |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_ProdVectorScalar'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdVectorScalar'"); |
| } | | } |
| | |
| void LinAlg_ProdVectorDouble(gVector *V1, double d, gVector *V2) | | void LinAlg_ProdVectorDouble(gVector *V1, double d, gVector *V2) |
| { | | { |
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
|
| if(V2 == V1){ | | if(V2 == V1) { |
| _try(VecScale(V1->V, tmp)); | | _try(VecScale(V1->V, tmp)); |
| _fillseq(V1); | | _fillseq(V1); |
| } | | } |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_ProdVectorDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdVectorDouble'"); |
| } | | } |
| | |
| void LinAlg_ProdVectorComplex(gVector *V1, double d1, double d2, gVector *V2) | | void LinAlg_ProdVectorComplex(gVector *V1, double d1, double d2, gVector *V2) |
| { | | { |
| Message::Error("ProdVectorComplex not yet implemented"); | | Message::Error("ProdVectorComplex not yet implemented"); |
| | |
| skipping to change at line 1074 | | skipping to change at line 1074 |
|---|
| *d = real(tmp); | | *d = real(tmp); |
| #else | | #else |
| *d = tmp; | | *d = tmp; |
| #endif | | #endif |
| } | | } |
| | |
| void LinAlg_ProdMatrixVector(gMatrix *M, gVector *V1, gVector *V2) | | void LinAlg_ProdMatrixVector(gMatrix *M, gVector *V1, gVector *V2) |
| { | | { |
| if(V2 == V1) | | if(V2 == V1) |
| Message::Error("Wrong arguments in 'LinAlg_ProdMatrixVector'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdMatrixVector'"); |
|
| else{ | | else { |
| _try(MatMult(M->M, V1->V, V2->V)); | | _try(MatMult(M->M, V1->V, V2->V)); |
| _fillseq(V2); | | _fillseq(V2); |
| } | | } |
| } | | } |
| | |
| void LinAlg_ProdMatrixScalar(gMatrix *M1, gScalar *S, gMatrix *M2) | | void LinAlg_ProdMatrixScalar(gMatrix *M1, gScalar *S, gMatrix *M2) |
| { | | { |
| if(M2 == M1) | | if(M2 == M1) |
| _try(MatScale(M1->M, S->s)); | | _try(MatScale(M1->M, S->s)); |
| else | | else |
| | |
| skipping to change at line 1100 | | skipping to change at line 1100 |
|---|
| PetscScalar tmp = d; | | PetscScalar tmp = d; |
| if(M2 == M1) | | if(M2 == M1) |
| _try(MatScale(M1->M, tmp)); | | _try(MatScale(M1->M, tmp)); |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_ProdMatrixDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdMatrixDouble'"); |
| } | | } |
| | |
| void LinAlg_ProdMatrixComplex(gMatrix *M1, double d1, double d2, gMatrix *M2) | | void LinAlg_ProdMatrixComplex(gMatrix *M1, double d1, double d2, gMatrix *M2) |
| { | | { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
|
| if(M2 == M1){ | | if(M2 == M1) { |
| PetscScalar tmp = d1 + (PETSC_i * d2); | | PetscScalar tmp = d1 + (PETSC_i * d2); |
| _try(MatScale(M1->M, tmp)); | | _try(MatScale(M1->M, tmp)); |
| } | | } |
| else | | else |
| Message::Error("Wrong arguments in 'LinAlg_ProdMatrixDouble'"); | | Message::Error("Wrong arguments in 'LinAlg_ProdMatrixDouble'"); |
| #else | | #else |
| Message::Error("ProdMatrixComplex not yet implemented"); | | Message::Error("ProdMatrixComplex not yet implemented"); |
| #endif | | #endif |
| } | | } |
| | |
| | |
| skipping to change at line 1153 | | skipping to change at line 1153 |
|---|
| { | | { |
| Message::Barrier(); | | Message::Barrier(); |
| _try(VecAssemblyBegin(V->V)); | | _try(VecAssemblyBegin(V->V)); |
| _try(VecAssemblyEnd(V->V)); | | _try(VecAssemblyEnd(V->V)); |
| _fillseq(V); | | _fillseq(V); |
| } | | } |
| | |
| #if defined(HAVE_ZITSOL) | | #if defined(HAVE_ZITSOL) |
| | |
| extern "C" { | | extern "C" { |
|
| int getdp_zitsol(int n, int nnz, int *row, int *col, double *valr, double *val | | int getdp_zitsol(int n, int nnz, int *row, int *col, double *valr, double *vali, |
| i, | | double *rhsr, double *rhsi, double *solr, double *soli, |
| double *rhsr, double *rhsi, double *solr, double *soli, | | int precond, int lfil, double tol0, double tol, int im, |
| int precond, int lfil, double tol0, double tol, int im, int ma | | int maxits); |
| xits); | | |
| } | | } |
| | |
| static void _zitsol(gMatrix *A, gVector *B, gVector *X) | | static void _zitsol(gMatrix *A, gVector *B, gVector *X) |
| { | | { |
| int precond = 1, lfil = 30, im = 100, maxits = 200; | | int precond = 1, lfil = 30, im = 100, maxits = 200; |
| double tol0 = 0.01, tol = 1e-10; | | double tol0 = 0.01, tol = 1e-10; |
| PetscTruth set; | | PetscTruth set; |
| PetscOptionsGetInt(PETSC_NULL, "-zitsol_precond", &precond, &set); | | PetscOptionsGetInt(PETSC_NULL, "-zitsol_precond", &precond, &set); |
| PetscOptionsGetInt(PETSC_NULL, "-zitsol_lfil", &lfil, &set); | | PetscOptionsGetInt(PETSC_NULL, "-zitsol_lfil", &lfil, &set); |
| PetscOptionsGetInt(PETSC_NULL, "-zitsol_im", &im, &set); | | PetscOptionsGetInt(PETSC_NULL, "-zitsol_im", &im, &set); |
| PetscOptionsGetInt(PETSC_NULL, "-zitsol_maxits", &maxits, &set); | | PetscOptionsGetInt(PETSC_NULL, "-zitsol_maxits", &maxits, &set); |
| PetscOptionsGetReal(PETSC_NULL, "-zitsol_tol0", &tol0, &set); | | PetscOptionsGetReal(PETSC_NULL, "-zitsol_tol0", &tol0, &set); |
| PetscOptionsGetReal(PETSC_NULL, "-zitsol_tol", &tol, &set); | | PetscOptionsGetReal(PETSC_NULL, "-zitsol_tol", &tol, &set); |
| | |
| MatInfo info; | | MatInfo info; |
| _try(MatGetInfo(A->M, MAT_LOCAL, &info)); | | _try(MatGetInfo(A->M, MAT_LOCAL, &info)); |
| int nnz = info.nz_used; | | int nnz = info.nz_used; |
|
| //int n = info.rows_local; | | // int n = info.rows_local; |
| PetscInt n; | | PetscInt n; |
| _try(VecGetLocalSize(B->V, &n)); | | _try(VecGetLocalSize(B->V, &n)); |
| | |
| Current.KSPSystemSize = n; | | Current.KSPSystemSize = n; |
| | |
|
| int *row = (int*)Malloc(nnz * sizeof(int)); | | int *row = (int *)Malloc(nnz * sizeof(int)); |
| int *col = (int*)Malloc(nnz * sizeof(int)); | | int *col = (int *)Malloc(nnz * sizeof(int)); |
| double *valr = (double*)Malloc(nnz * sizeof(double)); | | double *valr = (double *)Malloc(nnz * sizeof(double)); |
| double *vali = (double*)Malloc(nnz * sizeof(double)); | | double *vali = (double *)Malloc(nnz * sizeof(double)); |
| double *rhsr = (double*)Malloc(n * sizeof(double)); | | double *rhsr = (double *)Malloc(n * sizeof(double)); |
| double *rhsi = (double*)Malloc(n * sizeof(double)); | | double *rhsi = (double *)Malloc(n * sizeof(double)); |
| double *solr = (double*)Malloc(n * sizeof(double)); | | double *solr = (double *)Malloc(n * sizeof(double)); |
| double *soli = (double*)Malloc(n * sizeof(double)); | | double *soli = (double *)Malloc(n * sizeof(double)); |
| | |
| int k = 0; | | int k = 0; |
|
| for(int i = 0; i < n; i++){ | | for(int i = 0; i < n; i++) { |
| PetscInt ncols; | | PetscInt ncols; |
| const PetscInt *cols; | | const PetscInt *cols; |
| const PetscScalar *vals; | | const PetscScalar *vals; |
| _try(MatGetRow(A->M, i, &ncols, &cols, &vals)); | | _try(MatGetRow(A->M, i, &ncols, &cols, &vals)); |
|
| for(int j = 0; j < ncols; j++){ | | for(int j = 0; j < ncols; j++) { |
| if(k >= nnz){ | | if(k >= nnz) { |
| Message::Error("Something wrong in nnz: %d >= %d", k, nnz); | | Message::Error("Something wrong in nnz: %d >= %d", k, nnz); |
| return; | | return; |
| } | | } |
| row[k] = i; | | row[k] = i; |
| col[k] = cols[j]; | | col[k] = cols[j]; |
| Message::Debug("A[%d][%d] = ", row[k], col[k]); | | Message::Debug("A[%d][%d] = ", row[k], col[k]); |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| valr[k] = real(vals[j]); | | valr[k] = real(vals[j]); |
| vali[k] = imag(vals[j]); | | vali[k] = imag(vals[j]); |
| Message::Debug("%g+i*%g", valr[k], vali[k]); | | Message::Debug("%g+i*%g", valr[k], vali[k]); |
| #else | | #else |
| valr[k] = vals[j]; | | valr[k] = vals[j]; |
| | |
| skipping to change at line 1221 | | skipping to change at line 1222 |
|---|
| k++; | | k++; |
| } | | } |
| _try(MatRestoreRow(A->M, i, &ncols, &cols, &vals)); | | _try(MatRestoreRow(A->M, i, &ncols, &cols, &vals)); |
| } | | } |
| | |
| Message::Info("n = %d, nnz = %d (check k = %d)", n, nnz, k); | | Message::Info("n = %d, nnz = %d (check k = %d)", n, nnz, k); |
| | |
| PetscScalar *b, *x; | | PetscScalar *b, *x; |
| _try(VecGetArray(B->V, &b)); | | _try(VecGetArray(B->V, &b)); |
| _try(VecGetArray(X->V, &x)); | | _try(VecGetArray(X->V, &x)); |
|
| for(int i = 0; i < n; i++){ | | for(int i = 0; i < n; i++) { |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| rhsr[i] = real(b[i]); | | rhsr[i] = real(b[i]); |
| rhsi[i] = imag(b[i]); | | rhsi[i] = imag(b[i]); |
| solr[i] = real(x[i]); | | solr[i] = real(x[i]); |
| soli[i] = imag(x[i]); | | soli[i] = imag(x[i]); |
| #else | | #else |
| rhsr[i] = b[i]; | | rhsr[i] = b[i]; |
| rhsi[i] = 0.; | | rhsi[i] = 0.; |
| solr[i] = x[i]; | | solr[i] = x[i]; |
| soli[i] = 0.; | | soli[i] = 0.; |
| | |
| skipping to change at line 1246 | | skipping to change at line 1247 |
|---|
| | |
| int its = getdp_zitsol(n, nnz, row, col, valr, vali, rhsr, rhsi, solr, soli, | | int its = getdp_zitsol(n, nnz, row, col, valr, vali, rhsr, rhsi, solr, soli, |
| precond, lfil, tol0, tol, im, maxits); | | precond, lfil, tol0, tol, im, maxits); |
| if(its >= maxits) | | if(its >= maxits) |
| Message::Error("Did not converge in %d iterations", maxits); | | Message::Error("Did not converge in %d iterations", maxits); |
| else | | else |
| Message::Info("Converged in %d iterations", its); | | Message::Info("Converged in %d iterations", its); |
| | |
| Current.KSPIterations = its; | | Current.KSPIterations = its; |
| | |
|
| for(PetscInt i = 0; i < n; i++){ | | for(PetscInt i = 0; i < n; i++) { |
| PetscScalar d; | | PetscScalar d; |
| #if defined(PETSC_USE_COMPLEX) | | #if defined(PETSC_USE_COMPLEX) |
| d = solr[i] + PETSC_i * soli[i]; | | d = solr[i] + PETSC_i * soli[i]; |
| #else | | #else |
| d = solr[i]; | | d = solr[i]; |
| #endif | | #endif |
| _try(VecSetValues(X->V, 1, &i, &d, INSERT_VALUES)); | | _try(VecSetValues(X->V, 1, &i, &d, INSERT_VALUES)); |
| } | | } |
| | |
|
| Free(row); Free(col); | | Free(row); |
| Free(valr); Free(vali); | | Free(col); |
| Free(rhsr); Free(rhsi); | | Free(valr); |
| Free(solr); Free(soli); | | Free(vali); |
| | Free(rhsr); |
| | Free(rhsi); |
| | Free(solr); |
| | Free(soli); |
| } | | } |
| | |
| #endif | | #endif |
| | |
|
| static PetscErrorCode _myKspMonitor(KSP ksp, PetscInt it, PetscReal rnorm, void | | static PetscErrorCode _myKspMonitor(KSP ksp, PetscInt it, PetscReal rnorm, |
| *mctx) | | void *mctx) |
| { | | { |
| Message::Info("%3ld KSP Residual norm %14.12e", (long)it, rnorm); | | Message::Info("%3ld KSP Residual norm %14.12e", (long)it, rnorm); |
| Current.KSPIteration = it; | | Current.KSPIteration = it; |
| Current.KSPResidual = rnorm; | | Current.KSPResidual = rnorm; |
| return 0; | | return 0; |
| } | | } |
| | |
| static void _solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, | | static void _solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, |
| int precond, int kspIndex) | | int precond, int kspIndex) |
| { | | { |
| #if defined(HAVE_ZITSOL) | | #if defined(HAVE_ZITSOL) |
| // testing Yousef's new preconditioners and solvers | | // testing Yousef's new preconditioners and solvers |
| PetscTruth set, zitsol = PETSC_FALSE; | | PetscTruth set, zitsol = PETSC_FALSE; |
| PetscOptionsGetTruth(PETSC_NULL, "-zitsol", &zitsol, &set); | | PetscOptionsGetTruth(PETSC_NULL, "-zitsol", &zitsol, &set); |
|
| if(zitsol){ _zitsol(A, B, X); return; } | | if(zitsol) { |
| | _zitsol(A, B, X); |
| | return; |
| | } |
| #endif | | #endif |
|
| if(kspIndex < 0 || kspIndex > 9){ | | if(kspIndex < 0 || kspIndex > 9) { |
| Message::Error("Linear Solver index out of range (%d)", kspIndex); | | Message::Error("Linear Solver index out of range (%d)", kspIndex); |
| return; | | return; |
| } | | } |
| | |
| PetscInt i, j; | | PetscInt i, j; |
| _try(MatGetSize(A->M, &i, &j)); | | _try(MatGetSize(A->M, &i, &j)); |
| Current.KSPSystemSize = i; | | Current.KSPSystemSize = i; |
|
| if(!i){ | | if(!i) { |
| Message::Warning("Zero-size system: skipping solve!"); | | Message::Warning("Zero-size system: skipping solve!"); |
| return; | | return; |
| } | | } |
| | |
| int view = !Solver->ksp[kspIndex]; | | int view = !Solver->ksp[kspIndex]; |
| | |
|
| if(kspIndex != 0) | | if(kspIndex != 0) Message::Info("Using solver index %d", kspIndex); |
| Message::Info("Using solver index %d", kspIndex); | | |
| | |
| if(!Solver->ksp[kspIndex]) { | | if(!Solver->ksp[kspIndex]) { |
| _try(KSPCreate(MyComm, &Solver->ksp[kspIndex])); | | _try(KSPCreate(MyComm, &Solver->ksp[kspIndex])); |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5) |
| _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M)); | | _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M)); |
| #else | | #else |
|
| _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M, DIFFERENT_NONZERO_PA | | _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M, |
| TTERN)); | | DIFFERENT_NONZERO_PATTERN)); |
| #endif | | #endif |
|
| _try(KSPMonitorSet(Solver->ksp[kspIndex], _myKspMonitor, PETSC_NULL, PETSC_N | | _try(KSPMonitorSet(Solver->ksp[kspIndex], _myKspMonitor, PETSC_NULL, |
| ULL)); | | PETSC_NULL)); |
| PC pc; | | PC pc; |
| _try(KSPGetPC(Solver->ksp[kspIndex], &pc)); | | _try(KSPGetPC(Solver->ksp[kspIndex], &pc)); |
| | |
| // set some default options: use direct solver (PARDISO, MUMPS, UMFPACK, or | | // set some default options: use direct solver (PARDISO, MUMPS, UMFPACK, or |
| // native PETSc LU) | | // native PETSc LU) |
| _try(KSPSetType(Solver->ksp[kspIndex], "preonly")); | | _try(KSPSetType(Solver->ksp[kspIndex], "preonly")); |
| _try(PCSetType(pc, PCLU)); | | _try(PCSetType(pc, PCLU)); |
|
| #if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) | | #if(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) |
| _try(PCFactorSetMatSolverPackage(pc, "mumps")); | | _try(PCFactorSetMatSolverPackage(pc, "mumps")); |
|
| #elif (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MKL_PARDISO) | | #elif(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MKL_PARDISO) |
| _try(PCFactorSetMatSolverPackage(pc, "mkl_pardiso")); | | _try(PCFactorSetMatSolverPackage(pc, "mkl_pardiso")); |
|
| #elif (PETSC_VERSION_MAJOR > 2) && (defined(PETSC_HAVE_UMFPACK) || defined(PETSC | | #elif(PETSC_VERSION_MAJOR > 2) && \ |
| _HAVE_SUITESPARSE)) | | (defined(PETSC_HAVE_UMFPACK) || defined(PETSC_HAVE_SUITESPARSE)) |
| _try(PCFactorSetMatSolverPackage(pc, "umfpack")); | | _try(PCFactorSetMatSolverPackage(pc, "umfpack")); |
| #else | | #else |
| _try(PetscOptionsSetValue("-pc_factor_nonzeros_along_diagonal", "1e-12")); | | _try(PetscOptionsSetValue("-pc_factor_nonzeros_along_diagonal", "1e-12")); |
|
| #if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR == 3) && (PETSC_VERSION_S | | #if(PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR == 3) && \ |
| UBMINOR < 3) | | (PETSC_VERSION_SUBMINOR < 3) |
| _try(PCFactorSetMatOrdering(pc, MATORDERING_RCM)); | | _try(PCFactorSetMatOrdering(pc, MATORDERING_RCM)); |
| #else | | #else |
| _try(PCFactorSetMatOrderingType(pc, MATORDERINGRCM)); | | _try(PCFactorSetMatOrderingType(pc, MATORDERINGRCM)); |
| #endif | | #endif |
| #endif | | #endif |
| | |
| // override the default options with the ones from the option database (if | | // override the default options with the ones from the option database (if |
| // any) | | // any) |
| _try(KSPSetFromOptions(Solver->ksp[kspIndex])); | | _try(KSPSetFromOptions(Solver->ksp[kspIndex])); |
| | |
|
| if(view && (!Message::GetCommRank() || !Message::GetIsCommWorld())){ | | if(view && (!Message::GetCommRank() || !Message::GetIsCommWorld())) { |
| // either we are on parallel (!GetIsCommWorld) or in sequential with rank | | // either we are on parallel (!GetIsCommWorld) or in sequential with rank |
| // = 0 (GetIsCommWorld) | | // = 0 (GetIsCommWorld) |
| | |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 4) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 4) |
| const char *ksptype = ""; | | const char *ksptype = ""; |
| _try(KSPGetType(Solver->ksp[kspIndex], &ksptype)); | | _try(KSPGetType(Solver->ksp[kspIndex], &ksptype)); |
| const char *pctype = ""; | | const char *pctype = ""; |
| _try(PCGetType(pc, &pctype)); | | _try(PCGetType(pc, &pctype)); |
| #else | | #else |
| const KSPType ksptype; | | const KSPType ksptype; |
| _try(KSPGetType(Solver->ksp[kspIndex], &ksptype)); | | _try(KSPGetType(Solver->ksp[kspIndex], &ksptype)); |
| const PCType pctype; | | const PCType pctype; |
| _try(PCGetType(pc, &pctype)); | | _try(PCGetType(pc, &pctype)); |
| #endif | | #endif |
| | |
|
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 9) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 9) |
| MatSolverType stype; | | MatSolverType stype; |
| _try(PCFactorGetMatSolverType(pc, &stype)); | | _try(PCFactorGetMatSolverType(pc, &stype)); |
|
| #elif (PETSC_VERSION_MAJOR > 2) | | #elif(PETSC_VERSION_MAJOR > 2) |
| const MatSolverPackage stype; | | const MatSolverPackage stype; |
| _try(PCFactorGetMatSolverPackage(pc, &stype)); | | _try(PCFactorGetMatSolverPackage(pc, &stype)); |
| #else | | #else |
| const char *stype = ""; | | const char *stype = ""; |
| #endif | | #endif |
| Message::Info("N: %ld - %s %s %s", long(i), ksptype, pctype, stype); | | Message::Info("N: %ld - %s %s %s", long(i), ksptype, pctype, stype); |
| } | | } |
| } | | } |
|
| else if(precond){ | | else if(precond) { |
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5) |
| _try(KSPSetReusePreconditioner(Solver->ksp[kspIndex], PETSC_FALSE)); | | _try(KSPSetReusePreconditioner(Solver->ksp[kspIndex], PETSC_FALSE)); |
| _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M)); | | _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M)); |
| #else | | #else |
|
| _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M, DIFFERENT_NONZERO_PA | | _try(KSPSetOperators(Solver->ksp[kspIndex], A->M, A->M, |
| TTERN)); | | DIFFERENT_NONZERO_PATTERN)); |
| #endif | | #endif |
| } | | } |
|
| else{ | | else { |
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5) | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 5) |
| _try(KSPSetReusePreconditioner(Solver->ksp[kspIndex], PETSC_TRUE)); | | _try(KSPSetReusePreconditioner(Solver->ksp[kspIndex], PETSC_TRUE)); |
| #endif | | #endif |
| } | | } |
| | |
| _try(KSPSolve(Solver->ksp[kspIndex], B->V, X->V)); | | _try(KSPSolve(Solver->ksp[kspIndex], B->V, X->V)); |
| | |
| // copy result on all procs | | // copy result on all procs |
| _fillseq(X); | | _fillseq(X); |
| | |
| if(view && Message::GetVerbosity() > 5) | | if(view && Message::GetVerbosity() > 5) |
| _try(KSPView(Solver->ksp[kspIndex], MyPetscViewer)); | | _try(KSPView(Solver->ksp[kspIndex], MyPetscViewer)); |
| | |
| PetscInt its; | | PetscInt its; |
| _try(KSPGetIterationNumber(Solver->ksp[kspIndex], &its)); | | _try(KSPGetIterationNumber(Solver->ksp[kspIndex], &its)); |
|
| if(!Message::GetCommRank() || !Message::GetIsCommWorld()){ | | if(!Message::GetCommRank() || !Message::GetIsCommWorld()) { |
| if(its > 1) Message::Info("%d iterations", its); | | if(its > 1) Message::Info("%d iterations", its); |
| } | | } |
| Current.KSPIterations = its; | | Current.KSPIterations = its; |
| | |
| PetscTruth set, kspfree = PETSC_FALSE; | | PetscTruth set, kspfree = PETSC_FALSE; |
| PetscOptionsGetTruth(PETSC_NULL, "-kspfree", &kspfree, &set); | | PetscOptionsGetTruth(PETSC_NULL, "-kspfree", &kspfree, &set); |
|
| if(kspfree){ | | if(kspfree) { |
| Message::Info("Freeing KSP solver"); | | Message::Info("Freeing KSP solver"); |
| LinAlg_DestroySolver(Solver); | | LinAlg_DestroySolver(Solver); |
| } | | } |
|
| | | |
| } | | } |
| | |
| void LinAlg_Solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, | | void LinAlg_Solve(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, |
| int solverIndex) | | int solverIndex) |
| { | | { |
| _solve(A, B, Solver, X, 1, solverIndex); | | _solve(A, B, Solver, X, 1, solverIndex); |
| } | | } |
| | |
| void LinAlg_SolveAgain(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, | | void LinAlg_SolveAgain(gMatrix *A, gVector *B, gSolver *Solver, gVector *X, |
| int solverIndex) | | int solverIndex) |
| { | | { |
| _solve(A, B, Solver, X, 0, solverIndex); | | _solve(A, B, Solver, X, 0, solverIndex); |
| } | | } |
| | |
| void LinAlg_SetGlobalSolverOptions(const std::string &opt) | | void LinAlg_SetGlobalSolverOptions(const std::string &opt) |
| { | | { |
| _try(PetscOptionsInsertString(opt.c_str())); | | _try(PetscOptionsInsertString(opt.c_str())); |
| } | | } |
| | |
|
| extern void Generate_Residual (gVector *x, gVector *f) ; | | extern void Generate_Residual(gVector *x, gVector *f); |
| extern void Generate_FullJacobian (gVector *x, gMatrix *Jac) ; | | extern void Generate_FullJacobian(gVector *x, gMatrix *Jac); |
| | |
| static PetscErrorCode _NLFormFunction(SNES snes, Vec x, Vec f, void *mctx) | | static PetscErrorCode _NLFormFunction(SNES snes, Vec x, Vec f, void *mctx) |
| { | | { |
|
| gVector gx, gf ; | | gVector gx, gf; |
| gx.V = x ; | | gx.V = x; |
| gx.haveSeq = 0; | | gx.haveSeq = 0; |
|
| gf.V = f ; | | gf.V = f; |
| gf.haveSeq = 0; | | gf.haveSeq = 0; |
|
| Generate_Residual(&gx, &gf) ; | | Generate_Residual(&gx, &gf); |
| | |
|
| PetscScalar *ff ; | | PetscScalar *ff; |
| _try(VecGetArray(gf.V, &ff)) ; | | _try(VecGetArray(gf.V, &ff)); |
| PetscInt n; | | PetscInt n; |
|
| _try(VecGetSize(f, &n)) ; | | _try(VecGetSize(f, &n)); |
| for(PetscInt i = 0; i < n; i++) | | for(PetscInt i = 0; i < n; i++) |
| _try(VecSetValues(f, 1, &i, &ff[i], INSERT_VALUES)); | | _try(VecSetValues(f, 1, &i, &ff[i], INSERT_VALUES)); |
| _try(VecGetArray(f, &ff)); | | _try(VecGetArray(f, &ff)); |
| return 0; | | return 0; |
| } | | } |
| | |
|
| #if (PETSC_VERSION_MAJOR == 2) || ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_ | | #if(PETSC_VERSION_MAJOR == 2) || \ |
| MINOR < 5)) | | ((PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR < 5)) |
| static PetscErrorCode _NLFormJacobian(SNES snes, Vec x, Mat *J, Mat *PC, | | static PetscErrorCode _NLFormJacobian(SNES snes, Vec x, Mat *J, Mat *PC, |
|
| MatStructure *flag, | | MatStructure *flag, void *mctx) |
| void *mctx) | | |
| { | | { |
|
| gVector gx ; | | gVector gx; |
| gx.V = x ; | | gx.V = x; |
| gx.haveSeq = 0; | | gx.haveSeq = 0; |
|
| gMatrix gJ ; | | gMatrix gJ; |
| gJ.M = *J ; | | gJ.M = *J; |
| Generate_FullJacobian(&gx, &gJ); | | Generate_FullJacobian(&gx, &gJ); |
|
| *J = gJ.M ; | | *J = gJ.M; |
| *flag = DIFFERENT_NONZERO_PATTERN ; | | *flag = DIFFERENT_NONZERO_PATTERN; |
| Message::Barrier(); | | Message::Barrier(); |
| _try(MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY)); | | _try(MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY)); |
| _try(MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY)); | | _try(MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY)); |
|
| if (*PC != *J){ | | if(*PC != *J) { |
| _try(MatAssemblyBegin(*PC, MAT_FINAL_ASSEMBLY)); | | _try(MatAssemblyBegin(*PC, MAT_FINAL_ASSEMBLY)); |
| _try(MatAssemblyEnd(*PC, MAT_FINAL_ASSEMBLY)); | | _try(MatAssemblyEnd(*PC, MAT_FINAL_ASSEMBLY)); |
| } | | } |
| return 0; | | return 0; |
| } | | } |
| #else | | #else |
| static PetscErrorCode _NLFormJacobian(SNES snes, Vec x, Mat J, Mat PC, | | static PetscErrorCode _NLFormJacobian(SNES snes, Vec x, Mat J, Mat PC, |
| void *mctx) | | void *mctx) |
| { | | { |
|
| gVector gx ; | | gVector gx; |
| gx.V = x ; | | gx.V = x; |
| gx.haveSeq = 0; | | gx.haveSeq = 0; |
|
| gMatrix gJ ; | | gMatrix gJ; |
| Generate_FullJacobian(&gx, &gJ); | | Generate_FullJacobian(&gx, &gJ); |
|
| //J = gJ.M; | | // J = gJ.M; |
| MatCopy(gJ.M, J, SAME_NONZERO_PATTERN); | | MatCopy(gJ.M, J, SAME_NONZERO_PATTERN); |
|
| //Message::Barrier(); | | // Message::Barrier(); |
| return 0; | | return 0; |
| } | | } |
| #endif | | #endif |
| | |
|
| static PetscErrorCode _mySnesMonitor(SNES snes, PetscInt it, PetscReal rnorm, vo | | static PetscErrorCode _mySnesMonitor(SNES snes, PetscInt it, PetscReal rnorm, |
| id *mctx) | | void *mctx) |
| { | | { |
| Message::Info("%3ld SNES Residual norm %14.12e", (long)it, rnorm); | | Message::Info("%3ld SNES Residual norm %14.12e", (long)it, rnorm); |
| return 0; | | return 0; |
| } | | } |
| | |
|
| static void _solveNL(gMatrix *A, gVector *B, gMatrix *J, gVector *R, gSolver *So | | static void _solveNL(gMatrix *A, gVector *B, gMatrix *J, gVector *R, |
| lver, | | gSolver *Solver, gVector *X, int precond, int solverIndex) |
| gVector *X, int precond, int solverIndex) | | |
| { | | { |
|
| if(solverIndex < 0 || solverIndex > 9){ | | if(solverIndex < 0 || solverIndex > 9) { |
| Message::Error("NonLinear Solver index out of range (%d)", solverIndex); | | Message::Error("NonLinear Solver index out of range (%d)", solverIndex); |
| return; | | return; |
| } | | } |
| | |
| PetscInt n, m; | | PetscInt n, m; |
| _try(MatGetSize(J->M, &n, &m)); | | _try(MatGetSize(J->M, &n, &m)); |
|
| if(!n){ | | if(!n) { |
| Message::Warning("Zero-size jacobian: skipping solve!"); | | Message::Warning("Zero-size jacobian: skipping solve!"); |
| return; | | return; |
| } | | } |
| | |
| bool view = !Solver->snes[solverIndex]; | | bool view = !Solver->snes[solverIndex]; |
| | |
| // either we are on sequential (!GetIsCommWorld) or in parallel with rank = 0 | | // either we are on sequential (!GetIsCommWorld) or in parallel with rank = 0 |
| // (GetIsCommWorld) | | // (GetIsCommWorld) |
| if(view && (!Message::GetCommRank() || !Message::GetIsCommWorld())) | | if(view && (!Message::GetCommRank() || !Message::GetIsCommWorld())) |
| Message::Info("N: %ld", (long)n); | | Message::Info("N: %ld", (long)n); |
| | |
| if(solverIndex != 0) | | if(solverIndex != 0) |
| Message::Info("Using nonlinear solver index %d", solverIndex); | | Message::Info("Using nonlinear solver index %d", solverIndex); |
| | |
| // Setting nonlinear solver defaults | | // Setting nonlinear solver defaults |
| if(!Solver->snes[solverIndex]) { | | if(!Solver->snes[solverIndex]) { |
| _try(SNESCreate(MyComm, &Solver->snes[solverIndex])); | | _try(SNESCreate(MyComm, &Solver->snes[solverIndex])); |
|
| _try(SNESMonitorSet(Solver->snes[solverIndex], _mySnesMonitor, | | _try(SNESMonitorSet(Solver->snes[solverIndex], _mySnesMonitor, PETSC_NULL, |
| PETSC_NULL, PETSC_NULL)); | | PETSC_NULL)); |
| /* //kj+++ | | |
| PetscReal abstol = 1.e-4; //(PETSC_DEFAULT=1.e-15) | | |
| PetscReal rtol = 1.e-2; //(PETSC_DEFAULT=1.e-8) | | |
| PetscReal stol = 1.e-2; //(PETSC_DEFAULT=1.e-8) | | |
| PetscInt maxit = 50; //(PETSC_DEFAULT=50) | | |
| PetscInt maxf = 10000; //(PETSC_DEFAULT=10000) | | |
| _try(SNESSetTolerances(Solver->snes[solverIndex], abstol, rtol, | | |
| stol, maxit, maxf)); | | |
| */ | | |
| _try(SNESSetTolerances(Solver->snes[solverIndex], 1.e-12, PETSC_DEFAULT, | | _try(SNESSetTolerances(Solver->snes[solverIndex], 1.e-12, PETSC_DEFAULT, |
| PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); | | PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT)); |
| | |
| // override default options with those from database (if any) | | // override default options with those from database (if any) |
|
| //_try(SNESSetType(Solver->snes[solverIndex],SNESNEWTONLS)); // kj+++ | | |
| //_try(SNESQNSetType(Solver->snes[solverIndex], SNES_QN_LBFGS)); // kj+++ | | |
| _try(SNESSetFromOptions(Solver->snes[solverIndex])); | | _try(SNESSetFromOptions(Solver->snes[solverIndex])); |
| | |
|
| /* //kj+++ | | PetscTruth fd_jacobian = PETSC_FALSE, snes_fd = PETSC_FALSE; |
| SNESType mySNESType; | | |
| _try(SNESGetType(Solver->snes[solverIndex],&mySNESType)); | | |
| //Message::Info("Try to show Type"); | | |
| Message::Info("SNESType: %s", mySNESType); | | |
| */ | | |
| | | |
| PetscTruth fd_jacobian = PETSC_FALSE, snes_fd = PETSC_FALSE ; | | |
| PetscOptionsGetTruth(PETSC_NULL, "-fd_jacobian", &fd_jacobian, 0); | | PetscOptionsGetTruth(PETSC_NULL, "-fd_jacobian", &fd_jacobian, 0); |
| PetscOptionsGetTruth(PETSC_NULL, "-snes_fd", &snes_fd, 0); | | PetscOptionsGetTruth(PETSC_NULL, "-snes_fd", &snes_fd, 0); |
|
| if (fd_jacobian || snes_fd) { | | if(fd_jacobian || snes_fd) { |
| // Message::Error("Finite Difference Jacobian not yet implemented"); | | #if(PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 4) |
| | |
| #if (PETSC_VERSION_MAJOR == 3) && (PETSC_VERSION_MINOR >= 4) | | |
| _try(SNESSetJacobian(Solver->snes[solverIndex], J->M, J->M, | | _try(SNESSetJacobian(Solver->snes[solverIndex], J->M, J->M, |
| SNESComputeJacobianDefault, PETSC_NULL)); | | SNESComputeJacobianDefault, PETSC_NULL)); |
| #else | | #else |
| _try(SNESSetJacobian(Solver->snes[solverIndex], J->M, J->M, | | _try(SNESSetJacobian(Solver->snes[solverIndex], J->M, J->M, |
| SNESDefaultComputeJacobian, PETSC_NULL)); | | SNESDefaultComputeJacobian, PETSC_NULL)); |
| #endif | | #endif |
| } | | } |
| else { | | else { |
| Message::Info("Jacobian computed by GetDP"); | | Message::Info("Jacobian computed by GetDP"); |
|
| _try(SNESSetJacobian(Solver->snes[solverIndex], J->M, J->M, _NLFormJacobia | | _try(SNESSetJacobian(Solver->snes[solverIndex], J->M, J->M, |
| n, | | _NLFormJacobian, PETSC_NULL)); |
| PETSC_NULL)); | | |
| } | | } |
| _try(SNESSetFunction(Solver->snes[solverIndex], R->V, _NLFormFunction, | | _try(SNESSetFunction(Solver->snes[solverIndex], R->V, _NLFormFunction, |
| PETSC_NULL)); // R(x) = A(x)*x-b | | PETSC_NULL)); // R(x) = A(x)*x-b |
| } | | } |
| | |
| KSP ksp; | | KSP ksp; |
| SNESGetKSP(Solver->snes[solverIndex], &ksp); | | SNESGetKSP(Solver->snes[solverIndex], &ksp); |
| PC pc; | | PC pc; |
| _try(KSPGetPC(ksp, &pc)); | | _try(KSPGetPC(ksp, &pc)); |
| _try(KSPSetType(ksp, "preonly")); | | _try(KSPSetType(ksp, "preonly")); |
| _try(PCSetType(pc, PCLU)); | | _try(PCSetType(pc, PCLU)); |
|
| #if (PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) | | #if(PETSC_VERSION_MAJOR > 2) && defined(PETSC_HAVE_MUMPS) |
| _try(PCFactorSetMatSolverPackage(pc, "mumps")); | | _try(PCFactorSetMatSolverPackage(pc, "mumps")); |
| #endif | | #endif |
| _try(SNESSolve(Solver->snes[solverIndex], PETSC_NULL, X->V)); | | _try(SNESSolve(Solver->snes[solverIndex], PETSC_NULL, X->V)); |
| | |
| // copy result on all procs | | // copy result on all procs |
| _fillseq(X); | | _fillseq(X); |
| | |
| if(view && Message::GetVerbosity() > 5) | | if(view && Message::GetVerbosity() > 5) |
| _try(SNESView(Solver->snes[solverIndex], MyPetscViewer)); | | _try(SNESView(Solver->snes[solverIndex], MyPetscViewer)); |
| | |
|
| if(!Message::GetCommRank() || !Message::GetIsCommWorld()){ | | if(!Message::GetCommRank() || !Message::GetIsCommWorld()) { |
| PetscInt its; | | PetscInt its; |
| _try(SNESGetIterationNumber(Solver->snes[solverIndex], &its)); | | _try(SNESGetIterationNumber(Solver->snes[solverIndex], &its)); |
| Message::Info("Number of Newton iterations %d", its); | | Message::Info("Number of Newton iterations %d", its); |
| } | | } |
| } | | } |
| | |
| void LinAlg_SolveNL(gMatrix *A, gVector *B, gMatrix *J, gVector *R, | | void LinAlg_SolveNL(gMatrix *A, gVector *B, gMatrix *J, gVector *R, |
| gSolver *Solver, gVector *X, int solverIndex) | | gSolver *Solver, gVector *X, int solverIndex) |
| { | | { |
| _solveNL(A, B, J, R, Solver, X, 1, solverIndex); | | _solveNL(A, B, J, R, Solver, X, 1, solverIndex); |
| | | |
| End of changes. 158 change blocks. |
|---|
| 267 lines changed or deleted | | 236 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 