Actual source code: ex103.c

  1: static char help[] = "Test MatSetValues() by converting MATDENSE to MATELEMENTAL. \n\
  2: Modified from the code contributed by Yaning Liu @lbl.gov \n\n";
  3: /*
  4:  Example:
  5:    mpiexec -n <np> ./ex103
  6:    mpiexec -n <np> ./ex103 -mat_type elemental -mat_view
  7:    mpiexec -n <np> ./ex103 -mat_type aij
  8: */

 10: #include <petscmat.h>

 12: int main(int argc, char **argv)
 13: {
 14:   Mat             A, A_elemental;
 15:   PetscInt        i, j, M = 10, N = 5, nrows, ncols;
 16:   PetscMPIInt     rank, size;
 17:   IS              isrows, iscols;
 18:   const PetscInt *rows, *cols;
 19:   PetscScalar    *v;
 20:   MatType         type;
 21:   PetscBool       isDense, isAIJ, flg;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 25:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 26:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 28:   /* Creat a matrix */
 29:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 30:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
 31:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 32:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
 33:   PetscCall(MatSetType(A, MATDENSE));
 34:   PetscCall(MatSetFromOptions(A));
 35:   PetscCall(MatSetUp(A));

 37:   /* Set local matrix entries */
 38:   PetscCall(MatGetOwnershipIS(A, &isrows, &iscols));
 39:   PetscCall(ISGetLocalSize(isrows, &nrows));
 40:   PetscCall(ISGetIndices(isrows, &rows));
 41:   PetscCall(ISGetLocalSize(iscols, &ncols));
 42:   PetscCall(ISGetIndices(iscols, &cols));
 43:   PetscCall(PetscMalloc1(nrows * ncols, &v));

 45:   for (i = 0; i < nrows; i++) {
 46:     for (j = 0; j < ncols; j++) {
 47:       if (size == 1) {
 48:         v[i * ncols + j] = (PetscScalar)(i + j);
 49:       } else {
 50:         v[i * ncols + j] = (PetscScalar)rank + j * 0.1;
 51:       }
 52:     }
 53:   }
 54:   PetscCall(MatSetValues(A, nrows, rows, ncols, cols, v, INSERT_VALUES));
 55:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 56:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 57:   //PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%" PetscInt_FMT "] local nrows %" PetscInt_FMT ", ncols %" PetscInt_FMT "\n",rank,nrows,ncols));
 58:   //PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT));

 60:   /* Test MatSetValues() by converting A to A_elemental */
 61:   PetscCall(MatGetType(A, &type));
 62:   if (size == 1) {
 63:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &isDense));
 64:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &isAIJ));
 65:   } else {
 66:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIDENSE, &isDense));
 67:     PetscCall(PetscObjectTypeCompare((PetscObject)A, MATMPIAIJ, &isAIJ));
 68:   }

 70:   if (isDense || isAIJ) {
 71:     Mat Aexplicit;
 72:     PetscCall(MatConvert(A, MATELEMENTAL, MAT_INITIAL_MATRIX, &A_elemental));
 73:     PetscCall(MatComputeOperator(A_elemental, isAIJ ? MATAIJ : MATDENSE, &Aexplicit));
 74:     PetscCall(MatMultEqual(Aexplicit, A_elemental, 5, &flg));
 75:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Aexplicit != A_elemental.");
 76:     PetscCall(MatDestroy(&Aexplicit));

 78:     /* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */
 79:     PetscCall(MatConvert(A, MATELEMENTAL, MAT_INPLACE_MATRIX, &A));
 80:     PetscCall(MatMultEqual(A_elemental, A, 5, &flg));
 81:     PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A_elemental != A.");
 82:     PetscCall(MatDestroy(&A_elemental));
 83:   }

 85:   PetscCall(ISRestoreIndices(isrows, &rows));
 86:   PetscCall(ISRestoreIndices(iscols, &cols));
 87:   PetscCall(ISDestroy(&isrows));
 88:   PetscCall(ISDestroy(&iscols));
 89:   PetscCall(PetscFree(v));
 90:   PetscCall(MatDestroy(&A));
 91:   PetscCall(PetscFinalize());
 92:   return 0;
 93: }

 95: /*TEST

 97:    build:
 98:       requires: elemental

100:    test:
101:       nsize: 6

103:    test:
104:       suffix: 2
105:       nsize: 6
106:       args: -mat_type aij
107:       output_file: output/ex103_1.out

109:    test:
110:       suffix: 3
111:       nsize: 6
112:       args: -mat_type elemental
113:       output_file: output/ex103_1.out

115: TEST*/