Actual source code: ex266.c

  1: static char help[] = "Test MatDuplicate() with new nonzeros on the duplicate\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;
 10:   PetscScalar   *vals;
 11:   PetscMPIInt    rank;

 13:   // clang-format off
 14:   // i0/j0[] has a dense diagonal
 15:   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};
 16:   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};

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

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

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

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

 37:   PetscFunctionBeginUser;
 38:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 39:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 41:   mycoo = coos[rank / 3];

 43:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 44:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
 45:   PetscCall(MatSetFromOptions(A));

 47:   // Assemble matrix A with the full arrays
 48:   PetscCall(PetscMalloc1(mycoo.n, &vals));
 49:   for (k = 0; k < mycoo.n; k++) {
 50:     vals[k] = mycoo.j[k];
 51:     PetscCall(MatSetValue(A, mycoo.i[k], mycoo.j[k], vals[k], ADD_VALUES));
 52:   }
 53:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 54:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 56:   // Assemble matrix B with the 1st half of the arrays
 57:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 58:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, M, N));
 59:   PetscCall(MatSetFromOptions(B));
 60:   for (k = 0; k < mycoo.n / 2; k++) PetscCall(MatSetValue(B, mycoo.i[k], mycoo.j[k], vals[k], ADD_VALUES));
 61:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 62:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

 64:   // Duplicate B to C and continue adding nozeros to C with the 2nd half
 65:   PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &C));
 66:   for (k = mycoo.n / 2; k < mycoo.n; k++) PetscCall(MatSetValue(C, mycoo.i[k], mycoo.j[k], vals[k], ADD_VALUES));
 67:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 68:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 70:   // Test if A == C
 71:   PetscCall(MatMultEqual(A, C, 10, &equal));
 72:   if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatDuplicate() on regular matrices failed\n"));

 74:   PetscCall(PetscFree(vals));
 75:   PetscCall(MatDestroy(&A));
 76:   PetscCall(MatDestroy(&B));
 77:   PetscCall(MatDestroy(&C));
 78:   PetscCall(PetscFinalize());
 79:   return 0;
 80: }

 82: /*TEST

 84:   test:
 85:     nsize: 3
 86:     output_file: output/empty.out

 88: TEST*/