Actual source code: ex313.c
1: static const char
2: help[] = "Test MatNorm() on dense matrices, including submatrix views where the leading dimension exceeds the local row count.\nThe matrix under test is built with -mat_type, or through MatCreateDenseFromVecType() with -from_vec_type.\n\n";
4: // Contributed by: Steven Dargaville
6: #include <petscmat.h>
8: /*
9: Entries are small signed integers, so every norm is an exactly representable sum. The magnitude
10: 1 + (i % 5) + 2 * (j % 3) lies between 1 and 9 and varies with both the row and the column, so
11: the largest column sum and the largest row sum are attained at nontrivial positions. The sign
12: alternates in a checkerboard, so the absolute values taken by NORM_1 and NORM_INFINITY matter,
13: and complex builds add an imaginary part, so the modulus differs from |re| + |im| and
14: absolute-value shortcuts on the device are caught
15: */
16: static PetscScalar Entry(PetscInt i, PetscInt j)
17: {
18: PetscReal magnitude = (PetscReal)(1 + (i % 5) + 2 * (j % 3));
19: PetscReal sign = ((i + j) % 2) ? -1.0 : 1.0;
20: PetscScalar v = sign * magnitude;
22: #if PetscDefined(USE_COMPLEX)
23: v += PETSC_i * (PetscScalar)(i - 2 * j);
24: #endif
25: return v;
26: }
28: /*
29: Compare the norms of the matrix under test against those of a reference matrix holding the same entries
30: */
31: static PetscErrorCode CheckNorms(Mat A, Mat B)
32: {
33: const NormType types[] = {NORM_1, NORM_FROBENIUS, NORM_INFINITY};
34: PetscReal na, nb;
36: PetscFunctionBeginUser;
37: for (PetscInt t = 0; t < 3; t++) {
38: PetscCall(MatNorm(A, types[t], &na));
39: PetscCall(MatNorm(B, types[t], &nb));
40: PetscCheck(PetscAbsReal(na - nb) <= PETSC_SMALL + PETSC_SMALL * PetscAbsReal(nb), PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "NORM_%s mismatch: %g (tested) != %g (reference)", NormTypes[types[t]], (double)na, (double)nb);
41: }
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: int main(int argc, char **argv)
46: {
47: Mat A, B, C, Asub, Bsub;
48: PetscInt m = 12, n = 7, rstart, rend;
49: char vtype[64];
50: PetscBool from_vec_type = PETSC_FALSE;
52: PetscFunctionBeginUser;
53: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
54: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
55: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
56: PetscCall(PetscOptionsGetString(NULL, NULL, "-from_vec_type", vtype, sizeof(vtype), &from_vec_type));
57: PetscCheck(m > 2 && n > 2, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Need m > 2 and n > 2, got %" PetscInt_FMT " and %" PetscInt_FMT, m, n);
59: /* the matrix under test; MatCreateDenseFromVecType() is the construction route Kokkos codes use,
60: and it yields MATDENSECUDA or MATDENSEHIP on a device Kokkos build and MATDENSE on a host
61: Kokkos build */
62: if (from_vec_type) PetscCall(MatCreateDenseFromVecType(PETSC_COMM_WORLD, vtype, PETSC_DECIDE, PETSC_DECIDE, m, n, PETSC_DECIDE, NULL, &A));
63: else {
64: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
65: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
66: PetscCall(MatSetType(A, MATDENSE));
67: PetscCall(MatSetFromOptions(A));
68: PetscCall(MatSetUp(A));
69: }
71: /* the reference matrix, always on the host */
72: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
73: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
74: PetscCall(MatSetType(B, MATDENSE));
75: PetscCall(MatSetUp(B));
77: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
78: for (PetscInt i = rstart; i < rend; i++) {
79: for (PetscInt j = 0; j < n; j++) {
80: PetscCall(MatSetValue(A, i, j, Entry(i, j), INSERT_VALUES));
81: PetscCall(MatSetValue(B, i, j, Entry(i, j), INSERT_VALUES));
82: }
83: }
84: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
87: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
89: PetscCall(CheckNorms(A, B));
91: /* the entries of the interior block, in fresh contiguous storage where the leading dimension
92: equals the number of local rows, so it is unaffected by any leading dimension handling */
93: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
94: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m - 2, n - 2));
95: PetscCall(MatSetType(C, MATDENSE));
96: PetscCall(MatSetUp(C));
97: PetscCall(MatGetOwnershipRange(C, &rstart, &rend));
98: for (PetscInt i = rstart; i < rend; i++) {
99: for (PetscInt j = 0; j < n - 2; j++) PetscCall(MatSetValue(C, i, j, Entry(i + 1, j + 1), INSERT_VALUES));
100: }
101: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
102: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
104: /* an interior view keeps the leading dimension of its parent, so the leading dimension exceeds the number of local rows */
105: PetscCall(MatDenseGetSubMatrix(A, 1, m - 1, 1, n - 1, &Asub));
106: PetscCall(MatDenseGetSubMatrix(B, 1, m - 1, 1, n - 1, &Bsub));
107: PetscCall(CheckNorms(Asub, C));
108: PetscCall(CheckNorms(Bsub, C));
109: PetscCall(CheckNorms(Asub, Bsub));
110: PetscCall(MatDenseRestoreSubMatrix(A, &Asub));
111: PetscCall(MatDenseRestoreSubMatrix(B, &Bsub));
113: PetscCall(MatDestroy(&A));
114: PetscCall(MatDestroy(&B));
115: PetscCall(MatDestroy(&C));
116: PetscCall(PetscFinalize());
117: return 0;
118: }
120: /*TEST
122: testset:
123: nsize: {{1 2}}
124: output_file: output/empty.out
126: test:
127: suffix: cpu
129: test:
130: suffix: cuda
131: requires: cuda
132: args: -mat_type densecuda
134: test:
135: suffix: hip
136: requires: hip
137: args: -mat_type densehip
139: test:
140: suffix: from_vec_kokkos
141: requires: kokkos_kernels !sycl
142: args: -from_vec_type kokkos
144: TEST*/