Actual source code: ex254.c

  1: static char help[] = "Test MatSetValuesCOO() for MPIAIJ and its subclasses \n\n";

  3: #include <petscmat.h>
  4: int main(int argc, char **args)
  5: {
  6:   Mat            A, B, C;
  7:   const PetscInt M = 18, N = 18;
  8:   PetscBool      equal, isHypre;
  9:   PetscScalar   *vals;
 10:   PetscBool      flg = PETSC_FALSE, freecoo = PETSC_FALSE, missing_diagonal = PETSC_FALSE;
 11:   PetscInt       ncoos = 1;

 13:   // clang-format off
 14:   /* Construct 18 x 18 matrices, which are big enough to have complex communication patterns but still small enough for debugging */
 15:   // i0/j0[] has a dense diagonal
 16:   PetscInt i0[] = {7, 7, 8, 8, 9, 16, 17,  9, 10, 1, 1, -2, 2, 3, 3, 14, 4, 5, 10, 13,  9,  9, 10, 1, 0, 0, 5,  5,  6, 6, 13, 13, 14, -14, 4, 4, 5, 11, 11, 12, 15, 15, 16};
 17:   PetscInt j0[] = {7, 6, 8, 4, 9, 16, 17, 16, 10, 2, 1,  3, 2, 4, 3, 14, 4, 5, 15, 13, 10, 16, 11, 2, 0, 1, 5, -11, 0, 6, 15, 17, 11,  13, 4, 8, 2, 11, 17, 12,  3, 15,  9};

 19:   // i0/j0[] miss some diagonals
 20:   PetscInt i1[] = {8, 5, 15, 16, 6, 13, 4, 17, 8,  9, 9,  10, -6, 12, 7, 3, -4, 1, 1, 2, 5,  5, 6, 14, 17, 8,  9,  9, 10, 4,  5, 10, 11, 1, 2};
 21:   PetscInt j1[] = {2, 3, 16,  9, 5, 17, 1, 13, 4, 10, 16, 11, -5, 12, 1, 7, -1, 2, 7, 3, 6, 11, 0, 11, 13, 4, 10, 16, 11, 8, -2, 15, 12, 7, 3};

 23:   PetscInt i2[] = {3, 4, 1, 10, 0, 1, 1, 2, 1, 1, 2, 2, 3, 3, 4, 4, 1, 2, 5,  5, 6, 4, 17, 0, 1, 1, 8, 5,  5, 6, 4, 7, 8, 5};
 24:   PetscInt j2[] = {7, 1, 2, 11, 5, 2, 7, 3, 2, 7, 3, 8, 4, 9, 3, 5, 7, 3, 6, 11, 0, 1, 13, 5, 2, 7, 4, 6, 11, 0, 1, 3, 4, 2};
 25:   // clang-format on

 27:   typedef struct {
 28:     PetscInt *i, *j, n;
 29:   } coo_data;

 31:   coo_data coos[3] = {
 32:     {i0, j0, PETSC_STATIC_ARRAY_LENGTH(i0)},
 33:     {i1, j1, PETSC_STATIC_ARRAY_LENGTH(i1)},
 34:     {i2, j2, PETSC_STATIC_ARRAY_LENGTH(i2)}
 35:   };
 36:   coo_data mycoo;

 38:   PetscFunctionBeginUser;
 39:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 40:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-ignore_remote", &flg, NULL));
 41:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ncoos", &ncoos, NULL));
 42:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-missing_diagonal", &missing_diagonal, NULL));

 44:   mycoo.n = 0;
 45:   if (ncoos > 1) {
 46:     PetscLayout map;

 48:     freecoo = PETSC_TRUE;
 49:     PetscCall(PetscLayoutCreate(PETSC_COMM_WORLD, &map));
 50:     PetscCall(PetscLayoutSetSize(map, ncoos));
 51:     PetscCall(PetscLayoutSetUp(map));
 52:     PetscCall(PetscLayoutGetLocalSize(map, &ncoos));
 53:     for (PetscInt i = 0; i < ncoos; i++) mycoo.n += coos[i % 3].n;
 54:     PetscCall(PetscMalloc2(mycoo.n, &mycoo.i, mycoo.n, &mycoo.j));
 55:     mycoo.n = 0;
 56:     for (PetscInt i = 0; i < ncoos; i++) {
 57:       PetscCall(PetscArraycpy(mycoo.i + mycoo.n, coos[i % 3].i, coos[i % 3].n));
 58:       PetscCall(PetscArraycpy(mycoo.j + mycoo.n, coos[i % 3].j, coos[i % 3].n));
 59:       mycoo.n += coos[i % 3].n;
 60:     }
 61:     PetscCall(PetscLayoutDestroy(&map));
 62:   } else if (ncoos == 1 && PetscGlobalRank < 3) mycoo = coos[PetscGlobalRank];

 64:   if (missing_diagonal && PetscGlobalRank == 0) mycoo = coos[1];

 66:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 67:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
 68:   PetscCall(MatSetType(A, MATAIJ));
 69:   // Do not preallocate A to also test MatHash with MAT_IGNORE_OFF_PROC_ENTRIES
 70:   // PetscCall(MatSeqAIJSetPreallocation(A, 2, NULL));
 71:   // PetscCall(MatMPIAIJSetPreallocation(A, 2, NULL, 2, NULL));
 72:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 73:   PetscCall(MatSetOption(A, MAT_IGNORE_OFF_PROC_ENTRIES, flg));

 75:   PetscCall(PetscMalloc1(mycoo.n, &vals));
 76:   for (PetscInt k = 0; k < mycoo.n; k++) {
 77:     vals[k] = mycoo.j[k];
 78:     PetscCall(MatSetValue(A, mycoo.i[k], mycoo.j[k], vals[k], ADD_VALUES));
 79:   }
 80:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 81:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 82:   PetscCall(MatViewFromOptions(A, NULL, "-a_view"));

 84:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 85:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, N));
 86:   PetscCall(MatSetFromOptions(B));
 87:   PetscCall(MatSetOption(B, MAT_IGNORE_OFF_PROC_ENTRIES, flg));
 88:   PetscCall(MatSetPreallocationCOO(B, mycoo.n, mycoo.i, mycoo.j));

 90:   /* Test with ADD_VALUES on a zeroed matrix */
 91:   PetscCall(MatSetValuesCOO(B, vals, ADD_VALUES));
 92:   PetscCall(MatMultEqual(A, B, 10, &equal));
 93:   if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSetValuesCOO() failed\n"));
 94:   PetscCall(MatViewFromOptions(B, NULL, "-b_view"));

 96:   /* Test with MatDuplicate on a zeroed matrix */
 97:   PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &C));
 98:   PetscCall(MatSetValuesCOO(C, vals, ADD_VALUES));
 99:   PetscCall(MatMultEqual(A, C, 10, &equal));
