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*/