Actual source code: ex212.c

  1: static char help[] = "Test MatCreateSubMatrix with -mat_type nest and block sizes.\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat                    A, B, C, mats[6];
  8:   IS                     rows[2];
  9:   ISLocalToGlobalMapping cmap, rmap;
 10:   const PetscInt         indices[3] = {0, 1, 2};
 11:   PetscInt               i;
 12:   PetscMPIInt            size;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 16:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 17:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only coded for one process");
 18:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 19:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
 20:   PetscCall(MatSetBlockSizes(A, 2, 1));
 21:   PetscCall(MatSetType(A, MATAIJ));
 22:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 2, 1, indices, PETSC_COPY_VALUES, &rmap));
 23:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 1, indices, PETSC_COPY_VALUES, &cmap));
 24:   PetscCall(MatSetLocalToGlobalMapping(A, rmap, cmap));
 25:   PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
 26:   PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
 27:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 2, indices, PETSC_COPY_VALUES, &rmap));
 28:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, 3, indices, PETSC_COPY_VALUES, &cmap));
 29:   PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
 30:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, 2, 3));
 31:   PetscCall(MatSetBlockSizes(B, 1, 1));
 32:   PetscCall(MatSetType(B, MATAIJ));
 33:   PetscCall(MatSetLocalToGlobalMapping(B, rmap, cmap));
 34:   PetscCall(ISLocalToGlobalMappingDestroy(&rmap));
 35:   PetscCall(ISLocalToGlobalMappingDestroy(&cmap));
 36:   PetscCall(MatSetUp(A));
 37:   PetscCall(MatSetUp(B));
 38:   mats[0] = A;
 39:   mats[1] = B;
 40:   mats[2] = A;
 41:   mats[3] = NULL;
 42:   mats[4] = B;
 43:   mats[5] = A;
 44:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 3, NULL, mats, &C));
 45:   PetscCall(MatSetUp(C));
 46:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 47:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
 48:   PetscCall(MatView(C, NULL));
 49:   PetscCall(MatNestGetISs(C, rows, NULL));
 50:   for (i = 0; i < 2; i++) {
 51:     Mat      submat;
 52:     IS       cols[3];
 53:     PetscInt j;
 54:     PetscCall(ISView(rows[i], NULL));
 55:     PetscCall(MatCreateSubMatrix(C, rows[i], NULL, MAT_INITIAL_MATRIX, &submat));
 56:     PetscCall(MatView(submat, NULL));
 57:     PetscCall(MatNestGetISs(submat, NULL, cols));
 58:     for (j = 0; j < 3; j++) PetscCall(ISView(cols[j], NULL));
 59:     PetscCall(MatDestroy(&submat));
 60:   }
 61:   PetscCall(MatDestroy(&C));
 62:   PetscCall(MatDestroy(&B));
 63:   PetscCall(MatDestroy(&A));
 64:   PetscCall(PetscFinalize());
 65:   return 0;
 66: }

 68: /*TEST

 70:     test:

 72: TEST*/