100:   if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatSetValuesCOO() on duplicated matrix failed\n"));
101:   PetscCall(MatViewFromOptions(C, NULL, "-c_view"));

103:   /* Test aij->diag on COO matrix are correctly set up */
104:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &isHypre));
105:   if (!isHypre) { // TODO: MATHYPRE currently does not support MatSetValues
106:     PetscCall(MatShift(A, 2.0));
107:     PetscCall(MatShift(B, 2.0));
108:     PetscCall(MatMultEqual(A, B, 10, &equal));
109:     if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatShift() on a duplicated COO matrix failed\n"));
110:   }

112:   PetscCall(PetscFree(vals));
113:   if (freecoo) PetscCall(PetscFree2(mycoo.i, mycoo.j));
114:   PetscCall(MatDestroy(&A));
115:   PetscCall(MatDestroy(&B));
116:   PetscCall(MatDestroy(&C));

118:   PetscCall(PetscFinalize());
119:   return 0;
120: }

122: /*TEST

124:   testset:
125:     output_file: output/empty.out
126:     nsize: {{1 2 3}}
127:     args: -ignore_remote {{0 1}} -missing_diagonal {{0 1}}
128:     filter: grep -v type | grep -v "Mat Object"

130:     test:
131:       suffix: kokkos
132:       requires: kokkos_kernels
133:       args: -mat_type aijkokkos

135:     test:
136:       suffix: cuda
137:       requires: cuda
138:       args: -mat_type aijcusparse

140:     test:
141:       suffix: hip
142:       requires: hip
143:       args: -mat_type aijhipsparse

145:     test:
146:       suffix: aij
147:       args: -mat_type aij

149:     test:
150:       suffix: hypre
151:       requires: hypre
152:       args: -mat_type hypre

154:   testset:
155:     output_file: output/empty.out
156:     nsize: 1
157:     args: -ncoos 3
158:     filter: grep -v type | grep -v "Mat Object"

160:     test:
161:       suffix: 2_kokkos
162:       requires: kokkos_kernels
163:       args: -mat_type aijkokkos

165:     test:
166:       suffix: 2_cuda
167:       requires: cuda
168:       args: -mat_type aijcusparse

170:     test:
171:       suffix: 2_hip
172:       requires: hip
173:       args: -mat_type aijhipsparse

175:     test:
176:       suffix: 2_aij
177:       args: -mat_type aij

179:     test:
180:       suffix: 2_hypre
181:       requires: hypre
182:       args: -mat_type hypre

184: TEST*/