Actual source code: ex311.c
1: static const char help[] = "Test MatDiagonalScale() on dense matrices with scaling Vecs of any type\n\n";
3: // Contributed by: Steven Dargaville
5: #include <petscmat.h>
6: #include <petsc/private/petscimpl.h>
8: /* The scaling Vecs take their type from -vec_type, independently of the Mat type set
9: with -mat_type. Scaling is verified against a duplicate matrix scaled with VECSTANDARD
10: Vecs holding the same values: for a device matrix both scalings execute the same kernel
11: on identical data, so MatEqual() is exact, and no cross-type Vec or Mat copies are
12: needed */
13: int main(int argc, char **args)
14: {
15: Mat A, B;
16: Vec l, r, lstd, rstd;
17: PetscInt m = 5, n = 4, mloc, nloc, rstart, rend;
18: PetscBool equal = PETSC_FALSE, check_copies = PETSC_FALSE;
19: PetscLogEvent event;
20: PetscMPIInt size;
22: PetscFunctionBeginUser;
23: PetscCall(PetscInitialize(&argc, &args, NULL, help));
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
26: PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_copies", &check_copies, NULL));
27: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
28: /* PetscSF stages device buffers through the host when GPU-aware MPI is disabled. */
29: if (check_copies && size > 1 && !use_gpu_aware_mpi) {
30: PetscCall(PetscInfo(NULL, "Skipping copy checks because GPU-aware MPI is disabled; numerical checks remain enabled\n"));
31: check_copies = PETSC_FALSE;
32: }
33: if (check_copies) PetscCall(PetscLogDefaultBegin());
34: PetscCall(PetscLogEventRegister("ScaleCheck", MAT_CLASSID, &event));
36: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
37: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
38: PetscCall(MatSetType(A, MATDENSE));
39: PetscCall(MatSetFromOptions(A));
40: PetscCall(MatSetUp(A));
41: PetscCall(MatSetRandom(A, NULL));
42: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
43: PetscCall(MatGetLocalSize(A, &mloc, &nloc));
45: // l and lstd match A's row layout, r and rstd its column layout
46: PetscCall(VecCreate(PETSC_COMM_WORLD, &l));
47: PetscCall(VecSetSizes(l, mloc, m));
48: PetscCall(VecSetFromOptions(l));
49: PetscCall(VecCreate(PETSC_COMM_WORLD, &lstd));
50: PetscCall(VecSetSizes(lstd, mloc, m));
51: PetscCall(VecSetType(lstd, VECSTANDARD));
52: PetscCall(VecCreate(PETSC_COMM_WORLD, &r));
53: PetscCall(VecSetSizes(r, nloc, n));
54: PetscCall(VecSetFromOptions(r));
55: PetscCall(VecCreate(PETSC_COMM_WORLD, &rstd));
56: PetscCall(VecSetSizes(rstd, nloc, n));
57: PetscCall(VecSetType(rstd, VECSTANDARD));
59: PetscCall(VecGetOwnershipRange(l, &rstart, &rend));
60: for (PetscInt i = rstart; i < rend; i++) {
61: PetscCall(VecSetValue(l, i, (PetscScalar)(i + 2), INSERT_VALUES));
62: PetscCall(VecSetValue(lstd, i, (PetscScalar)(i + 2), INSERT_VALUES));
63: }
64: PetscCall(VecGetOwnershipRange(r, &rstart, &rend));
65: for (PetscInt j = rstart; j < rend; j++) {
66: PetscCall(VecSetValue(r, j, (PetscScalar)(j + 3), INSERT_VALUES));
67: PetscCall(VecSetValue(rstd, j, (PetscScalar)(j + 3), INSERT_VALUES));
68: }
69: PetscCall(VecAssemblyBegin(l));
70: PetscCall(VecAssemblyEnd(l));
71: PetscCall(VecAssemblyBegin(lstd));
72: PetscCall(VecAssemblyEnd(lstd));
73: PetscCall(VecAssemblyBegin(r));
74: PetscCall(VecAssemblyEnd(r));
75: PetscCall(VecAssemblyBegin(rstd));
76: PetscCall(VecAssemblyEnd(rstd));
78: // left scaling only
79: PetscCall(MatDiagonalScale(A, l, NULL));
80: PetscCall(MatDiagonalScale(B, lstd, NULL));
81: PetscCall(MatEqual(A, B, &equal));
82: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Left scaling gives the wrong result");
84: // right scaling on top of the left scaling
85: PetscCall(MatDiagonalScale(A, NULL, r));
86: PetscCall(MatDiagonalScale(B, NULL, rstd));
87: PetscCall(MatEqual(A, B, &equal));
88: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Right scaling gives the wrong result");
90: // both sides in one call; l and r were already consumed by the calls above, so any copy logged here is the matrix
91: PetscCall(PetscLogEventBegin(event, 0, 0, 0, 0));
92: PetscCall(MatDiagonalScale(A, l, r));
93: PetscCall(PetscLogEventEnd(event, 0, 0, 0, 0));
94: PetscCall(MatDiagonalScale(B, lstd, rstd));
95: PetscCall(MatEqual(A, B, &equal));
96: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Two-sided scaling gives the wrong result");
98: #if PetscDefined(HAVE_DEVICE)
99: /* a device matrix must never come back to the host to be scaled, and a device-resident scaling Vec must be
100: consumed in place; a host scaling Vec is copied to the device by design, so only the first check applies then */
101: if (check_copies) {
102: PetscEventPerfInfo info;
103: const PetscScalar *array;
104: PetscMemType mtype;
106: PetscCall(VecGetArrayReadAndMemType(l, &array, &mtype));
107: PetscCall(VecRestoreArrayReadAndMemType(l, &array));
108: PetscCall(PetscLogEventGetPerfInfo(PETSC_DETERMINE, event, &info));
109: PetscCheck(info.GpuToCpuCount == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%g unexpected GPU to CPU copies (%g bytes) in MatDiagonalScale()", info.GpuToCpuCount, info.GpuToCpuSize);
110: PetscCheck(!PetscMemTypeDevice(mtype) || info.CpuToGpuCount == 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "%g unexpected CPU to GPU copies (%g bytes) in MatDiagonalScale() with device scaling Vecs", info.CpuToGpuCount, info.CpuToGpuSize);
111: }
112: #endif
114: PetscCall(VecDestroy(&l));
115: PetscCall(VecDestroy(&lstd));
116: PetscCall(VecDestroy(&r));
117: PetscCall(VecDestroy(&rstd));
118: PetscCall(MatDestroy(&A));
119: PetscCall(MatDestroy(&B));
120: PetscCall(PetscFinalize());
121: return 0;
122: }
124: /*TEST
126: test:
127: suffix: cpu
128: nsize: {{1 2}}
129: output_file: output/empty.out
131: test:
132: suffix: kokkos
133: nsize: {{1 2}}
134: requires: kokkos_kernels
135: args: -vec_type kokkos
136: output_file: output/empty.out
138: # -check_copies requires GPU-aware MPI: without it PetscSF stages the device buffers of the
139: # right scaling through the host, and those copies are logged inside MatDiagonalScale()
141: testset:
142: nsize: {{1 2}}
143: requires: cuda
144: args: -mat_type densecuda -vec_type {{cuda standard}}
145: output_file: output/empty.out
147: test:
148: suffix: cuda
150: test:
151: suffix: cuda_copies
152: requires: defined(PETSC_HAVE_MPI_GPU_AWARE) defined(PETSC_USE_LOG)
153: args: -check_copies
155: testset:
156: nsize: {{1 2}}
157: requires: cuda kokkos_kernels
158: args: -mat_type densecuda -vec_type kokkos
159: output_file: output/empty.out
161: test:
162: suffix: densecuda_vec_kokkos
164: test:
165: suffix: densecuda_vec_kokkos_copies
166: requires: defined(PETSC_HAVE_MPI_GPU_AWARE) defined(PETSC_USE_LOG)
167: args: -check_copies
169: testset:
170: nsize: {{1 2}}
171: requires: hip
172: args: -mat_type densehip -vec_type {{hip standard}}
173: output_file: output/empty.out
175: test:
176: suffix: hip
178: test:
179: suffix: hip_copies
180: requires: defined(PETSC_HAVE_MPI_GPU_AWARE) defined(PETSC_USE_LOG)
181: args: -check_copies
183: testset:
184: nsize: {{1 2}}
185: requires: hip kokkos_kernels
186: args: -mat_type densehip -vec_type kokkos
187: output_file: output/empty.out
189: test:
190: suffix: densehip_vec_kokkos
192: test:
193: suffix: densehip_vec_kokkos_copies
194: requires: defined(PETSC_HAVE_MPI_GPU_AWARE) defined(PETSC_USE_LOG)
195: args: -check_copies
197: TEST*/