1 static char help[] = "Tests MatSetValuesBlockedStencil() in 3d.\n\n"; 2 3 #include <petscmat.h> 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 7 int main(int argc, char **argv) 8 { 9 PetscInt M = 3, N = 4, P = 2, s = 1, w = 2, i, m = PETSC_DECIDE, n = PETSC_DECIDE, p = PETSC_DECIDE; 10 DM da; 11 Mat mat; 12 DMDAStencilType stencil_type = DMDA_STENCIL_BOX; 13 PetscBool flg = PETSC_FALSE; 14 MatStencil idx[2], idy[2]; 15 PetscScalar *values; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL)); 20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL)); 21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-P", &P, NULL)); 22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 23 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 24 PetscCall(PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL)); 25 PetscCall(PetscOptionsGetInt(NULL, NULL, "-s", &s, NULL)); 26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-w", &w, NULL)); 27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-star", &flg, NULL)); 28 if (flg) stencil_type = DMDA_STENCIL_STAR; 29 30 /* Create distributed array and get vectors */ 31 PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, stencil_type, M, N, P, m, n, p, w, s, 0, 0, 0, &da)); 32 PetscCall(DMSetFromOptions(da)); 33 PetscCall(DMSetUp(da)); 34 PetscCall(DMCreateMatrix(da, &mat)); 35 PetscCall(MatSetFromOptions(mat)); 36 37 idx[0].i = 1; 38 idx[0].j = 1; 39 idx[0].k = 0; 40 idx[1].i = 2; 41 idx[1].j = 1; 42 idx[1].k = 0; 43 idy[0].i = 1; 44 idy[0].j = 2; 45 idy[0].k = 0; 46 idy[1].i = 2; 47 idy[1].j = 2; 48 idy[1].k = 0; 49 PetscCall(PetscMalloc1(2 * 2 * w * w, &values)); 50 for (i = 0; i < 2 * 2 * w * w; i++) values[i] = i; 51 PetscCall(MatSetValuesBlockedStencil(mat, 2, idx, 2, idy, values, INSERT_VALUES)); 52 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 53 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 54 55 flg = PETSC_FALSE; 56 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_matdiagonalscalelocal", &flg, NULL)); 57 if (flg) { 58 Vec vec, scaling; 59 PetscInt size; 60 PetscMPIInt rank; 61 62 PetscCall(DMGetLocalVector(da, &vec)); 63 /* We need to duplicate vec, since we are not allowed to write to its ghost values, yet we need a vector in which 64 we can write to these locations for testing MatDiagonalScaleLocal(). */ 65 PetscCall(VecDuplicate(vec, &scaling)); 66 PetscCall(VecGetLocalSize(scaling, &size)); 67 PetscCall(PetscFree(values)); 68 PetscCall(VecGetArrayWrite(scaling, &values)); 69 PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank)); 70 for (i = 0; i < size; i++) values[i] = (PetscScalar)(rank + 1); 71 PetscCall(VecRestoreArrayWrite(scaling, &values)); 72 PetscCall(MatDiagonalScaleLocal(mat, scaling)); 73 PetscCall(DMRestoreLocalVector(da, &vec)); 74 PetscCall(DMRestoreLocalVector(da, &scaling)); 75 } 76 77 /* Free memory */ 78 PetscCall(PetscFree(values)); 79 PetscCall(MatDestroy(&mat)); 80 PetscCall(DMDestroy(&da)); 81 PetscCall(PetscFinalize()); 82 return 0; 83 } 84 85 /*TEST 86 87 test: 88 suffix: baij 89 nsize: 2 90 args: -mat_type baij 91 output_file: output/empty.out 92 93 test: 94 suffix: aij 95 nsize: 2 96 args: -mat_type aij 97 output_file: output/empty.out 98 99 test: 100 suffix: baij-diagonalscalelocal 101 nsize: 2 102 args: -mat_type baij -test_matdiagonalscalelocal 103 output_file: output/empty.out 104 105 test: 106 suffix: aij-diagonalscalelocal 107 nsize: 2 108 args: -mat_type aij -test_matdiagonalscalelocal 109 output_file: output/empty.out 110 111 TEST*/ 112