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 */
9
10 #include <petscmat.h>
11
main(int argc,char ** argv)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;
22
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));
27
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));
36
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));
44
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));
59
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 }
69
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));
77
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 }
84
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 }
94
95 /*TEST
96
97 build:
98 requires: elemental
99
100 test:
101 nsize: 6
102 output_file: output/empty.out
103
104 test:
105 suffix: 2
106 nsize: 6
107 args: -mat_type aij
108 output_file: output/empty.out
109
110 test:
111 suffix: 3
112 nsize: 6
113 args: -mat_type elemental
114 output_file: output/empty.out
115
116 TEST*/
117