Actual source code: ex108.c

  1: static char help[] = "Testing MatCreateSeqBAIJWithArrays() and MatCreateSeqSBAIJWithArrays().\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat             A, B, As;
  8:   const PetscInt *ai, *aj;
  9:   PetscInt        i, j, k, nz, n, asi[] = {0, 2, 3, 4, 6, 7};
 10:   PetscInt        asj[] = {0, 4, 1, 2, 3, 4, 4};
 11:   PetscScalar     asa[7], *aa;
 12:   PetscRandom     rctx;
 13:   PetscMPIInt     size;
 14:   PetscBool       flg;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 18:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 19:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

 21:   /* Create a aij matrix for checking */
 22:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, 5, 5, 2, NULL, &A));
 23:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
 24:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
 25:   PetscCall(PetscRandomSetFromOptions(rctx));

 27:   k = 0;
 28:   for (i = 0; i < 5; i++) {
 29:     nz = asi[i + 1] - asi[i]; /* length of i_th row of A */
 30:     for (j = 0; j < nz; j++) {
 31:       PetscCall(PetscRandomGetValue(rctx, &asa[k]));
 32:       PetscCall(MatSetValues(A, 1, &i, 1, &asj[k], &asa[k], INSERT_VALUES));
 33:       if (i != asj[k]) { /* insert symmetric entry */
 34:         PetscCall(MatSetValues(A, 1, &asj[k], 1, &i, &asa[k], INSERT_VALUES));
 35:       }
 36:       k++;
 37:     }
 38:   }
 39:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 40:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

 42:   /* Create a baij matrix using MatCreateSeqBAIJWithArrays() */
 43:   PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ai, &aj, &flg));
 44:   PetscCall(MatSeqAIJGetArray(A, &aa));
 45:   /* WARNING: This sharing is dangerous if either A or B is later assembled */
 46:   PetscCall(MatCreateSeqBAIJWithArrays(PETSC_COMM_SELF, 1, 5, 5, (PetscInt *)ai, (PetscInt *)aj, aa, &B));
 47:   PetscCall(MatSeqAIJRestoreArray(A, &aa));
 48:   PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ai, &aj, &flg));
 49:   PetscCall(MatMultEqual(A, B, 10, &flg));
 50:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "MatMult(A,B) are NOT equal");

 52:   /* Create a sbaij matrix using MatCreateSeqSBAIJWithArrays() */
 53:   PetscCall(MatCreateSeqSBAIJWithArrays(PETSC_COMM_SELF, 1, 5, 5, asi, asj, asa, &As));
 54:   PetscCall(MatMultEqual(A, As, 10, &flg));
 55:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "MatMult(A,As) are NOT equal");

 57:   /* Free spaces */
 58:   PetscCall(PetscRandomDestroy(&rctx));
 59:   PetscCall(MatDestroy(&A));
 60:   PetscCall(MatDestroy(&B));
 61:   PetscCall(MatDestroy(&As));
 62:   PetscCall(PetscFinalize());
 63:   return 0;
 64: }

 66: /*TEST

 68:   test:
 69:     output_file: output/empty.out

 71: TEST*/