Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Comparing changes

Choose two branches to see what's changed or to start a new pull request. If you need to, you can also compare across forks.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also compare across forks.
  • 6 commits
  • 37 files changed
  • 0 commit comments
  • 2 contributors
Commits on Jan 30, 2014
@BarrySmith BarrySmith merge petscsnesfas.h into petscsnes.h a055c8c
@BarrySmith BarrySmith Merge branch 'barry/rm-petscsnesfas_h' into next bcd0eba
@BarrySmith BarrySmith removed PCASA and supporting DMADDA
they don't work and if needed need to be rewritten
649e3ec
@BarrySmith BarrySmith Merge branch 'barry/rm-pcasa' into next
Conflicts:
	src/dm/impls/makefile
	src/ksp/ksp/examples/tutorials/ex38.c
	src/ksp/ksp/examples/tutorials/ex39.c
	src/ksp/pc/impls/makefile
78855ac
@jedbrown jedbrown PetscMalloc: allow ptr=malloc(0) and free(ptr)
  ISO C90 §7.10.3: "If the size of the space requested is zero, the
  behavior is implementation-defined; the value returned shall be either
  a null pointer or a unique pointer."

We believe that ptr=malloc(0) and free(ptr) is now handled properly, so
it's no longer necessary to short-circuit our PetscMalloc.  This allows
memory checking tools to match size-0 allocations and frees if they so
choose.
49d7da5
@jedbrown jedbrown Merge branch 'jed/malloc-zero' into next
* jed/malloc-zero:
  PetscMalloc: allow ptr=malloc(0) and free(ptr)
30fbb49
Showing with 82 additions and 3,957 deletions.
  1. +1 −1  include/finclude/makefile
  2. +0 −1  include/finclude/petsc.h
  3. +0 −1  include/finclude/petsc.h90
  4. +0 −11 include/finclude/petscdmadda.h90
  5. +0 −1  include/finclude/petscpcdef.h
  6. +0 −2  include/petsc.h
  7. +0 −24 include/petscdmadda.h
  8. +0 −1  include/petscpc.h
  9. +0 −13 include/petscpcasa.h
  10. +63 −0 include/petscsnes.h
  11. +0 −70 include/petscsnesfas.h
  12. +8 −6 include/petscsys.h
  13. +0 −4 src/dm/f90-mod/petscdmmod.F
  14. +0 −823 src/dm/impls/adda/adda.c
  15. +0 −30 src/dm/impls/adda/addaimpl.h
  16. +0 −11 src/dm/impls/adda/examples/makefile
  17. +0 −24 src/dm/impls/adda/examples/tests/ex1.c
  18. +0 −19 src/dm/impls/adda/examples/tests/makefile
  19. +0 −16 src/dm/impls/adda/makefile
  20. +1 −1  src/dm/impls/makefile
  21. +0 −2  src/dm/interface/dmregall.c
  22. +0 −1  src/docs/website/miscellaneous/acknwldg.html
  23. +0 −109 src/ksp/ksp/examples/tutorials/ex38.c
  24. +0 −256 src/ksp/ksp/examples/tutorials/ex39.c
  25. +0 −302 src/ksp/ksp/examples/tutorials/ex40.c
  26. +3 −21 src/ksp/ksp/examples/tutorials/makefile
  27. +0 −2,047 src/ksp/pc/impls/asa/asa.c
  28. +0 −130 src/ksp/pc/impls/asa/asa.h
  29. +0 −16 src/ksp/pc/impls/asa/makefile
  30. +1 −1  src/ksp/pc/impls/makefile
  31. +0 −2  src/ksp/pc/interface/pcregis.c
  32. +1 −1  src/snes/impls/fas/fas.c
  33. +1 −1  src/snes/impls/fas/fasfunc.c
  34. +1 −1  src/snes/impls/fas/fasgalerkin.c
  35. +0 −1  src/snes/impls/fas/fasimpls.h
  36. +1 −6 src/sys/memory/mtr.c
  37. +1 −1  src/ts/impls/implicit/theta/theta.c
