Actual source code: ex58.c

  1: static char help[] = "Tests MatTranspose() and MatEqual() for MPIAIJ matrices.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat         A, B;
  8:   PetscInt    m = 7, n, i, rstart, rend, cols[3];
  9:   PetscScalar v[3];
 10:   PetscBool   equal;
 11:   const char *eq[2];

 13:   PetscFunctionBeginUser;
 14:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 15:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 16:   n = m;

 18:   /* ------- Assemble matrix, --------- */

 20:   PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, n, 0, 0, 0, 0, &A));
 21:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
 22:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 23:   if (!rstart) {
 24:     cols[0] = 0;
 25:     cols[1] = 1;
 26:     v[0]    = 2.0;
 27:     v[1]    = -1.0;
 28:     PetscCall(MatSetValues(A, 1, &rstart, 2, cols, v, INSERT_VALUES));
 29:     rstart++;
 30:   }
 31:   if (rend == m) {
 32:     rend--;
 33:     cols[0] = rend - 1;
 34:     cols[1] = rend;
 35:     v[0]    = -1.0;
 36:     v[1]    = 2.0;
 37:     PetscCall(MatSetValues(A, 1, &rend, 2, cols, v, INSERT_VALUES));
 38:   }
 39:   v[0] = -1.0;
 40:   v[1] = 2.0;
 41:   v[2] = -1.0;
 42:   for (i = rstart; i < rend; i++) {
 43:     cols[0] = i - 1;
 44:     cols[1] = i;
 45:     cols[2] = i + 1;
 46:     PetscCall(MatSetValues(A, 1, &i, 3, cols, v, INSERT_VALUES));
 47:   }
 48:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 49:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 51:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));

 53:   PetscCall(MatEqual(A, B, &equal));

 55:   eq[0] = "not equal";
 56:   eq[1] = "equal";
 57:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Matrices are %s\n", eq[equal]));

 59:   PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &B));
 60:   PetscCall(MatEqual(A, B, &equal));
 61:   if (!equal) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatTranspose with MAT_REUSE_MATRIX failed"));

 63:   /* Free data structures */
 64:   PetscCall(MatDestroy(&A));
 65:   PetscCall(MatDestroy(&B));

 67:   PetscCall(PetscFinalize());
 68:   return 0;
 69: }

 71: /*TEST

 73:     test:

 75:     test:
 76:       suffix: 2
 77:       nsize: 2
 78:       output_file: output/ex58_1.out

 80: TEST*/