Actual source code: ex13f.F90
1: !
2: ! Tests PCASMWeightedGetScaling() from Fortran with an output pointer that starts disassociated,
3: ! and PCASMWeightedSetComputeScaling() with a Fortran callback and context
4: !
5: ! -----------------------------------------------------------------------
6: #include <petsc/finclude/petscksp.h>
8: ! Fills the weights of the only local subdomain with the value passed as context
9: subroutine FillScaling(pc, local, scaling, value, ierr)
10: use petscksp
11: implicit none
13: PC pc
14: PetscInt local
15: Vec scaling
16: PetscScalar value
17: PetscErrorCode ierr
19: PetscCheck(local == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Expected a single local subdomain')
20: PetscCall(VecSet(scaling, value, ierr))
21: end subroutine
23: program main
24: use petscksp
25: implicit none
27: PC pc
28: Mat A
29: Vec x, y, supplied(1)
30: Vec, pointer :: scaling(:) => null()
31: PetscInt n, m, i, istart, iend
32: PetscInt, parameter :: nlocal = 4
33: PetscReal norm
34: PetscScalar total, three
35: PetscScalar, parameter :: one = 1.0, two = 2.0
36: PetscErrorCode ierr
37: external FillScaling
39: PetscCallA(PetscInitialize(ierr))
40: PetscCallA(PCCreate(PETSC_COMM_WORLD, pc, ierr))
41: PetscCallA(PCSetType(pc, PCASM, ierr))
42: PetscCallA(PCASMSetType(pc, PC_ASM_WEIGHTED, ierr))
44: ! No weights yet: the count is zero and the pointer stays disassociated
45: PetscCallA(PCASMWeightedGetScaling(pc, n, scaling, ierr))
46: PetscCheckA(n == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Expected no scaling vectors before PCASMWeightedSetScaling()')
47: PetscCheckA(.not. associated(scaling), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Expected a disassociated pointer before PCASMWeightedSetScaling()')
49: ! Identity operator with one subdomain per process and no overlap
50: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, nlocal, nlocal, PETSC_DETERMINE, PETSC_DETERMINE, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 0_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, A, ierr))
51: PetscCallA(MatGetOwnershipRange(A, istart, iend, ierr))
52: do i = istart, iend - 1
53: PetscCallA(MatSetValue(A, i, i, one, INSERT_VALUES, ierr))
54: end do
55: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
56: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
57: PetscCallA(PCSetOperators(pc, A, A, ierr))
58: PetscCallA(PCASMSetOverlap(pc, 0_PETSC_INT_KIND, ierr))
59: PetscCallA(PCSetUp(pc, ierr))
61: ! Supply one scaling vector; pc keeps its own reference
62: PetscCallA(MatGetLocalSize(A, m, PETSC_NULL_INTEGER, ierr))
63: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, m, supplied(1), ierr))
64: PetscCallA(VecSet(supplied(1), two, ierr))
65: PetscCallA(PCASMWeightedSetScaling(pc, 1_PETSC_INT_KIND, supplied, ierr))
66: PetscCallA(VecDestroy(supplied(1), ierr))
68: ! The getter must associate a pointer that has never been associated
69: PetscCallA(PCASMWeightedGetScaling(pc, n, scaling, ierr))
70: PetscCheckA(n == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Expected one scaling vector')
71: PetscCheckA(associated(scaling), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'PCASMWeightedGetScaling() left the output pointer disassociated')
72: PetscCheckA(size(scaling) == 1 .and. lbound(scaling, 1) == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Wrong scaling array bounds')
73: PetscCallA(VecSum(scaling(1), total, ierr))
74: PetscCheckA(abs(total - two*m) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Wrong scaling vector contents')
76: ! The count may be omitted
77: PetscCallA(PCASMWeightedGetScaling(pc, PETSC_NULL_INTEGER, scaling, ierr))
78: PetscCheckA(associated(scaling) .and. size(scaling) == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'PCASMWeightedGetScaling() failed with PETSC_NULL_INTEGER')
80: ! With an identity operator and no overlap the preconditioner is the scaling itself
81: PetscCallA(MatCreateVecs(A, x, y, ierr))
82: PetscCallA(VecSet(x, one, ierr))
83: PetscCallA(PCApply(pc, x, y, ierr))
84: PetscCallA(VecAXPY(y, -two, x, ierr))
85: PetscCallA(VecNorm(y, NORM_INFINITY, norm, ierr))
86: PetscCheckA(norm < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'PCApply() did not use the supplied weights')
88: ! PCReset() releases the weights and the getter clears the associated pointer
89: PetscCallA(PCReset(pc, ierr))
90: PetscCallA(PCASMWeightedGetScaling(pc, n, scaling, ierr))
91: PetscCheckA(n == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Expected no scaling vectors after PCReset()')
92: PetscCheckA(.not. associated(scaling), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Expected a disassociated pointer after PCReset()')
94: ! A callback registered before setup computes the weights, with its context passed through
95: three = 3.0
96: PetscCallA(PCASMWeightedSetComputeScaling(pc, FillScaling, three, ierr))
97: PetscCallA(PCSetOperators(pc, A, A, ierr))
98: PetscCallA(PCSetUp(pc, ierr))
99: PetscCallA(PCASMWeightedGetScaling(pc, n, scaling, ierr))
100: PetscCheckA(n == 1 .and. associated(scaling), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'PCASMWeightedSetComputeScaling() did not create the weights')
101: PetscCallA(PCApply(pc, x, y, ierr))
102: PetscCallA(VecAXPY(y, -three, x, ierr))
103: PetscCallA(VecNorm(y, NORM_INFINITY, norm, ierr))
104: PetscCheckA(norm < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'PCApply() did not use the computed weights')
106: ! PETSC_NULL_FUNCTION disables the callback and keeps the computed weights
107: PetscCallA(PCASMWeightedSetComputeScaling(pc, PETSC_NULL_FUNCTION, 0, ierr))
108: PetscCallA(PCASMWeightedGetScaling(pc, n, scaling, ierr))
109: PetscCheckA(n == 1 .and. associated(scaling), PETSC_COMM_SELF, PETSC_ERR_PLIB, 'Disabling the callback discarded the weights')
111: PetscCallA(VecDestroy(x, ierr))
112: PetscCallA(VecDestroy(y, ierr))
113: PetscCallA(MatDestroy(A, ierr))
114: PetscCallA(PCDestroy(pc, ierr))
115: PetscCallA(PetscFinalize(ierr))
116: end
118: !/*TEST
119: !
120: ! test:
121: ! nsize: {{1 2}}
122: ! output_file: output/empty.out
123: !TEST*/