View
2  include/finclude/makefile
@@ -8,7 +8,7 @@ SOURCEC =
SOURCEF =
SOURCEH = petsc.h petscsys.h petsclog.h petscvec.h petscsnes.h petscdm.h petscdmda.h petscdraw.h petscmat.h \
petscksp.h petscpc.h petscviewer.h petscis.h petscao.h \
- petscts.h petscis.h90 petscvec.h90 petscmat.h90 petscdm.h90 petscdmda.h90 petscdmadda.h90 petscdmcomposite.h90 \
+ petscts.h petscis.h90 petscvec.h90 petscmat.h90 petscdm.h90 petscdmda.h90 petscdmcomposite.h90 \
petscdmredundant.h90 \
petscdef.h petscsysdef.h petsclogdef.h petscvecdef.h petscsnesdef.h petscdmdef.h petscdmdadef.h \
petscdrawdef.h petscmatdef.h petsckspdef.h petscpcdef.h petscviewerdef.h petscisdef.h petscaodef.h \
View
1  include/finclude/petsc.h
@@ -15,7 +15,6 @@
#include "finclude/petscdm.h"
#include "finclude/petscdmda.h"
#include "finclude/petscdmplex.h"
-!#include "finclude/petscdmadda.h"
!#include "finclude/petscdmcomposite.h"
!#include "finclude/petscdmsliced.h"
#include "finclude/petscts.h"
View
1  include/finclude/petsc.h90
@@ -9,7 +9,6 @@
#include "finclude/petscdm.h90"
#include "finclude/petscdmda.h90"
#include "finclude/petscdmplex.h90"
-#include "finclude/petscdmadda.h90"
#include "finclude/petscdmcomposite.h90"
#include "finclude/petscdmredundant.h90"
#include "finclude/petscpc.h90"
View
11 include/finclude/petscdmadda.h90
@@ -1,11 +0,0 @@
-!
-!
-! Additional DM include file for use of PETSc with Fortran 90/HPF
-!
-#include "finclude/ftn-custom/petscdm.h90"
-#if defined(PETSC_USE_FORTRAN_INTERFACES)
- interface
-#include "finclude/ftn-auto/petscdmadda.h90"
- end interface
-#endif
-
View
1  include/finclude/petscpcdef.h
@@ -58,7 +58,6 @@
#define PCGALERKIN 'galerkin'
#define PCEXOTIC 'exotic'
#define PCSUPPORTGRAPH 'supportgraph'
-#define PCASA 'asa'
#define PCCP 'cp'
#define PCBFBT 'bfbt'
#define PCLSC 'lsc'
View
2  include/petsc.h
@@ -11,7 +11,6 @@
#include <petscdraw.h>
#include <petscdmda.h>
-#include <petscdmadda.h>
#include <petscdmcomposite.h>
#include <petscdmpatch.h>
#include <petscdmplex.h>
@@ -23,7 +22,6 @@
#include <petsccharacteristic.h>
-#include <petscsnesfas.h>
#include <petscts.h>
#include <taosolver.h>
View
24 include/petscdmadda.h
@@ -1,24 +0,0 @@
-#if !defined(__PETSCDMDA_H)
-#define __PETSCDMDA_H
-
-#include <petscdm.h>
-
-PetscErrorCode DMADDACreate(MPI_Comm,PetscInt,PetscInt*,PetscInt*,PetscInt,PetscBool *,DM*);
-PetscErrorCode DMADDASetParameters(DM,PetscInt,PetscInt*,PetscInt*,PetscInt,PetscBool*);
-PetscErrorCode DMADDASetRefinement(DM, PetscInt *,PetscInt);
-PetscErrorCode DMADDAGetCorners(DM, PetscInt **, PetscInt **);
-PetscErrorCode DMADDAGetGhostCorners(DM, PetscInt **, PetscInt **);
-PetscErrorCode DMADDAGetMatrixNS(DM, DM, MatType , Mat *);
-
-/* functions to set values in vectors and matrices */
-struct _ADDAIdx_s {
- PetscInt *x; /* the coordinates, user has to make sure it is the correct size! */
- PetscInt d; /* indexes the dof */
-};
-typedef struct _ADDAIdx_s ADDAIdx;
-
-PetscErrorCode DMADDAMatSetValues(Mat, DM, PetscInt, const ADDAIdx[], DM, PetscInt, const ADDAIdx[], const PetscScalar[], InsertMode);
-PetscBool ADDAHCiterStartup(const PetscInt, const PetscInt *const, const PetscInt *const, PetscInt *const);
-PetscBool ADDAHCiter(const PetscInt, const PetscInt *const, const PetscInt *const, PetscInt *const);
-
-#endif
View
1  include/petscpc.h
@@ -64,7 +64,6 @@ typedef const char* PCType;
#define PCML "ml"
#define PCGALERKIN "galerkin"
#define PCEXOTIC "exotic"
-#define PCASA "asa"
#define PCCP "cp"
#define PCBFBT "bfbt"
#define PCLSC "lsc"
View
13 include/petscpcasa.h
@@ -1,13 +0,0 @@
-/*
- Defines the interface functions for the ASA geometric preconditioner
-*/
-#ifndef __PETSCPCASA_H
-#define __PETSCPCASA_H
-#include <petscksp.h>
-#include <petscdm.h>
-
-
-PETSC_EXTERN PetscErrorCode PCASASetTolerances(PC, PetscReal,PetscReal, PetscReal, PetscInt);
-PETSC_EXTERN PetscErrorCode PCSolve_ASA(PC, Vec, Vec);
-
-#endif
View
63 include/petscsnes.h
@@ -733,4 +733,67 @@ PETSC_EXTERN PetscErrorCode SNESCompositeAddSNES(SNES,SNESType);
PETSC_EXTERN PetscErrorCode SNESCompositeGetSNES(SNES,PetscInt,SNES *);
PETSC_EXTERN PetscErrorCode SNESCompositeSetDamping(SNES,PetscInt,PetscReal);
+/*E
+ SNESFASType - Determines the type of nonlinear multigrid method that is run.
+
+ Level: beginner
+
+ Values:
++ SNES_FAS_MULTIPLICATIVE (default) - traditional V or W cycle as determined by SNESFASSetCycles()
+. SNES_FAS_ADDITIVE - additive FAS cycle
+. SNES_FAS_FULL - full FAS cycle
+- SNES_FAS_KASKADE - Kaskade FAS cycle
+.seealso: PCMGSetType(), PCMGType
+
+E*/
+typedef enum { SNES_FAS_MULTIPLICATIVE, SNES_FAS_ADDITIVE, SNES_FAS_FULL, SNES_FAS_KASKADE } SNESFASType;
+PETSC_EXTERN const char *const SNESFASTypes[];
+
+/* called on the finest level FAS instance*/
+PETSC_EXTERN PetscErrorCode SNESFASSetType(SNES, SNESFASType);
+PETSC_EXTERN PetscErrorCode SNESFASGetType(SNES, SNESFASType*);
+PETSC_EXTERN PetscErrorCode SNESFASSetLevels(SNES, PetscInt, MPI_Comm *);
+PETSC_EXTERN PetscErrorCode SNESFASGetLevels(SNES, PetscInt *);
+PETSC_EXTERN PetscErrorCode SNESFASGetCycleSNES(SNES, PetscInt, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASSetNumberSmoothUp(SNES, PetscInt);
+PETSC_EXTERN PetscErrorCode SNESFASSetNumberSmoothDown(SNES, PetscInt);
+PETSC_EXTERN PetscErrorCode SNESFASSetCycles(SNES, PetscInt);
+PETSC_EXTERN PetscErrorCode SNESFASSetMonitor(SNES, PetscBool);
+PETSC_EXTERN PetscErrorCode SNESFASSetLog(SNES, PetscBool);
+
+PETSC_EXTERN PetscErrorCode SNESFASSetGalerkin(SNES, PetscBool);
+PETSC_EXTERN PetscErrorCode SNESFASGetGalerkin(SNES, PetscBool*);
+
+/* called on any level -- "Cycle" FAS instance */
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetSmoother(SNES, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetSmootherUp(SNES, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetSmootherDown(SNES, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetCorrection(SNES, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetInterpolation(SNES, Mat*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetRestriction(SNES, Mat*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetInjection(SNES, Mat*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleGetRScale(SNES, Vec*);
+PETSC_EXTERN PetscErrorCode SNESFASCycleSetCycles(SNES, PetscInt);
+PETSC_EXTERN PetscErrorCode SNESFASCycleIsFine(SNES, PetscBool*);
+
+/* called on the (outer) finest level FAS to set/get parameters on any level instance */
+PETSC_EXTERN PetscErrorCode SNESFASSetInterpolation(SNES, PetscInt, Mat);
+PETSC_EXTERN PetscErrorCode SNESFASGetInterpolation(SNES, PetscInt, Mat*);
+PETSC_EXTERN PetscErrorCode SNESFASSetRestriction(SNES, PetscInt, Mat);
+PETSC_EXTERN PetscErrorCode SNESFASGetRestriction(SNES, PetscInt, Mat*);
+PETSC_EXTERN PetscErrorCode SNESFASSetInjection(SNES, PetscInt, Mat);
+PETSC_EXTERN PetscErrorCode SNESFASGetInjection(SNES, PetscInt, Mat*);
+PETSC_EXTERN PetscErrorCode SNESFASSetRScale(SNES, PetscInt, Vec);
+PETSC_EXTERN PetscErrorCode SNESFASGetRScale(SNES, PetscInt, Vec*);
+
+PETSC_EXTERN PetscErrorCode SNESFASGetSmoother(SNES, PetscInt, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASGetSmootherUp(SNES, PetscInt, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASGetSmootherDown(SNES, PetscInt, SNES*);
+PETSC_EXTERN PetscErrorCode SNESFASGetCoarseSolve(SNES, SNES*);
+
+/* parameters for full FAS */
+PETSC_EXTERN PetscErrorCode SNESFASFullSetDownSweep(SNES,PetscBool);
+PETSC_EXTERN PetscErrorCode SNESFASCreateCoarseVec(SNES,Vec*);
+PETSC_EXTERN PetscErrorCode SNESFASRestrict(SNES,Vec,Vec);
+
#endif
View
70 include/petscsnesfas.h
@@ -1,70 +0,0 @@
-#if !(defined __PETSCSNESFAS_H)
-#define __PETSCSNESFAS_H
-#include <petscsnes.h>
-
-
-/*E
- SNESFASType - Determines the type of nonlinear multigrid method that is run.
-
- Level: beginner
-
- Values:
-+ SNES_FAS_MULTIPLICATIVE (default) - traditional V or W cycle as determined by SNESFASSetCycles()
-. SNES_FAS_ADDITIVE - additive FAS cycle
-. SNES_FAS_FULL - full FAS cycle
-- SNES_FAS_KASKADE - Kaskade FAS cycle
-.seealso: PCMGSetType(), PCMGType
-
-E*/
-typedef enum { SNES_FAS_MULTIPLICATIVE, SNES_FAS_ADDITIVE, SNES_FAS_FULL, SNES_FAS_KASKADE } SNESFASType;
-PETSC_EXTERN const char *const SNESFASTypes[];
-
-/* called on the finest level FAS instance*/
-PETSC_EXTERN PetscErrorCode SNESFASSetType(SNES, SNESFASType);
-PETSC_EXTERN PetscErrorCode SNESFASGetType(SNES, SNESFASType*);
-PETSC_EXTERN PetscErrorCode SNESFASSetLevels(SNES, PetscInt, MPI_Comm *);
-PETSC_EXTERN PetscErrorCode SNESFASGetLevels(SNES, PetscInt *);
-PETSC_EXTERN PetscErrorCode SNESFASGetCycleSNES(SNES, PetscInt, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASSetNumberSmoothUp(SNES, PetscInt);
-PETSC_EXTERN PetscErrorCode SNESFASSetNumberSmoothDown(SNES, PetscInt);
-PETSC_EXTERN PetscErrorCode SNESFASSetCycles(SNES, PetscInt);
-PETSC_EXTERN PetscErrorCode SNESFASSetMonitor(SNES, PetscBool);
-PETSC_EXTERN PetscErrorCode SNESFASSetLog(SNES, PetscBool);
-
-
-PETSC_EXTERN PetscErrorCode SNESFASSetGalerkin(SNES, PetscBool);
-PETSC_EXTERN PetscErrorCode SNESFASGetGalerkin(SNES, PetscBool*);
-
-/* called on any level -- "Cycle" FAS instance */
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetSmoother(SNES, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetSmootherUp(SNES, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetSmootherDown(SNES, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetCorrection(SNES, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetInterpolation(SNES, Mat*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetRestriction(SNES, Mat*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetInjection(SNES, Mat*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleGetRScale(SNES, Vec*);
-PETSC_EXTERN PetscErrorCode SNESFASCycleSetCycles(SNES, PetscInt);
-PETSC_EXTERN PetscErrorCode SNESFASCycleIsFine(SNES, PetscBool*);
-
-/* called on the (outer) finest level FAS to set/get parameters on any level instance */
-PETSC_EXTERN PetscErrorCode SNESFASSetInterpolation(SNES, PetscInt, Mat);
-PETSC_EXTERN PetscErrorCode SNESFASGetInterpolation(SNES, PetscInt, Mat*);
-PETSC_EXTERN PetscErrorCode SNESFASSetRestriction(SNES, PetscInt, Mat);
-PETSC_EXTERN PetscErrorCode SNESFASGetRestriction(SNES, PetscInt, Mat*);
-PETSC_EXTERN PetscErrorCode SNESFASSetInjection(SNES, PetscInt, Mat);
-PETSC_EXTERN PetscErrorCode SNESFASGetInjection(SNES, PetscInt, Mat*);
-PETSC_EXTERN PetscErrorCode SNESFASSetRScale(SNES, PetscInt, Vec);
-PETSC_EXTERN PetscErrorCode SNESFASGetRScale(SNES, PetscInt, Vec*);
-
-PETSC_EXTERN PetscErrorCode SNESFASGetSmoother(SNES, PetscInt, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASGetSmootherUp(SNES, PetscInt, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASGetSmootherDown(SNES, PetscInt, SNES*);
-PETSC_EXTERN PetscErrorCode SNESFASGetCoarseSolve(SNES, SNES*);
-
-/* parameters for full FAS */
-PETSC_EXTERN PetscErrorCode SNESFASFullSetDownSweep(SNES,PetscBool);
-PETSC_EXTERN PetscErrorCode SNESFASCreateCoarseVec(SNES,Vec*);
-PETSC_EXTERN PetscErrorCode SNESFASRestrict(SNES,Vec,Vec);
-
-#endif
View
14 include/petscsys.h
@@ -516,17 +516,17 @@ PETSC_EXTERN PetscErrorCode PetscCommDestroy(MPI_Comm*);
Level: beginner
- Notes: Memory is always allocated at least double aligned
+ Notes:
+ Memory is always allocated at least double aligned
- If you request memory of zero size it will allocate no space and assign the pointer to 0; PetscFree() will
- properly handle not freeing the null pointer.
+ It is safe to allocate size 0 and pass the resulting pointer (which may or may not be NULL) to PetscFree().
.seealso: PetscFree(), PetscNew()
Concepts: memory allocation
M*/
-#define PetscMalloc(a,b) ((a != 0) ? (*PetscTrMalloc)((a),__LINE__,PETSC_FUNCTION_NAME,__FILE__,(void**)(b)) : (*(b) = 0,0) )
+#define PetscMalloc(a,b) ((*PetscTrMalloc)((a),__LINE__,PETSC_FUNCTION_NAME,__FILE__,(void**)(b)))
/*MC
PetscAddrAlign - Rounds up an address to PETSC_MEMALIGN alignment
@@ -1072,14 +1072,16 @@ M*/
Level: beginner
- Notes: Memory must have been obtained with PetscNew() or PetscMalloc()
+ Notes:
+ Memory must have been obtained with PetscNew() or PetscMalloc().
+ It is safe to call PetscFree() on a NULL pointer.
.seealso: PetscNew(), PetscMalloc(), PetscFreeVoid()
Concepts: memory allocation
M*/
-#define PetscFree(a) ((a) && ((*PetscTrFree)((void*)(a),__LINE__,PETSC_FUNCTION_NAME,__FILE__) || ((a) = 0,0)))
+#define PetscFree(a) ((*PetscTrFree)((void*)(a),__LINE__,PETSC_FUNCTION_NAME,__FILE__) || ((a) = 0,0))
/*MC
PetscFreeVoid - Frees memory
View
4 src/dm/f90-mod/petscdmmod.F
@@ -41,7 +41,3 @@ module petscdmcomposite
#include <finclude/petscdmcomposite.h90>
end module
- module petscdmadda
- use petscdm
-#include <finclude/petscdmadda.h90>
- end module
View
823 src/dm/impls/adda/adda.c
@@ -1,823 +0,0 @@
-/*
-
- Contributed by Arvid Bessen, Columbia University, June 2007
-
- Extension of DMDA object to any number of dimensions.
-
-*/
-#include <../src/dm/impls/adda/addaimpl.h> /*I "petscdmadda.h" I*/
-
-
-#undef __FUNCT__
-#define __FUNCT__ "DMDestroy_ADDA"
-PetscErrorCode DMDestroy_ADDA(DM dm)
-{
- PetscErrorCode ierr;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
-
- PetscFunctionBegin;
- /* destroy the allocated data */
- ierr = PetscFree(dd->nodes);CHKERRQ(ierr);
- ierr = PetscFree(dd->procs);CHKERRQ(ierr);
- ierr = PetscFree(dd->lcs);CHKERRQ(ierr);
- ierr = PetscFree(dd->lce);CHKERRQ(ierr);
- ierr = PetscFree(dd->lgs);CHKERRQ(ierr);
- ierr = PetscFree(dd->lge);CHKERRQ(ierr);
- ierr = PetscFree(dd->refine);CHKERRQ(ierr);
-
- ierr = VecDestroy(&dd->global);CHKERRQ(ierr);
- /* This was originally freed in DMDestroy(), but that prevents reference counting of backend objects */
- ierr = PetscFree(dd);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMView_ADDA"
-PetscErrorCode DMView_ADDA(DM dm, PetscViewer v)
-{
- PetscFunctionBegin;
- SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP, "Not implemented yet");
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreateGlobalVector_ADDA"
-PetscErrorCode DMCreateGlobalVector_ADDA(DM dm, Vec *vec)
-{
- PetscErrorCode ierr;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm,DM_CLASSID,1);
- PetscValidPointer(vec,2);
- ierr = VecDuplicate(dd->global, vec);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreateColoring_ADDA"
-PetscErrorCode DMCreateColoring_ADDA(DM dm, ISColoringType ctype,ISColoring *coloring)
-{
- PetscFunctionBegin;
- SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP, "Not implemented yet");
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreateMatrix_ADDA"
-PetscErrorCode DMCreateMatrix_ADDA(DM dm, Mat *mat)
-{
- PetscErrorCode ierr;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
- ierr = MatCreate(PetscObjectComm((PetscObject)dm), mat);CHKERRQ(ierr);
- ierr = MatSetSizes(*mat, dd->lsize, dd->lsize, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr);
- ierr = MatSetType(*mat, dm->mattype);CHKERRQ(ierr);
- ierr = MatSetUp(*mat);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDAGetMatrixNS"
-/*@
- DMADDAGetMatrixNS - Creates matrix compatiable with two distributed arrays
-
- Collective on ADDA
-
- Input Parameter:
-. addar - the distributed array for which we create the matrix, which indexes the rows
-. addac - the distributed array for which we create the matrix, which indexes the columns
-- mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, or
- any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).
-
- Output Parameter:
-. mat - the empty Jacobian
-
- Level: beginner
-
-.keywords: distributed array, matrix
-
-.seealso: DMCreateMatrix()
-@*/
-PetscErrorCode DMADDAGetMatrixNS(DM dm, DM dmc, MatType mtype, Mat *mat)
-{
- PetscErrorCode ierr;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
- DM_ADDA *ddc = (DM_ADDA*)dmc->data;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
- PetscValidHeaderSpecific(dmc, DM_CLASSID, 2);
- PetscCheckSameComm(dm, 1, dmc, 2);
- ierr = MatCreate(PetscObjectComm((PetscObject)dm), mat);CHKERRQ(ierr);
- ierr = MatSetSizes(*mat, dd->lsize, ddc->lsize, PETSC_DECIDE, PETSC_DECIDE);CHKERRQ(ierr);
- ierr = MatSetType(*mat, mtype);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreateInterpolation_ADDA"
-PetscErrorCode DMCreateInterpolation_ADDA(DM dm1,DM dm2,Mat *mat,Vec *vec)
-{
- PetscFunctionBegin;
- SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP, "Not implemented yet");
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMRefine_ADDA"
-PetscErrorCode DMRefine_ADDA(DM dm, MPI_Comm comm, DM *dmf)
-{
- PetscFunctionBegin;
- SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP, "Not implemented yet");
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCoarsen_ADDA"
-PetscErrorCode DMCoarsen_ADDA(DM dm, MPI_Comm comm,DM *dmc)
-{
- PetscErrorCode ierr;
- PetscInt *nodesc;
- PetscInt dofc;
- PetscInt i;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
- PetscValidPointer(dmc, 3);
- ierr = PetscMalloc1(dd->dim, &nodesc);CHKERRQ(ierr);
- for (i=0; i<dd->dim; i++) {
- nodesc[i] = (dd->nodes[i] % dd->refine[i]) ? dd->nodes[i] / dd->refine[i] + 1 : dd->nodes[i] / dd->refine[i];
- }
- dofc = (dd->dof % dd->dofrefine) ? dd->dof / dd->dofrefine + 1 : dd->dof / dd->dofrefine;
- ierr = DMADDACreate(PetscObjectComm((PetscObject)dm), dd->dim, nodesc, dd->procs, dofc, dd->periodic, dmc);CHKERRQ(ierr);
- ierr = PetscFree(nodesc);CHKERRQ(ierr);
- /* copy refinement factors */
- ierr = DMADDASetRefinement(*dmc, dd->refine, dd->dofrefine);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreateInjection_ADDA"
-PetscErrorCode DMCreateInjection_ADDA(DM dm1,DM dm2, VecScatter *ctx)
-{
- PetscFunctionBegin;
- SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP, "Not implemented yet");
- PetscFunctionReturn(0);
-}
-
-/*@C
- ADDAHCiterStartup - performs the first check for an iteration through a hypercube
- lc, uc, idx all have to be valid arrays of size dim
- This function sets idx to lc and then checks, whether the lower corner (lc) is less
- than thre upper corner (uc). If lc "<=" uc in all coordinates, it returns PETSC_TRUE,
- and PETSC_FALSE otherwise.
-
- Input Parameters:
-+ dim - the number of dimension
-. lc - the "lower" corner
-- uc - the "upper" corner
-
- Output Parameters:
-. idx - the index that this function increases
-
- Developer Notes: This code is crap! You cannot return a value and NO ERROR code in PETSc!
-
- Level: developer
-@*/
-PetscBool ADDAHCiterStartup(const PetscInt dim, const PetscInt *const lc, const PetscInt *const uc, PetscInt *const idx)
-{
- PetscErrorCode ierr;
- PetscInt i;
-
- ierr = PetscMemcpy(idx, lc, sizeof(PetscInt)*dim);
- if (ierr) {
- PetscError(PETSC_COMM_SELF,__LINE__,__FUNCT__,__FILE__,ierr,PETSC_ERROR_REPEAT," ");
- return PETSC_FALSE;
- }
- for (i=0; i<dim; i++) {
- if (lc[i] > uc[i]) return PETSC_FALSE;
- }
- return PETSC_TRUE;
-}
-
-/*@C
- ADDAHCiter - iterates through a hypercube
- lc, uc, idx all have to be valid arrays of size dim
- This function return PETSC_FALSE, if idx exceeds uc, PETSC_TRUE otherwise.
- There are no guarantees on what happens if idx is not in the hypercube
- spanned by lc, uc, this should be checked with ADDAHCiterStartup.
-
- Use this code as follows:
- if (ADDAHCiterStartup(dim, lc, uc, idx)) {
- do {
- ...
- } while (ADDAHCiter(dim, lc, uc, idx));
- }
-
- Input Parameters:
-+ dim - the number of dimension
-. lc - the "lower" corner
-- uc - the "upper" corner
-
- Output Parameters:
-. idx - the index that this function increases
-
- Level: developer
-@*/
-PetscBool ADDAHCiter(const PetscInt dim, const PetscInt *const lc, const PetscInt *const uc, PetscInt *const idx)
-{
- PetscInt i;
- for (i=dim-1; i>=0; i--) {
- idx[i] += 1;
- if (uc[i] > idx[i]) {
- return PETSC_TRUE;
- } else {
- idx[i] -= uc[i] - lc[i];
- }
- }
- return PETSC_FALSE;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreateAggregates_ADDA"
-PetscErrorCode DMCreateAggregates_ADDA(DM dmc,DM dmf,Mat *rest)
-{
- PetscErrorCode ierr=0;
- PetscInt i;
- PetscInt dim;
- PetscInt dofc, doff;
- PetscInt *lcs_c, *lce_c;
- PetscInt *lcs_f, *lce_f;
- PetscInt *fgs, *fge;
- PetscInt fgdofs, fgdofe;
- ADDAIdx iter_c, iter_f;
- PetscInt max_agg_size;
- PetscMPIInt comm_size;
- ADDAIdx *fine_nodes;
- PetscInt fn_idx;
- PetscScalar *one_vec;
- DM_ADDA *ddc = (DM_ADDA*)dmc->data;
- DM_ADDA *ddf = (DM_ADDA*)dmf->data;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dmc, DM_CLASSID, 1);
- PetscValidHeaderSpecific(dmf, DM_CLASSID, 2);
- PetscValidPointer(rest,3);
- if (ddc->dim != ddf->dim) SETERRQ2(PetscObjectComm((PetscObject)dmf),PETSC_ERR_ARG_INCOMP,"Dimensions of ADDA do not match %D %D", ddc->dim, ddf->dim);CHKERRQ(ierr);
-/* if (dmc->dof != dmf->dof) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"DOF of ADDA do not match %D %D", dmc->dof, dmf->dof);CHKERRQ(ierr); */
- dim = ddc->dim;
- dofc = ddc->dof;
- doff = ddf->dof;
-
- ierr = DMADDAGetCorners(dmc, &lcs_c, &lce_c);CHKERRQ(ierr);
- ierr = DMADDAGetCorners(dmf, &lcs_f, &lce_f);CHKERRQ(ierr);
-
- /* compute maximum size of aggregate */
- max_agg_size = 1;
- for (i=0; i<dim; i++) {
- max_agg_size *= ddf->nodes[i] / ddc->nodes[i] + 1;
- }
- max_agg_size *= doff / dofc + 1;
-
- /* create the matrix that will contain the restriction operator */
- ierr = MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);CHKERRQ(ierr);
-
- /* construct matrix */
- if (comm_size == 1) {
- ierr = DMADDAGetMatrixNS(dmc, dmf, MATSEQAIJ, rest);CHKERRQ(ierr);
- ierr = MatSeqAIJSetPreallocation(*rest, max_agg_size, NULL);CHKERRQ(ierr);
- } else {
- ierr = DMADDAGetMatrixNS(dmc, dmf, MATMPIAIJ, rest);CHKERRQ(ierr);
- ierr = MatMPIAIJSetPreallocation(*rest, max_agg_size, NULL, max_agg_size, NULL);CHKERRQ(ierr);
- }
- /* store nodes in the fine grid here */
- ierr = PetscMalloc(sizeof(ADDAIdx)*max_agg_size, &fine_nodes);CHKERRQ(ierr);
- /* these are the values to set to, a collection of 1's */
- ierr = PetscMalloc(sizeof(PetscScalar)*max_agg_size, &one_vec);CHKERRQ(ierr);
- /* initialize */
- for (i=0; i<max_agg_size; i++) {
- ierr = PetscMalloc(sizeof(PetscInt)*dim, &(fine_nodes[i].x));CHKERRQ(ierr);
- one_vec[i] = 1.0;
- }
-
- /* get iterators */
- ierr = PetscMalloc(sizeof(PetscInt)*dim, &(iter_c.x));CHKERRQ(ierr);
- ierr = PetscMalloc(sizeof(PetscInt)*dim, &(iter_f.x));CHKERRQ(ierr);
-
- /* the fine grid node corner for each coarse grid node */
- ierr = PetscMalloc(sizeof(PetscInt)*dim, &fgs);CHKERRQ(ierr);
- ierr = PetscMalloc(sizeof(PetscInt)*dim, &fge);CHKERRQ(ierr);
-
- /* loop over all coarse nodes */
- ierr = PetscMemcpy(iter_c.x, lcs_c, sizeof(PetscInt)*dim);CHKERRQ(ierr);
- if (ADDAHCiterStartup(dim, lcs_c, lce_c, iter_c.x)) {
- do {
- /* find corresponding fine grid nodes */
- for (i=0; i<dim; i++) {
- fgs[i] = iter_c.x[i]*ddf->nodes[i]/ddc->nodes[i];
- fge[i] = PetscMin((iter_c.x[i]+1)*ddf->nodes[i]/ddc->nodes[i], ddf->nodes[i]);
- }
- /* treat all dof of the coarse grid */
- for (iter_c.d=0; iter_c.d<dofc; iter_c.d++) {
- /* find corresponding fine grid dof's */
- fgdofs = iter_c.d*doff/dofc;
- fgdofe = PetscMin((iter_c.d+1)*doff/dofc, doff);
- /* we now know the "box" of all the fine grid nodes that are mapped to one coarse grid node */
- fn_idx = 0;
- /* loop over those corresponding fine grid nodes */
- if (ADDAHCiterStartup(dim, fgs, fge, iter_f.x)) {
- do {
- /* loop over all corresponding fine grid dof */
- for (iter_f.d=fgdofs; iter_f.d<fgdofe; iter_f.d++) {
- ierr = PetscMemcpy(fine_nodes[fn_idx].x, iter_f.x, sizeof(PetscInt)*dim);CHKERRQ(ierr);
-
- fine_nodes[fn_idx].d = iter_f.d;
- fn_idx++;
- }
- } while (ADDAHCiter(dim, fgs, fge, iter_f.x));
- }
- /* add all these points to one aggregate */
- ierr = DMADDAMatSetValues(*rest, dmc, 1, &iter_c, dmf, fn_idx, fine_nodes, one_vec, INSERT_VALUES);CHKERRQ(ierr);
- }
- } while (ADDAHCiter(dim, lcs_c, lce_c, iter_c.x));
- }
-
- /* free memory */
- ierr = PetscFree(fgs);CHKERRQ(ierr);
- ierr = PetscFree(fge);CHKERRQ(ierr);
- ierr = PetscFree(iter_c.x);CHKERRQ(ierr);
- ierr = PetscFree(iter_f.x);CHKERRQ(ierr);
- ierr = PetscFree(lcs_c);CHKERRQ(ierr);
- ierr = PetscFree(lce_c);CHKERRQ(ierr);
- ierr = PetscFree(lcs_f);CHKERRQ(ierr);
- ierr = PetscFree(lce_f);CHKERRQ(ierr);
- ierr = PetscFree(one_vec);CHKERRQ(ierr);
- for (i=0; i<max_agg_size; i++) {
- ierr = PetscFree(fine_nodes[i].x);CHKERRQ(ierr);
- }
- ierr = PetscFree(fine_nodes);CHKERRQ(ierr);
-
- ierr = MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
- ierr = MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDASetRefinement"
-/*@
- DMADDASetRefinement - Sets the refinement factors of the distributed arrays.
-
- Collective on ADDA
-
- Input Parameter:
-+ adda - the ADDA object
-. refine - array of refinement factors
-- dofrefine - the refinement factor for the dof, usually just 1
-
- Level: developer
-
-.keywords: distributed array, refinement
-@*/
-PetscErrorCode DMADDASetRefinement(DM dm, PetscInt *refine, PetscInt dofrefine)
-{
- DM_ADDA *dd = (DM_ADDA*)dm->data;
- PetscErrorCode ierr;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
- PetscValidPointer(refine,3);
- ierr = PetscMemcpy(dd->refine, refine, dd->dim*sizeof(PetscInt));CHKERRQ(ierr);
- dd->dofrefine = dofrefine;
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDAGetCorners"
-/*@
- DMADDAGetCorners - Gets the corners of the local area
-
- Not Collective
-
- Input Parameter:
-. adda - the ADDA object
-
- Output Parameter:
-+ lcorner - the "lower" corner
-- ucorner - the "upper" corner
-
- Both lcorner and ucorner are allocated by this procedure and will point to an
- array of size dd->dim.
-
- Level: beginner
-
-.keywords: distributed array, refinement
-@*/
-PetscErrorCode DMADDAGetCorners(DM dm, PetscInt **lcorner, PetscInt **ucorner)
-{
- DM_ADDA *dd = (DM_ADDA*)dm->data;
- PetscErrorCode ierr;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
- PetscValidPointer(lcorner,2);
- PetscValidPointer(ucorner,3);
- ierr = PetscMalloc1(dd->dim, lcorner);CHKERRQ(ierr);
- ierr = PetscMalloc1(dd->dim, ucorner);CHKERRQ(ierr);
- ierr = PetscMemcpy(*lcorner, dd->lcs, dd->dim*sizeof(PetscInt));CHKERRQ(ierr);
- ierr = PetscMemcpy(*ucorner, dd->lce, dd->dim*sizeof(PetscInt));CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDAGetGhostCorners"
-/*@
- DMADDAGetGhostCorners - Gets the ghost corners of the local area
-
- Note Collective
-
- Input Parameter:
-. adda - the ADDA object
-
- Output Parameter:
-+ lcorner - the "lower" corner of the ghosted area
-- ucorner - the "upper" corner of the ghosted area
-
- Both lcorner and ucorner are allocated by this procedure and will point to an
- array of size dd->dim.
-
- Level: beginner
-
-.keywords: distributed array, refinement
-@*/
-PetscErrorCode DMADDAGetGhostCorners(DM dm, PetscInt **lcorner, PetscInt **ucorner)
-{
- DM_ADDA *dd = (DM_ADDA*)dm->data;
- PetscErrorCode ierr;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
- PetscValidPointer(lcorner,2);
- PetscValidPointer(ucorner,3);
- ierr = PetscMalloc1(dd->dim, lcorner);CHKERRQ(ierr);
- ierr = PetscMalloc1(dd->dim, ucorner);CHKERRQ(ierr);
- ierr = PetscMemcpy(*lcorner, dd->lgs, dd->dim*sizeof(PetscInt));CHKERRQ(ierr);
- ierr = PetscMemcpy(*ucorner, dd->lge, dd->dim*sizeof(PetscInt));CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDAMatSetValues"
-/*@C
- DMADDAMatSetValues - Inserts or adds a block of values into a matrix. The values
- are indexed geometrically with the help of the ADDA data structure.
- These values may be cached, so MatAssemblyBegin() and MatAssemblyEnd()
- MUST be called after all calls to ADDAMatSetValues() have been completed.
-
- Not Collective
-
- Input Parameters:
-+ mat - the matrix
-. addam - the ADDA geometry information for the rows
-. m - the number of rows
-. idxm - the row indices, each of the a proper ADDAIdx
-+ addan - the ADDA geometry information for the columns
-. n - the number of columns
-. idxn - the column indices, each of the a proper ADDAIdx
-. v - a logically two-dimensional array of values of size m*n
-- addv - either ADD_VALUES or INSERT_VALUES, where
- ADD_VALUES adds values to any existing entries, and
- INSERT_VALUES replaces existing entries with new values
-
- Notes:
- By default the values, v, are row-oriented and unsorted.
- See MatSetOption() for other options.
-
- Calls to ADDAMatSetValues() (and MatSetValues()) with the INSERT_VALUES and ADD_VALUES
- options cannot be mixed without intervening calls to the assembly
- routines.
-
- Efficiency Alert:
- The routine ADDAMatSetValuesBlocked() may offer much better efficiency
- for users of block sparse formats (MATSEQBAIJ and MATMPIBAIJ).
-
- Level: beginner
-
- Concepts: matrices^putting entries in
-
-.seealso: MatSetOption(), MatAssemblyBegin(), MatAssemblyEnd(), MatSetValues(), ADDAMatSetValuesBlocked(),
- InsertMode, INSERT_VALUES, ADD_VALUES
-@*/
-PetscErrorCode DMADDAMatSetValues(Mat mat, DM dmm, PetscInt m, const ADDAIdx idxm[],DM dmn, PetscInt n, const ADDAIdx idxn[],const PetscScalar v[], InsertMode addv)
-{
- DM_ADDA *ddm = (DM_ADDA*)dmm->data;
- DM_ADDA *ddn = (DM_ADDA*)dmn->data;
- PetscErrorCode ierr;
- PetscInt *nodemult;
- PetscInt i, j;
- PetscInt *matidxm, *matidxn;
- PetscInt *x, d;
- PetscInt idx;
-
- PetscFunctionBegin;
- /* find correct multiplying factors */
- ierr = PetscMalloc1(ddm->dim, &nodemult);CHKERRQ(ierr);
-
- nodemult[ddm->dim-1] = 1;
- for (j=ddm->dim-2; j>=0; j--) {
- nodemult[j] = nodemult[j+1]*(ddm->nodes[j+1]);
- }
- /* convert each coordinate in idxm to the matrix row index */
- ierr = PetscMalloc1(m, &matidxm);CHKERRQ(ierr);
- for (i=0; i<m; i++) {
- x = idxm[i].x; d = idxm[i].d;
- idx = 0;
- for (j=ddm->dim-1; j>=0; j--) {
- if (x[j] < 0) { /* "left", "below", etc. of boundary */
- if (ddm->periodic[j]) { /* periodic wraps around */
- x[j] += ddm->nodes[j];
- } else { /* non-periodic get discarded */
- matidxm[i] = -1; /* entries with -1 are ignored by MatSetValues() */
- goto endofloop_m;
- }
- }
- if (x[j] >= ddm->nodes[j]) { /* "right", "above", etc. of boundary */
- if (ddm->periodic[j]) { /* periodic wraps around */
- x[j] -= ddm->nodes[j];
- } else { /* non-periodic get discarded */
- matidxm[i] = -1; /* entries with -1 are ignored by MatSetValues() */
- goto endofloop_m;
- }
- }
- idx += x[j]*nodemult[j];
- }
- matidxm[i] = idx*(ddm->dof) + d;
- endofloop_m:
- ;
- }
- ierr = PetscFree(nodemult);CHKERRQ(ierr);
-
- /* find correct multiplying factors */
- ierr = PetscMalloc1(ddn->dim, &nodemult);CHKERRQ(ierr);
-
- nodemult[ddn->dim-1] = 1;
- for (j=ddn->dim-2; j>=0; j--) {
- nodemult[j] = nodemult[j+1]*(ddn->nodes[j+1]);
- }
- /* convert each coordinate in idxn to the matrix colum index */
- ierr = PetscMalloc1(n, &matidxn);CHKERRQ(ierr);
- for (i=0; i<n; i++) {
- x = idxn[i].x; d = idxn[i].d;
- idx = 0;
- for (j=ddn->dim-1; j>=0; j--) {
- if (x[j] < 0) { /* "left", "below", etc. of boundary */
- if (ddn->periodic[j]) { /* periodic wraps around */
- x[j] += ddn->nodes[j];
- } else { /* non-periodic get discarded */
- matidxn[i] = -1; /* entries with -1 are ignored by MatSetValues() */
- goto endofloop_n;
- }
- }
- if (x[j] >= ddn->nodes[j]) { /* "right", "above", etc. of boundary */
- if (ddn->periodic[j]) { /* periodic wraps around */
- x[j] -= ddn->nodes[j];
- } else { /* non-periodic get discarded */
- matidxn[i] = -1; /* entries with -1 are ignored by MatSetValues() */
- goto endofloop_n;
- }
- }
- idx += x[j]*nodemult[j];
- }
- matidxn[i] = idx*(ddn->dof) + d;
-endofloop_n:
- ;
- }
- /* call original MatSetValues() */
- ierr = MatSetValues(mat, m, matidxm, n, matidxn, v, addv);CHKERRQ(ierr);
- /* clean up */
- ierr = PetscFree(nodemult);CHKERRQ(ierr);
- ierr = PetscFree(matidxm);CHKERRQ(ierr);
- ierr = PetscFree(matidxn);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDASetParameters"
-PetscErrorCode DMADDASetParameters(DM dm,PetscInt dim, PetscInt *nodes,PetscInt *procs,PetscInt dof,PetscBool *periodic)
-{
- PetscErrorCode ierr;
- PetscMPIInt rank,size;
- MPI_Comm comm;
- PetscInt i;
- PetscInt nodes_total;
- PetscInt nodesleft;
- PetscInt procsleft;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
-
- PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
- ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
- ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
-
- /* total number of nodes */
- nodes_total = 1;
- for (i=0; i<dim; i++) nodes_total *= nodes[i];
- dd->dim = dim;
- dd->dof = dof;
- dd->periodic = periodic;
-
- ierr = PetscMalloc1(dim, &(dd->nodes));CHKERRQ(ierr);
- ierr = PetscMemcpy(dd->nodes, nodes, dim*sizeof(PetscInt));CHKERRQ(ierr);
-
- /* procs */
- ierr = PetscMalloc1(dim, &(dd->procs));CHKERRQ(ierr);
- /* create distribution of nodes to processors */
- if (procs == NULL) {
- procs = dd->procs;
- nodesleft = nodes_total;
- procsleft = size;
- /* figure out a good way to split the array to several processors */
- for (i=0; i<dim; i++) {
- if (i==dim-1) {
- procs[i] = procsleft;
- } else {
- /* calculate best partition */
- procs[i] = (PetscInt)(((PetscReal) nodes[i])*PetscPowReal(((PetscReal) procsleft)/((PetscReal) nodesleft),1./((PetscReal)(dim-i)))+0.5);
- if (procs[i]<1) procs[i]=1;
- while (procs[i] > 0) {
- if (procsleft % procs[i]) procs[i]--;
- else break;
- }
- nodesleft /= nodes[i];
- procsleft /= procs[i];
- }
- }
- } else {
- /* user provided the number of processors */
- ierr = PetscMemcpy(dd->procs, procs, dim*sizeof(PetscInt));CHKERRQ(ierr);
- }
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMSetUp_ADDA"
-PetscErrorCode DMSetUp_ADDA(DM dm)
-{
- PetscErrorCode ierr;
- PetscInt s=1; /* stencil width, fixed to 1 at the moment */
- PetscMPIInt rank,size;
- PetscInt i;
- PetscInt procsleft;
- PetscInt procsdimi;
- PetscInt ranki;
- PetscInt rpq;
- DM_ADDA *dd = (DM_ADDA*)dm->data;
- MPI_Comm comm;
- PetscInt *nodes,*procs,dim,dof;
- PetscBool *periodic;
-
- PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
- ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
- ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
- procs = dd->procs;
- nodes = dd->nodes;
- dim = dd->dim;
- dof = dd->dof;
- periodic = dd->periodic;
-
- /* check for validity */
- procsleft = 1;
- for (i=0; i<dim; i++) {
- if (nodes[i] < procs[i]) SETERRQ3(comm,PETSC_ERR_ARG_OUTOFRANGE,"Partition in direction %d is too fine! %D nodes, %D processors", i, nodes[i], procs[i]);
- procsleft *= procs[i];
- }
- if (procsleft != size) SETERRQ(comm,PETSC_ERR_PLIB, "Created or was provided with inconsistent distribution of processors");
-
-
- /* find out local region */
- ierr = PetscMalloc1(dim, &(dd->lcs));CHKERRQ(ierr);
- ierr = PetscMalloc1(dim, &(dd->lce));CHKERRQ(ierr);
- procsdimi = size;
- ranki = rank;
- for (i=0; i<dim; i++) {
- /* What is the number of processor for dimensions i+1, ..., dim-1? */
- procsdimi /= procs[i];
- /* these are all nodes that come before our region */
- rpq = ranki / procsdimi;
- dd->lcs[i] = rpq * (nodes[i]/procs[i]);
- if (rpq + 1 < procs[i]) {
- dd->lce[i] = (rpq + 1) * (nodes[i]/procs[i]);
- } else {
- /* last one gets all the rest */
- dd->lce[i] = nodes[i];
- }
- ranki = ranki - rpq*procsdimi;
- }
-
- /* compute local size */
- dd->lsize=1;
- for (i=0; i<dim; i++) {
- dd->lsize *= (dd->lce[i]-dd->lcs[i]);
- }
- dd->lsize *= dof;
-
- /* find out ghost points */
- ierr = PetscMalloc1(dim, &(dd->lgs));CHKERRQ(ierr);
- ierr = PetscMalloc1(dim, &(dd->lge));CHKERRQ(ierr);
- for (i=0; i<dim; i++) {
- if (periodic[i]) {
- dd->lgs[i] = dd->lcs[i] - s;
- dd->lge[i] = dd->lce[i] + s;
- } else {
- dd->lgs[i] = PetscMax(dd->lcs[i] - s, 0);
- dd->lge[i] = PetscMin(dd->lce[i] + s, nodes[i]);
- }
- }
-
- /* compute local size with ghost points */
- dd->lgsize=1;
- for (i=0; i<dim; i++) {
- dd->lgsize *= (dd->lge[i]-dd->lgs[i]);
- }
- dd->lgsize *= dof;
-
- /* create global and local prototype vector */
- ierr = VecCreateMPIWithArray(comm,dd->dof,dd->lsize,PETSC_DECIDE,0,&(dd->global));CHKERRQ(ierr);
-#if ADDA_NEEDS_LOCAL_VECTOR
- /* local includes ghost points */
- ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,dof,dd->lgsize,0,&(dd->local));CHKERRQ(ierr);
-#endif
-
- ierr = PetscMalloc1(dim, &(dd->refine));CHKERRQ(ierr);
- for (i=0; i<dim; i++) dd->refine[i] = 3;
- dd->dofrefine = 1;
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMCreate_ADDA"
-PETSC_EXTERN PetscErrorCode DMCreate_ADDA(DM dm)
-{
- PetscErrorCode ierr;
- DM_ADDA *dd;
-
- PetscFunctionBegin;
- ierr = PetscNewLog(dm,&dd);CHKERRQ(ierr);
- dm->data = (void*)dd;
-
- dm->ops->view = DMView;
- dm->ops->createglobalvector = DMCreateGlobalVector_ADDA;
- dm->ops->getcoloring = DMCreateColoring_ADDA;
- dm->ops->creatematrix = DMCreateMatrix_ADDA;
- dm->ops->createinterpolation = DMCreateInterpolation_ADDA;
- dm->ops->refine = DMRefine_ADDA;
- dm->ops->coarsen = DMCoarsen_ADDA;
- dm->ops->getinjection = DMCreateInjection_ADDA;
- dm->ops->getaggregates = DMCreateAggregates_ADDA;
- dm->ops->setup = DMSetUp_ADDA;
- dm->ops->destroy = DMDestroy_ADDA;
- PetscFunctionReturn(0);
-}
-
-
-#undef __FUNCT__
-#define __FUNCT__ "DMADDACreate"
-/*@C
- DMADDACreate - Creates and ADDA object that translate between coordinates
- in a geometric grid of arbitrary dimension and data in a PETSc vector
- distributed on several processors.
-
- Collective on MPI_Comm
-
- Input Parameters:
-+ comm - MPI communicator
-. dim - the dimension of the grid
-. nodes - array with d entries that give the number of nodes in each dimension
-. procs - array with d entries that give the number of processors in each dimension
- (or NULL if to be determined automatically)
-. dof - number of degrees of freedom per node
-- periodic - array with d entries that, i-th entry is set to true iff dimension i is periodic
-
- Output Parameters:
-. adda - pointer to ADDA data structure that is created
-
- Level: intermediate
-
-@*/
-PetscErrorCode DMADDACreate(MPI_Comm comm, PetscInt dim, PetscInt *nodes,PetscInt *procs,PetscInt dof, PetscBool *periodic,DM *dm_p)
-{
- PetscErrorCode ierr;
-
- PetscFunctionBegin;
- ierr = DMCreate(comm,dm_p);CHKERRQ(ierr);
- ierr = DMSetType(*dm_p,DMADDA);CHKERRQ(ierr);
- ierr = DMADDASetParameters(*dm_p,dim,nodes,procs,dof,periodic);CHKERRQ(ierr);
- ierr = DMSetUp(*dm_p);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
View
30 src/dm/impls/adda/addaimpl.h
@@ -1,30 +0,0 @@
-#if !defined(_ADDA_H)
-#define _ADDA_H
-
-#include <petscdmadda.h>
-#include <petsc-private/dmimpl.h>
-
-/* vector was allocated and never referenced, clearly some task was not finished */
-#define ADDA_HAS_LOCAL_VECTOR 0
-
-
-typedef struct {
- PetscInt dim; /* dimension of lattice */
- PetscInt dof; /* degrees of freedom per node */
- PetscInt *nodes; /* array of number of nodes in each dimension */
- PetscInt *procs; /* processor layout */
- PetscBool *periodic; /* true, if that dimension is periodic */
- PetscInt *lcs, *lce; /* corners of the locally stored portion of the grid */
- PetscInt *lgs, *lge; /* corners of the local portion of the grid
- including the ghost points */
- PetscInt lsize; /* number of nodes in local region */
- PetscInt lgsize; /* number of nodes in local region including ghost points */
- Vec global; /* global prototype vector */
-#if ADDA_HAS_LOCAL_VECTOR
- Vec local; /* local prototype vector */
-#endif
- PetscInt *refine; /* refinement factors for each dimension */
- PetscInt dofrefine; /* refinement factor for the dof */
-} DM_ADDA;
-
-#endif
View
11 src/dm/impls/adda/examples/makefile
@@ -1,11 +0,0 @@
-
-ALL: lib
-
-DIRS = tests
-LOCDIR = src/dm/impls/adda/examples/
-MANSEC = DM
-SUBMANSEC= DMADDA
-
-include ${PETSC_DIR}/conf/variables
-include ${PETSC_DIR}/conf/rules
-include ${PETSC_DIR}/conf/test
View
24 src/dm/impls/adda/examples/tests/ex1.c
@@ -1,24 +0,0 @@
-
-static char help[] = "Tests various ADDA routines.\n\n";
-
-#include <petscdmadda.h>
-
-#undef __FUNCT__
-#define __FUNCT__ "main"
-int main(int argc,char **argv)
-{
- PetscErrorCode ierr;
- ADDA adda;
- PetscInt nodes[4] = {20, 20, 10, 10};
- PetscBool periodic[4] = {PETSC_TRUE, PETSC_TRUE, PETSC_FALSE, PETSC_FALSE};
-
- ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);
-
- /* Create distributed array and get vectors */
- ierr = DMADDACreate(PETSC_COMM_WORLD, 4, nodes, NULL, 2, periodic, &adda);CHKERRQ(ierr);
-
- /* Free memory */
- ierr = DMDestroy(&adda);CHKERRQ(ierr);
- ierr = PetscFinalize();
- return 0;
-}
View
19 src/dm/impls/adda/examples/tests/makefile
@@ -1,19 +0,0 @@
-
-CFLAGS =
-FFLAGS =
-CPPFLAGS =
-FPPFLAGS =
-LOCDIR = src/dm/impls/adda/examples/tests/
-EXAMPLESC = ex1.c
-EXAMPLESF =
-MANSEC = DM
-SUBMANSEC = DMADDA
-
-include ${PETSC_DIR}/conf/variables
-include ${PETSC_DIR}/conf/rules
-
-ex1: ex1.o chkopts
- -${CLINKER} -o ex1 ex1.o ${PETSC_DM_LIB}
- ${RM} -f ex1.o
-
-include ${PETSC_DIR}/conf/test
View
16 src/dm/impls/adda/makefile
@@ -1,16 +0,0 @@
-
-ALL: lib
-
-CFLAGS =
-FFLAGS =
-SOURCEC = adda.c
-SOURCEF =
-LIBBASE = libpetscdm
-DIRS = examples
-LOCDIR = src/dm/impls/adda/
-MANSEC = DM
-SUBMANSEC= DMADDA
-
-include ${PETSC_DIR}/conf/variables
-include ${PETSC_DIR}/conf/rules
-include ${PETSC_DIR}/conf/test
View
2  src/dm/impls/makefile
@@ -1,6 +1,6 @@
ALL: lib
-DIRS = da adda sliced composite redundant plex shell patch moab circuit
+DIRS = da sliced composite redundant plex shell patch moab circuit
LOCDIR = src/dm/impls/
MANSEC = DM
View
2  src/dm/interface/dmregall.c
@@ -4,7 +4,6 @@ PETSC_EXTERN PetscErrorCode DMCreate_DA(DM);
PETSC_EXTERN PetscErrorCode DMCreate_Composite(DM);
PETSC_EXTERN PetscErrorCode DMCreate_Sliced(DM);
PETSC_EXTERN PetscErrorCode DMCreate_Shell(DM);
-PETSC_EXTERN PetscErrorCode DMCreate_ADDA(DM);
PETSC_EXTERN PetscErrorCode DMCreate_Redundant(DM);
PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM);
PETSC_EXTERN PetscErrorCode DMCreate_Patch(DM);
@@ -37,7 +36,6 @@ PetscErrorCode DMRegisterAll()
ierr = DMRegister(DMCOMPOSITE, DMCreate_Composite);CHKERRQ(ierr);
ierr = DMRegister(DMSLICED, DMCreate_Sliced);CHKERRQ(ierr);
ierr = DMRegister(DMSHELL, DMCreate_Shell);CHKERRQ(ierr);
- ierr = DMRegister(DMADDA, DMCreate_ADDA);CHKERRQ(ierr);
ierr = DMRegister(DMREDUNDANT, DMCreate_Redundant);CHKERRQ(ierr);
ierr = DMRegister(DMPLEX, DMCreate_Plex);CHKERRQ(ierr);
ierr = DMRegister(DMPATCH, DMCreate_Patch);CHKERRQ(ierr);
View
1  src/docs/website/miscellaneous/acknwldg.html
@@ -56,7 +56,6 @@
<li>Asbjorn Hoiland Aarrestad, (the explicit Runge-Kutta implementations)</li>
<li>G. Anciaux and J. Roman, (the interfaces to the partitioning packages PTScotch, Chaco, and Party)</li>
<li>Allison Baker, (the flexible GMRES and the LGMRES code)</li>
- <li>Arvid Bessen, (the adaptive smoothed aggregation algorithm (PCASA) and the ADDA object)</li>
<li>Chad Carroll, (the Win32 graphics),</li>
<li>Ethan Coon, (the PetscBag and many bug fixes),</li>
<li>Cameron Cooper, (portions of the VecScatter routines)</li>
View
109 src/ksp/ksp/examples/tutorials/ex38.c
@@ -1,109 +0,0 @@
-
-static char help[] = "Tests the aSA multigrid code.\n"
- "Parameters:\n"
- "-n n to use a matrix size of n\n";
-
-#include <petscdm.h>
-#include <petscdmda.h>
-#include <petscksp.h>
-#include <petscpcasa.h>
-
-PetscErrorCode Create1dLaplacian(PetscInt,Mat*);
-PetscErrorCode CalculateRhs(Vec);
-
-#undef __FUNCT__
-#define __FUNCT__ "main"
-int main(int Argc,char **Args)
-{
- PetscInt n = 60;
- PetscErrorCode ierr;
- Mat cmat;
- Vec b,x;
- KSP kspmg;
- PC pcmg;
- DM da;
-
- PetscInitialize(&Argc,&Args,(char*)0,help);
-
- ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
- ierr = Create1dLaplacian(n,&cmat);CHKERRQ(ierr);
- ierr = MatGetVecs(cmat,&x,0);CHKERRQ(ierr);
- ierr = MatGetVecs(cmat,&b,0);CHKERRQ(ierr);
- ierr = CalculateRhs(b);CHKERRQ(ierr);
- ierr = VecSet(x,0.0);CHKERRQ(ierr);
-
- ierr = KSPCreate(PETSC_COMM_WORLD,&kspmg);CHKERRQ(ierr);
- ierr = KSPSetType(kspmg, KSPCG);CHKERRQ(ierr);
-
- ierr = KSPGetPC(kspmg,&pcmg);CHKERRQ(ierr);
- ierr = PCSetType(pcmg,PCASA);CHKERRQ(ierr);
-
- /* maybe user wants to override some of the choices */
- ierr = KSPSetFromOptions(kspmg);CHKERRQ(ierr);
-
- ierr = KSPSetOperators(kspmg,cmat,cmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
-
- ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, n, 1, 1, 0, &da);CHKERRQ(ierr);
- ierr = DMDASetRefinementFactor(da, 3, 3, 3);CHKERRQ(ierr);
- ierr = PCSetDM(pcmg, (DM) da);CHKERRQ(ierr);
-
- ierr = PCASASetTolerances(pcmg, 1.e-10, 1.e-10, PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
- ierr = KSPSolve(kspmg,b,x);CHKERRQ(ierr);
- ierr = KSPDestroy(&kspmg);CHKERRQ(ierr);
- ierr = VecDestroy(&x);CHKERRQ(ierr);
- ierr = VecDestroy(&b);CHKERRQ(ierr);
- ierr = MatDestroy(&cmat);CHKERRQ(ierr);
- ierr = DMDestroy(&da);CHKERRQ(ierr);
- ierr = PetscFinalize();
- return 0;
-}
-
-/* --------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "Create1dLaplacian"
-PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
-{
- PetscScalar mone = -1.0,two = 2.0;
- PetscInt i,j,loc_start,loc_end;
- PetscErrorCode ierr;
-
- PetscFunctionBeginUser;
- ierr = MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n,PETSC_DECIDE, NULL, PETSC_DECIDE, NULL, mat);CHKERRQ(ierr);
-
- ierr = MatGetOwnershipRange(*mat,&loc_start,&loc_end);CHKERRQ(ierr);
- for (i=loc_start; i<loc_end; i++) {
- if (i>0) { j=i-1; ierr = MatSetValues(*mat,1,&i,1,&j,&mone,INSERT_VALUES);CHKERRQ(ierr); }
- ierr = MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);CHKERRQ(ierr);
- if (i<n-1) { j=i+1; ierr = MatSetValues(*mat,1,&i,1,&j,&mone,INSERT_VALUES);CHKERRQ(ierr); }
- }
- ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
- ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-/* --------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "CalculateRhs"
-PetscErrorCode CalculateRhs(Vec u)
-{
- PetscErrorCode ierr;
- PetscInt i,n,loc_start,loc_end;
- PetscReal h;
- PetscScalar uu;
-
- PetscFunctionBeginUser;
- ierr = VecGetSize(u,&n);CHKERRQ(ierr);
- ierr = VecGetOwnershipRange(u,&loc_start,&loc_end);CHKERRQ(ierr);
- h = 1.0/((PetscReal)(n+1));
- uu = 2.0*h*h;
- for (i=loc_start; i<loc_end; i++) {
- ierr = VecSetValues(u,1,&i,&uu,INSERT_VALUES);CHKERRQ(ierr);
- }
- ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
- ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
View
256 src/ksp/ksp/examples/tutorials/ex39.c
@@ -1,256 +0,0 @@
-
-static char help[] = "Lattice Gauge 2D model.\n"
- "Parameters:\n"
- "-size n to use a grid size of n, i.e n space and n time steps\n"
- "-beta b controls the randomness of the gauge field\n"
- "-rho r the quark mass (?)";
-
-#include <petscksp.h>
-#include <petscpcasa.h>
-#include <petscdm.h>
-#include <petscdmda.h>
-
-PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig);
-
-#undef __FUNCT__
-#define __FUNCT__ "main"
-int main(int Argc,char **Args)
-{
- PetscBool flg;
- PetscInt n = -6;
- PetscScalar rho = 1.0;
- PetscReal h;
- PetscReal beta = 1.0;
- DM da;
- PetscRandom rctx;
- PetscMPIInt comm_size;
- Mat H,HtH;
- PetscInt x, y, xs, ys, xm, ym;
- PetscReal r1, r2;
- PetscScalar uxy1, uxy2;
- MatStencil sxy, sxy_m;
- PetscScalar val, valconj;
- Vec b, Htb,xvec;
- KSP kspmg;
- PC pcmg;
- PetscErrorCode ierr;
- PetscInt ix[1] = {0};
- PetscScalar vals[1] = {1.0};
-
- PetscInitialize(&Argc,&Args,(char*)0,help);
- ierr = PetscOptionsGetInt(NULL,"-size",&n,&flg);CHKERRQ(ierr);
- ierr = PetscOptionsGetReal(NULL,"-beta",&beta,&flg);CHKERRQ(ierr);
- ierr = PetscOptionsGetScalar(NULL,"-rho",&rho,&flg);CHKERRQ(ierr);
-
- /* Set the fudge parameters, we scale the whole thing by 1/(2*h) later */
- h = 1.;
- rho *= 1./(2.*h);
-
- /* Geometry info */
- ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, n, n,
- PETSC_DECIDE, PETSC_DECIDE, 2 /* this is the # of dof's */,
- 1, NULL, NULL, &da);CHKERRQ(ierr);
-
- /* Random numbers */
- ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
- ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
-
- /* Single or multi processor ? */
- ierr = MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);CHKERRQ(ierr);
-
- /* construct matrix */
- ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);
- ierr = DMCreateMatrix(da, &H);CHKERRQ(ierr);
-
- /* get local corners for this processor */
- ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
-
- /* Assemble the matrix */
- for (x=xs; x<xs+xm; x++) {
- for (y=ys; y<ys+ym; y++) {
- /* each lattice point sets only the *forward* pointing parameters (right, down),
- i.e. Nabla_1^+ and Nabla_2^+.
- In this way we can use only local random number creation. That means
- we also have to set the corresponding backward pointing entries. */
- /* Compute some normally distributed random numbers via Box-Muller */
- ierr = PetscRandomGetValueReal(rctx, &r1);CHKERRQ(ierr);
- r1 = 1.-r1; /* to change from [0,1) to (0,1], which we need for the log */
- ierr = PetscRandomGetValueReal(rctx, &r2);CHKERRQ(ierr);
- PetscReal R = PetscSqrtReal(-2.*PetscLogReal(r1));
- PetscReal c = PetscCosReal(2.*PETSC_PI*r2);
- PetscReal s = PetscSinReal(2.*PETSC_PI*r2);
-
- /* use those to set the field */
- uxy1 = PetscExpScalar(((PetscScalar) (R*c/beta))*PETSC_i);
- uxy2 = PetscExpScalar(((PetscScalar) (R*s/beta))*PETSC_i);
-
- sxy.i = x; sxy.j = y; /* the point where we are */
-
- /* center action */
- sxy.c = 0; /* spin 0, 0 */
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy, &rho, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 1; /* spin 1, 1 */
- val = -rho;
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
-
- sxy_m.i = x+1; sxy_m.j = y; /* right action */
- sxy.c = 0; sxy_m.c = 0; /* spin 0, 0 */
- val = -uxy1; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 0; sxy_m.c = 1; /* spin 0, 1 */
- val = -uxy1; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
- val = uxy1; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
- val = uxy1; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy_m.i = x; sxy_m.j = y+1; /* down action */
- sxy.c = 0; sxy_m.c = 0; /* spin 0, 0 */
- val = -uxy2; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 0; sxy_m.c = 1; /* spin 0, 1 */
- val = -PETSC_i*uxy2; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 1; sxy_m.c = 0; /* spin 1, 0 */
- val = -PETSC_i*uxy2; valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- sxy.c = 1; sxy_m.c = 1; /* spin 1, 1 */
- val = PetscConj(uxy2); valconj = PetscConj(val);
- ierr = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- }
- }
-
- ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
- ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-
- /* scale H */
- ierr = MatScale(H, 1./(2.*h));CHKERRQ(ierr);
-
- /* it looks like H is Hermetian */
- /* construct normal equations */
- ierr = MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);CHKERRQ(ierr);
-
- /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
-/* Mat perm; */
-/* ierr = DMCreateMatrix(da, &perm);CHKERRQ(ierr); */
-/* PetscInt row, col; */
-/* PetscScalar one = 1.0; */
-/* for (PetscInt i=0; i<n; i++) { */
-/* for (PetscInt j=0; j<n; j++) { */
-/* row = (i*n+j)*2; col = i*n+j; */
-/* ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
-/* row = (i*n+j)*2+1; col = i*n+j + n*n; */
-/* ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
-/* } */
-/* } */
-/* ierr = MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
-/* ierr = MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
-
-/* Mat Hperm; */
-/* ierr = MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm);CHKERRQ(ierr); */
-/* ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n");CHKERRQ(ierr); */
-/* ierr = MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
-
-/* Mat HtHperm; */
-/* ierr = MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm);CHKERRQ(ierr); */
-/* ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n");CHKERRQ(ierr); */
-/* ierr = MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
-
- /* right hand side */
- ierr = DMCreateGlobalVector(da, &b);CHKERRQ(ierr);
- ierr = VecSet(b,0.0);CHKERRQ(ierr);
- ierr = VecSetValues(b, 1, ix, vals, INSERT_VALUES);CHKERRQ(ierr);
- ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
- ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
-/* ierr = VecSetRandom(b, rctx);CHKERRQ(ierr); */
- ierr = VecDuplicate(b, &Htb);CHKERRQ(ierr);
- ierr = MatMultTranspose(H, b, Htb);CHKERRQ(ierr);
-
- /* construct solver */
- ierr = KSPCreate(PETSC_COMM_WORLD,&kspmg);CHKERRQ(ierr);
- ierr = KSPSetType(kspmg, KSPCG);CHKERRQ(ierr);
-
- ierr = KSPGetPC(kspmg,&pcmg);CHKERRQ(ierr);
- ierr = PCSetType(pcmg,PCASA);CHKERRQ(ierr);
-
- /* maybe user wants to override some of the choices */
- ierr = KSPSetFromOptions(kspmg);CHKERRQ(ierr);
-
- ierr = KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
-
- ierr = DMDASetRefinementFactor(da, 3, 3, 3);CHKERRQ(ierr);
- ierr = PCSetDM(pcmg,da);CHKERRQ(ierr);
-
- ierr = PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
-
- ierr = VecDuplicate(b, &xvec);CHKERRQ(ierr);
- ierr = VecSet(xvec, 0.0);CHKERRQ(ierr);
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
- ierr = KSPSolve(kspmg, Htb, xvec);CHKERRQ(ierr);
-
-/* ierr = VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
-
- ierr = KSPDestroy(&kspmg);CHKERRQ(ierr);
- ierr = VecDestroy(&xvec);CHKERRQ(ierr);
-
- /* seems to be destroyed by KSPDestroy */
- ierr = VecDestroy(&b);CHKERRQ(ierr);
- ierr = VecDestroy(&Htb);CHKERRQ(ierr);
- ierr = MatDestroy(&HtH);CHKERRQ(ierr);
- ierr = MatDestroy(&H);CHKERRQ(ierr);
-
- ierr = DMDestroy(&da);CHKERRQ(ierr);
- ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
- ierr = PetscFinalize();
- return 0;
-}
-
-/* --------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "computeMaxEigVal"
-PetscErrorCode computeMaxEigVal(Mat A, PetscInt its, PetscScalar *eig)
-{
- PetscErrorCode ierr;
- PetscRandom rctx; /* random number generator context */
- Vec x0, x, x_1, tmp;
- PetscScalar lambda_its, lambda_its_1;
- PetscInt i;
-
- PetscFunctionBeginUser;
- ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
- ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
- ierr = MatGetVecs(A, &x_1, &x);CHKERRQ(ierr);
- ierr = VecSetRandom(x, rctx);CHKERRQ(ierr);
- ierr = VecDuplicate(x, &x0);CHKERRQ(ierr);
- ierr = VecCopy(x, x0);CHKERRQ(ierr);
-
- ierr = MatMult(A, x, x_1);CHKERRQ(ierr);
- for (i=0; i<its; i++) {
- tmp = x; x = x_1; x_1 = tmp;
- ierr = MatMult(A, x, x_1);CHKERRQ(ierr);
- }
- ierr = VecDot(x0, x, &lambda_its);CHKERRQ(ierr);
- ierr = VecDot(x0, x_1, &lambda_its_1);CHKERRQ(ierr);
-
- *eig = lambda_its_1/lambda_its;
-
- ierr = VecDestroy(&x0);CHKERRQ(ierr);
- ierr = VecDestroy(&x);CHKERRQ(ierr);
- ierr = VecDestroy(&x_1);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
View
302 src/ksp/ksp/examples/tutorials/ex40.c
@@ -1,302 +0,0 @@
-
-static char help[] = "Lattice Gauge 2D model.\n"
- "Parameters:\n"
- "-size n to use a grid size of n, i.e n space and n time steps\n"
- "-beta b controls the randomness of the gauge field\n"
- "-rho r the quark mass (?)";
-
-#include <petscksp.h>
-#include <petscpcasa.h>
-#include <petscdm.h>
-#include <petscdmadda.h>
-
-PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig);
-
-#undef __FUNCT__
-#define __FUNCT__ "main"
-int main(int Argc,char **Args)
-{
- PetscBool flg;
- PetscInt n = 6,i;
- PetscScalar rho = 1.0;
- PetscReal h;
- PetscReal beta = 1.0;
- DM adda;
- PetscInt nodes[2];
- PetscBool periodic[2];
- PetscInt refine[2];
- PetscRandom rctx;
- PetscMPIInt comm_size;
- Mat H;
- PetscInt *lcs, *lce;
- PetscInt x, y;
- PetscReal r1, r2;
- PetscScalar uxy1, uxy2;
- ADDAIdx sxy, sxy_m;
- PetscScalar val, valconj;
- Mat HtH;
- Vec b, Htb;
- Vec xvec;
- KSP kspmg;
- PC pcmg;
- PetscErrorCode ierr;
-
- PetscInitialize(&Argc,&Args,(char*)0,help);
- ierr = PetscOptionsGetInt(NULL,"-size",&n,&flg);CHKERRQ(ierr);
- ierr = PetscOptionsGetReal(NULL,"-beta",&beta,&flg);CHKERRQ(ierr);
- ierr = PetscOptionsGetScalar(NULL,"-rho",&rho,&flg);CHKERRQ(ierr);
-
- /* Set the fudge parameters, we scale the whole thing by 1/(2*h) later */
- h = 1.;
- rho *= 1./(2.*h);
-
- /* Geometry info */
- for (i=0; i<2; i++) {
- nodes[i] = n;
- periodic[i] = PETSC_TRUE;
- refine[i] = 3;
- }
- ierr = DMADDACreate(PETSC_COMM_WORLD, 2, nodes, NULL, 2,periodic, &adda);CHKERRQ(ierr);
- ierr = DMADDASetRefinement(adda, refine, 2);CHKERRQ(ierr);
-
- /* Random numbers */
- ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
- ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
-
- /* Single or multi processor ? */
- ierr = MPI_Comm_size(PETSC_COMM_WORLD,&comm_size);CHKERRQ(ierr);
-
- /* construct matrix */
- ierr = DMSetMatType(adda,MATAIJ);CHKERRQ(ierr);
- ierr = DMCreateMatrix(adda, &H);CHKERRQ(ierr);
-
- /* get local corners for this processor, user is responsible for freeing lcs,lce */
- ierr = DMADDAGetCorners(adda, &lcs, &lce);CHKERRQ(ierr);
-
- /* Allocate space for the indices that we use to construct the matrix */
- ierr = PetscMalloc1(2, &(sxy.x));CHKERRQ(ierr);
- ierr = PetscMalloc1(2, &(sxy_m.x));CHKERRQ(ierr);
-
- /* Assemble the matrix */
- for (x=lcs[0]; x<lce[0]; x++) {
- for (y=lcs[1]; y<lce[1]; y++) {
- /* each lattice point sets only the *forward* pointing parameters (right, down),
- i.e. Nabla_1^+ and Nabla_2^+.
- In this way we can use only local random number creation. That means
- we also have to set the corresponding backward pointing entries. */
- /* Compute some normally distributed random numbers via Box-Muller */
- ierr = PetscRandomGetValueReal(rctx, &r1);CHKERRQ(ierr);
- r1 = 1.-r1; /* to change from [0,1) to (0,1], which we need for the log */
- ierr = PetscRandomGetValueReal(rctx, &r2);CHKERRQ(ierr);
- PetscReal R = PetscSqrtReal(-2.*log(r1));
- PetscReal c = PetscCosReal(2.*PETSC_PI*r2);
- PetscReal s = PetscSinReal(2.*PETSC_PI*r2);
-
- /* use those to set the field */
- uxy1 = PetscExpScalar(((PetscScalar) (R*c/beta))*PETSC_i);
- uxy2 = PetscExpScalar(((PetscScalar) (R*s/beta))*PETSC_i);
-
- sxy.x[0] = x; sxy.x[1] = y; /* the point where we are */
-
- /* center action */
- sxy.d = 0; /* spin 0, 0 */
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy, &rho, ADD_VALUES);CHKERRQ(ierr);
- sxy.d = 1; /* spin 1, 1 */
- val = -rho;
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
-
- sxy_m.x[0] = x+1; sxy_m.x[1] = y; /* right action */
- sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
- val = -uxy1; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
- val = -uxy1; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
- val = uxy1; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
- val = uxy1; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy_m.x[0] = x; sxy_m.x[1] = y+1; /* down action */
- sxy.d = 0; sxy_m.d = 0; /* spin 0, 0 */
- val = -uxy2; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy.d = 0; sxy_m.d = 1; /* spin 0, 1 */
- val = -PETSC_i*uxy2; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy.d = 1; sxy_m.d = 0; /* spin 1, 0 */
- val = -PETSC_i*uxy2; valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
-
- sxy.d = 1; sxy_m.d = 1; /* spin 1, 1 */
- val = PetscConj(uxy2); valconj = PetscConj(val);
-
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy_m, adda, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
- ierr = DMADDAMatSetValues(H, adda, 1, &sxy, adda, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
- }
- }
-
- ierr = PetscFree(sxy.x);CHKERRQ(ierr);
- ierr = PetscFree(sxy_m.x);CHKERRQ(ierr);
-
- ierr = PetscFree(lcs);CHKERRQ(ierr);
- ierr = PetscFree(lce);CHKERRQ(ierr);
-
- ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
- ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
-
- /* scale H */
- ierr = MatScale(H, 1./(2.*h));CHKERRQ(ierr);
-
- /* construct normal equations */
- ierr = MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);CHKERRQ(ierr);
-
- PetscScalar mineval;
- ierr = computeMinEigVal(HtH, 1000, &mineval);CHKERRQ(ierr);
- ierr = PetscPrintf(PETSC_COMM_WORLD, "Minimum eigenvalue of H^{dag} H is %f\n", (double)PetscAbsScalar(mineval));CHKERRQ(ierr);
-
- /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
-/* Mat perm; */
-/* ierr = ADDACreatematrix(adda, MATSEQAIJ, &perm);CHKERRQ(ierr); */
-/* PetscInt row, col; */
-/* PetscScalar one = 1.0; */
-/* for (PetscInt i=0; i<n; i++) { */
-/* for (PetscInt j=0; j<n; j++) { */
-/* row = (i*n+j)*2; col = i*n+j; */
-/* ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
-/* row = (i*n+j)*2+1; col = i*n+j + n*n; */
-/* ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
-/* } */
-/* } */
-/* ierr = MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
-/* ierr = MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
-
-/* Mat Hperm; */
-/* ierr = MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm);CHKERRQ(ierr); */
-/* ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n");CHKERRQ(ierr); */
-/* ierr = MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
-
-/* Mat HtHperm; */
-/* ierr = MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm);CHKERRQ(ierr); */
-/* ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n");CHKERRQ(ierr); */
-/* ierr = MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
-
- /* right hand side */
- ierr = DMCreateGlobalVector(adda, &b);CHKERRQ(ierr);
- ierr = VecSet(b,0.0);CHKERRQ(ierr);
- PetscInt ix[1] = {0};
- PetscScalar vals[1] = {1.0};
- ierr = VecSetValues(b, 1, ix, vals, INSERT_VALUES);CHKERRQ(ierr);
- ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
- ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
-/* ierr = VecSetRandom(b, rctx);CHKERRQ(ierr); */
- ierr = VecDuplicate(b, &Htb);CHKERRQ(ierr);
- ierr = MatMultTranspose(H, b, Htb);CHKERRQ(ierr);
-
- /* construct solver */
- ierr = KSPCreate(PETSC_COMM_WORLD,&kspmg);CHKERRQ(ierr);
- ierr = KSPSetType(kspmg, KSPCG);CHKERRQ(ierr);
-
- ierr = KSPGetPC(kspmg,&pcmg);CHKERRQ(ierr);
- ierr = PCSetType(pcmg,PCASA);CHKERRQ(ierr);
-
- /* maybe user wants to override some of the choices */
- ierr = KSPSetFromOptions(kspmg);CHKERRQ(ierr);
-
- ierr = KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
-
- ierr = PCSetDM(pcmg,adda);CHKERRQ(ierr);
- ierr = DMDestroy(&adda);CHKERRQ(ierr);
-
- ierr = PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
-
- ierr = VecDuplicate(b, &xvec);CHKERRQ(ierr);
- ierr = VecSet(xvec, 0.0);CHKERRQ(ierr);
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
- ierr = KSPSolve(kspmg, Htb, xvec);CHKERRQ(ierr);
-
-/* ierr = VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */
-
- ierr = KSPDestroy(&kspmg);CHKERRQ(ierr);
-
- ierr = VecDestroy(&xvec);CHKERRQ(ierr);
- ierr = VecDestroy(&b);CHKERRQ(ierr);
- ierr = VecDestroy(&Htb);CHKERRQ(ierr);
- ierr = MatDestroy(&H);CHKERRQ(ierr);
- ierr = MatDestroy(&HtH);CHKERRQ(ierr);
-
- ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
- ierr = PetscFinalize();
- return 0;
-}
-
-/* --------------------------------------------------------------------- */
-#undef __FUNCT__
-#define __FUNCT__ "computeMinEigVal"
-PetscErrorCode computeMinEigVal(Mat A, PetscInt its, PetscScalar *eig)
-{
- PetscErrorCode ierr;
- PetscRandom rctx; /* random number generator context */
- Vec x0, x, x_1, tmp;
- PetscScalar lambda_its, lambda_its_1;
- PetscReal norm;
- Mat G;
- PetscInt i;
-
- PetscFunctionBeginUser;
- ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
- ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
-
- /* compute G = I-1/norm(A)*A */
- ierr = MatNorm(A, NORM_1, &norm);CHKERRQ(ierr);
- ierr = MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &G);CHKERRQ(ierr);
- ierr = MatShift(G, -norm);CHKERRQ(ierr);
- ierr = MatScale(G, -1./norm);CHKERRQ(ierr);
-
- ierr = MatGetVecs(G, &x_1, &x);CHKERRQ(ierr);
- ierr = VecSetRandom(x, rctx);CHKERRQ(ierr);
- ierr = VecDuplicate(x, &x0);CHKERRQ(ierr);
- ierr = VecCopy(x, x0);CHKERRQ(ierr);
-
- ierr = MatMult(G, x, x_1);CHKERRQ(ierr);
- for (i=0; i<its; i++) {
- tmp = x; x = x_1; x_1 = tmp;
- ierr = MatMult(G, x, x_1);CHKERRQ(ierr);
- }
- ierr = VecDot(x0, x, &lambda_its);CHKERRQ(ierr);
- ierr = VecDot(x0, x_1, &lambda_its_1);CHKERRQ(ierr);
-
- *eig = norm*(1.-lambda_its_1/lambda_its);
-
- ierr = VecDestroy(&x0);CHKERRQ(ierr);
- ierr = VecDestroy(&x);CHKERRQ(ierr);
- ierr = VecDestroy(&x_1);CHKERRQ(ierr);
- ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
- ierr = MatDestroy(&G);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
View
24 src/ksp/ksp/examples/tutorials/makefile
@@ -7,7 +7,7 @@ LOCDIR = src/ksp/ksp/examples/tutorials/
EXAMPLESC = ex1.c ex2.c ex3.c ex4.c ex5.c ex7.c ex8.c ex8g.c ex9.c \
ex10.c ex11.c ex12.c ex13.c ex15.c ex16.c ex18.c ex23.c \
ex25.c ex27.c ex28.c ex29.c ex30.c ex31.c ex32.c ex34.c \
- ex38.c ex39.c ex40.c ex41.c ex42.c ex43.c \
+ ex41.c ex42.c ex43.c \
ex45.c ex46.c ex49.c ex50.c ex51.c ex52.c ex53.c \
ex54.c ex55.c ex56.c ex58.c
EXAMPLESF = ex1f.F ex2f.F ex6f.F ex11f.F ex13f90.F ex14f.F ex15f.F ex21f.F ex22f.F ex44f.F90 ex45f.F ex54f.F
@@ -922,24 +922,6 @@ runex34:
else printf "${PWD}\nPossible problem with with ex34_1, diffs above\n=========================================\n"; fi; \
${RM} -f ex34_1.tmp
-runex38:
- -@${MPIEXEC} -n 1 ./ex38 -ksp_monitor_short > ex38_1.tmp 2>&1; \
- if (${DIFF} output/ex38_1.out ex38_1.tmp) then true; \
- else printf "${PWD}\nPossible problem with with ex38_1, diffs above\n=========================================\n"; fi; \
- ${RM} -f ex38_1.tmp
-
-runex39:
- -@${MPIEXEC} -n 1 ./ex39 -mat_no_inode -ksp_monitor_short > ex39_1.tmp 2>&1; \
- if (${DIFF} output/ex39_1.out ex39_1.tmp) then true; \
- else printf "${PWD}\nPossible problem with with ex39_1, diffs above\n=========================================\n"; fi; \
- ${RM} -f ex39_1.tmp
-
-runex40:
- -@${MPIEXEC} -n 1 ./ex40 -mat_no_inode -ksp_monitor_short > ex40_1.tmp 2>&1; \
- if (${DIFF} output/ex40_1.out ex40_1.tmp) then true; \
- else printf "${PWD}\nPossible problem with with ex40_1, diffs above\n=========================================\n"; fi; \
- ${RM} -f ex40_1.tmp
-
runex43:
-@${MPIEXEC} -n 1 ./ex43 -stokes_ksp_type fgmres -stokes_ksp_rtol 1e-8 -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short > ex43_1.tmp 2>&1; \
${DIFF} output/ex43_1.out ex43_1.tmp || printf "${PWD}\nPossible problem with with ex43_1, diffs above\n=========================================\n"; \
@@ -1132,7 +1114,7 @@ TESTEXAMPLES_C = ex1.PETSc runex1 runex1_2 runex1_3 ex1.rm ex2.PETSc run
ex15.PETSc runex15 ex15.rm ex16.PETSc runex16 ex16.rm \
ex23.PETSc runex23 runex23_2 ex23.rm ex25.PETSc runex25 runex25_2 ex25.rm \
ex27.PETSc ex27.rm ex28.PETSc ex28.rm ex29.PETSc ex29.rm \
- ex31.PETSc ex31.rm ex32.PETSc runex32 ex32.rm ex34.PETSc runex34 ex34.rm ex38.PETSc runex38 ex38.rm \
+ ex31.PETSc ex31.rm ex32.PETSc runex32 ex32.rm ex34.PETSc runex34 ex34.rm \
ex43.PETSc runex43 runex43_2 runex43_3 runex43_bjacobi runex43_bjacobi_baij ex43.rm \
ex45.PETSc runex45 runex45_2 ex45.rm \
ex49.PETSc runex49 runex49_2 runex49_3 runex49_5 ex49.rm ex53.PETSc runex53 ex53.rm ex55.PETSc runex55_Classical ex55.rm\
@@ -1146,7 +1128,7 @@ TESTEXAMPLES_FORTRAN_MPIUNI = ex1f.PETSc runex1f ex1f.rm ex6f.PETSc runex6f e
TESTEXAMPLES_C_X_MPIUNI = ex1.PETSc runex1 runex1_2 runex1_3 ex1.rm ex2.PETSc runex2 runex2_3 ex2.rm \
ex7.PETSc ex7.rm ex5.PETSc ex5.rm ex9.PETSc runex9 ex9.rm \
ex23.PETSc runex23 ex23.rm
-TESTEXAMPLES_C_COMPLEX = ex10.PETSc ex10.rm ex11.PETSc runex11 ex11.rm ex39.PETSc runex39 ex39.rm ex40.PETSc runex40 ex40.rm
+TESTEXAMPLES_C_COMPLEX = ex10.PETSc ex10.rm ex11.PETSc runex11 ex11.rm
TESTEXAMPLES_DATAFILESPATH = ex10.PETSc runex10_2 runex10_3 runex10_4 runex10_5 runex10_6 runex10_7 runex10_8 \
runex10_9 runex10_10 runex10_19 runex10_ILU runex10_ILUBAIJ runex10_cg \
runex10_cg_singlereduction runex10_seqaijcrl runex10_mpiaijcrl runex10_seqaijperm runex10_mpiaijperm ex10.rm
View
2,047 src/ksp/pc/impls/asa/asa.c
@@ -1,2047 +0,0 @@
-
-/* --------------------------------------------------------------------
-
- Contributed by Arvid Bessen, Columbia University, June 2007
-
- This file implements a ASA preconditioner in PETSc as part of PC.
-
- The adaptive smoothed aggregation algorithm is described in the paper
- "Adaptive Smoothed Aggregation (ASA)", M. Brezina, R. Falgout, S. MacLachlan,
- T. Manteuffel, S. McCormick, and J. Ruge, SIAM Journal on Scientific Computing,
- SISC Volume 25 Issue 6, Pages 1896-1920.
-
- For an example usage of this preconditioner, see, e.g.
- $PETSC_DIR/src/ksp/ksp/examples/tutorials/ex38.c ex39.c
- and other files in that directory.
-
- This code is still somewhat experimental. A number of improvements would be
- - keep vectors allocated on each level, instead of destroying them
- (see mainly PCApplyVcycleOnLevel_ASA)
- - in PCCreateTransferOp_ASA we get all of the submatrices at once, this could
- be optimized by differentiating between local and global matrices
- - the code does not handle it gracefully if there is just one level
- - if relaxation is sufficient, exit of PCInitializationStage_ASA is not
- completely clean
- - default values could be more reasonable, especially for parallel solves,
- where we need a parallel LU or similar
- - the richardson scaling parameter is somewhat special, should be treated in a
- good default way
- - a number of parameters for smoother (sor_omega, etc.) that we store explicitly
- could be kept in the respective smoothers themselves
- - some parameters have to be set via command line options, there are no direct
- function calls available
- - numerous other stuff
-
- Example runs in parallel would be with parameters like
- mpiexec ./program -pc_asa_coarse_pc_factor_mat_solver_package mumps -pc_asa_direct_solver 200
- -pc_asa_max_cand_vecs 4 -pc_asa_mu_initial 50 -pc_asa_richardson_scale 1.0
- -pc_asa_rq_improve 0.9 -asa_smoother_pc_type asm -asa_smoother_sub_pc_type sor
-
- -------------------------------------------------------------------- */
-
-/*
- This defines the data structures for the smoothed aggregation procedure
-*/
-#include <../src/ksp/pc/impls/asa/asa.h>
-#include <petscblaslapack.h>
-
-/* -------------------------------------------------------------------------- */
-
-/* Event logging */
-
-PetscLogEvent PC_InitializationStage_ASA, PC_GeneralSetupStage_ASA;
-PetscLogEvent PC_CreateTransferOp_ASA, PC_CreateVcycle_ASA;
-PetscBool asa_events_registered = PETSC_FALSE;
-
-#undef __FUNCT__
-#define __FUNCT__ "PCASASetTolerances"
-/*@C
- PCASASetTolerances - Sets the convergence thresholds for ASA algorithm
-
- Collective on PC
-
- Input Parameter:
-+ pc - the context
-. rtol - the relative convergence tolerance
- (relative decrease in the residual norm)
-. abstol - the absolute convergence tolerance
- (absolute size of the residual norm)
-. dtol - the divergence tolerance
- (amount residual can increase before KSPConvergedDefault()
- concludes that the method is diverging)
-- maxits - maximum number of iterations to use
-
- Notes:
- Use PETSC_DEFAULT to retain the default value of any of the tolerances.
-
- Level: advanced
-@*/
-PetscErrorCode PCASASetTolerances(PC pc, PetscReal rtol, PetscReal abstol,PetscReal dtol, PetscInt maxits)
-{
- PetscErrorCode ierr;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(pc,PC_CLASSID,1);
- ierr = PetscTryMethod(pc,"PCASASetTolerances_C",(PC,PetscReal,PetscReal,PetscReal,PetscInt),(pc,rtol,abstol,dtol,maxits));CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCASASetTolerances_ASA"
-static PetscErrorCode PCASASetTolerances_ASA(PC pc, PetscReal rtol, PetscReal abstol,PetscReal dtol, PetscInt maxits)
-{
- PC_ASA *asa = (PC_ASA*) pc->data;
-
- PetscFunctionBegin;
- PetscValidHeaderSpecific(pc,PC_CLASSID,1);
- if (rtol != PETSC_DEFAULT) asa->rtol = rtol;
- if (abstol != PETSC_DEFAULT) asa->abstol = abstol;
- if (dtol != PETSC_DEFAULT) asa->divtol = dtol;
- if (maxits != PETSC_DEFAULT) asa->max_it = maxits;
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PCCreateLevel_ASA"
-/*
- PCCreateLevel_ASA - Creates one level for the ASA algorithm
-
- Input Parameters:
-+ level - current level
-. comm - MPI communicator object
-. next - pointer to next level
-. prev - pointer to previous level
-. ksptype - the KSP type for the smoothers on this level
-- pctype - the PC type for the smoothers on this level
-
- Output Parameters:
-. new_asa_lev - the newly created level
-
-.keywords: ASA, create, levels, multigrid
-*/
-PetscErrorCode PCCreateLevel_ASA(PC_ASA_level **new_asa_lev, int level,MPI_Comm comm, PC_ASA_level *prev,PC_ASA_level *next,KSPType ksptype, PCType pctype)
-{
- PetscErrorCode ierr;
- PC_ASA_level *asa_lev;
-
- PetscFunctionBegin;
- ierr = PetscMalloc(sizeof(PC_ASA_level), &asa_lev);CHKERRQ(ierr);
-
- asa_lev->level = level;
- asa_lev->size = 0;
-
- asa_lev->A = 0;
- asa_lev->B = 0;
- asa_lev->x = 0;
- asa_lev->b = 0;
- asa_lev->r = 0;
-
- asa_lev->dm = 0;
- asa_lev->aggnum = 0;
- asa_lev->agg = 0;
- asa_lev->loc_agg_dofs = 0;
- asa_lev->agg_corr = 0;
- asa_lev->bridge_corr = 0;
-
- asa_lev->P = 0;
- asa_lev->Pt = 0;
- asa_lev->smP = 0;
- asa_lev->smPt = 0;
-
- asa_lev->comm = comm;
-
- asa_lev->smoothd = 0;
- asa_lev->smoothu = 0;
-
- asa_lev->prev = prev;
- asa_lev->next = next;
-
- *new_asa_lev = asa_lev;
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PrintResNorm"
-PetscErrorCode PrintResNorm(Mat A, Vec x, Vec b, Vec r)
-{
- PetscErrorCode ierr;
- PetscBool destroyr = PETSC_FALSE;
- PetscReal resnorm;
- MPI_Comm Acomm;
-
- PetscFunctionBegin;
- if (!r) {
- ierr = MatGetVecs(A, NULL, &r);CHKERRQ(ierr);
- destroyr = PETSC_TRUE;
- }
- ierr = MatMult(A, x, r);CHKERRQ(ierr);
- ierr = VecAYPX(r, -1.0, b);CHKERRQ(ierr);
- ierr = VecNorm(r, NORM_2, &resnorm);CHKERRQ(ierr);
- ierr = PetscObjectGetComm((PetscObject) A, &Acomm);CHKERRQ(ierr);
- ierr = PetscPrintf(Acomm, "Residual norm is %f.\n", (double)resnorm);CHKERRQ(ierr);
-
- if (destroyr) {
- ierr = VecDestroy(&r);CHKERRQ(ierr);
- }
- PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "PrintEnergyNorm"
-PetscErrorCode PrintEnergyNormOfDiff(Mat A, Vec x, Vec y)
-{
- PetscErrorCode ierr;
- Vec vecdiff, Avecdiff;
- PetscScalar dotprod;
- PetscReal dotabs;
- MPI_Comm Acomm;
-
- PetscFunctionBegin;
- ierr = VecDuplicate(x, &vecdiff);CHKERRQ(ierr);
- ierr = VecWAXPY(vecdiff, -1.0, x, y);CHKERRQ(ierr);
- ierr = MatGetVecs(A, NULL, &Avecdiff);CHKERRQ(ierr);
- ierr = MatMult(A, vecdiff, Avecdiff);CHKERRQ(ierr);
- ierr = VecDot(vecdiff, Avecdiff, &dotprod);CHKERRQ(ierr);
- dotabs = PetscAbsScalar(dotprod);
- ierr = PetscObjectGetComm((PetscObject) A, &Acomm);CHKERRQ(ierr);
- ierr = PetscPrintf(Acomm, "Energy norm %f.\n", (double)dotabs);CHKERRQ(ierr);
- ierr = VecDestroy(&vecdiff);CHKERRQ(ierr);
- ierr = VecDestroy(&Avecdiff);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-/* -------------------------------------------------------------------------- */
-/*
- PCDestroyLevel_ASA - Destroys one level of the ASA preconditioner
-
- Input Parameter:
-. asa_lev - pointer to level that should be destroyed
-
-*/
-#undef __FUNCT__
-#define __FUNCT__ "PCDestroyLevel_ASA"
-PetscErrorCode PCDestroyLevel_ASA(PC_ASA_level *asa_lev)
-{
- PetscErrorCode ierr;
-
- PetscFunctionBegin;
- ierr = MatDestroy(&(asa_lev->A));CHKERRQ(ierr);
- ierr = MatDestroy(&(asa_lev->B));CHKERRQ(ierr);
- ierr = VecDestroy(&(asa_lev->x));CHKERRQ(ierr);
- ierr = VecDestroy(&(asa_lev->b));CHKERRQ(ierr);
- ierr = VecDestroy(&(asa_lev->r));CHKERRQ(ierr);
-
- if (asa_lev->dm) {ierr = DMDestroy(&asa_lev->dm);CHKERRQ(ierr);}
-
- ierr = MatDestroy(&(asa_lev->agg));CHKERRQ(ierr);
- ierr = PetscFree(asa_lev->loc_agg_dofs);CHKERRQ(ierr);
- ierr = MatDestroy(&(asa_lev->agg_corr));CHKERRQ(ierr);
- ierr = MatDestroy(&(asa_lev->bridge_corr));CHKERRQ(ierr);
-
- ierr = MatDestroy(&(asa_lev->P));CHKERRQ(ierr);
- ierr = MatDestroy(&(asa_lev->Pt));CHKERRQ(ierr);
- ierr = MatDestroy(&(asa_lev->smP));CHKERRQ(ierr);
- ierr = MatDestroy(&(asa_lev->smPt));CHKERRQ(ierr);
-
- if (asa_lev->smoothd != asa_lev->smoothu) {
- if (asa_lev->smoothd) {ierr = KSPDestroy(&asa_lev->smoothd);CHKERRQ(ierr);}
- }
- if (asa_lev->smoothu) {ierr = KSPDestroy(&asa_lev->smoothu);CHKERRQ(ierr);}
-
- ierr = PetscFree(asa_lev);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-/* -------------------------------------------------------------------------- */
-/*
- PCComputeSpectralRadius_ASA - Computes the spectral radius of asa_lev->A
- and stores it it asa_lev->spec_rad
-
- Input Parameters:
-. asa_lev - the level we are treating
-
- Compute spectral radius with sqrt(||A||_1 ||A||_inf) >= ||A||_2 >= rho(A)
-
-*/
-#undef __FUNCT__
-#define __FUNCT__ "PCComputeSpectralRadius_ASA"
-PetscErrorCode PCComputeSpectralRadius_ASA(PC_ASA_level *asa_lev)
-{
- PetscErrorCode ierr;
- PetscReal norm_1, norm_inf;
-
- PetscFunctionBegin;
- ierr = MatNorm(asa_lev->A, NORM_1, &norm_1);CHKERRQ(ierr);