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:   PetscInt       k;
  8:   const PetscInt M = 18, N = 18;
  9:   PetscBool      equal, isHypre;
 10:   PetscScalar   *vals;
 11:   PetscBool      flg = PETSC_FALSE, freecoo = PETSC_FALSE, missing_diagonal = PETSC_FALSE;
 12:   PetscInt       ncoos = 1;

 14:   // clang-format off
 15:   /* Construct 18 x 18 matrices, which are big enough to have complex communication patterns but still small enough for debugging */
 16:   // i0/j0[] has a dense diagonal
 17:   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};
 18:   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};

 20:   // i0/j0[] miss some diagonals
 21:   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};
 22:   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};

 24:   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};
 25:   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};
 26:   // clang-format on

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

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

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

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

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

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

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

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

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

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

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

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

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

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

123: /*TEST

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

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

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

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

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

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

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

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

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

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

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

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

185: TEST*/