1 static char help[] = "Tests matrix ordering routines.\n\n";
2
3 #include <petscmat.h>
4 extern PetscErrorCode MatGetOrdering_myordering(Mat, MatOrderingType, IS *, IS *);
5
main(int argc,char ** args)6 int main(int argc, char **args)
7 {
8 Mat C, Cperm;
9 PetscInt i, j, m = 5, n = 5, Ii, J, ncols;
10 PetscScalar v;
11 PetscMPIInt size;
12 IS rperm, cperm, icperm;
13 const PetscInt *rperm_ptr, *cperm_ptr, *cols;
14 const PetscScalar *vals;
15 PetscBool TestMyorder = PETSC_FALSE;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &args, NULL, help));
19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
20 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
21
22 /* create the matrix for the five point stencil, YET AGAIN */
23 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m * n, m * n, 5, NULL, &C));
24 PetscCall(MatSetUp(C));
25 for (i = 0; i < m; i++) {
26 for (j = 0; j < n; j++) {
27 v = -1.0;
28 Ii = j + n * i;
29 if (i > 0) {
30 J = Ii - n;
31 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
32 }
33 if (i < m - 1) {
34 J = Ii + n;
35 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
36 }
37 if (j > 0) {
38 J = Ii - 1;
39 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
40 }
41 if (j < n - 1) {
42 J = Ii + 1;
43 PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
44 }
45 v = 4.0;
46 PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
47 }
48 }
49 PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
50 PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
51
52 PetscCall(MatGetOrdering(C, MATORDERINGND, &rperm, &cperm));
53 PetscCall(ISView(rperm, PETSC_VIEWER_STDOUT_SELF));
54 PetscCall(ISDestroy(&rperm));
55 PetscCall(ISDestroy(&cperm));
56
57 PetscCall(MatGetOrdering(C, MATORDERINGRCM, &rperm, &cperm));
58 PetscCall(ISView(rperm, PETSC_VIEWER_STDOUT_SELF));
59 PetscCall(ISDestroy(&rperm));
60 PetscCall(ISDestroy(&cperm));
61
62 PetscCall(MatGetOrdering(C, MATORDERINGQMD, &rperm, &cperm));
63 PetscCall(ISView(rperm, PETSC_VIEWER_STDOUT_SELF));
64 PetscCall(ISDestroy(&rperm));
65 PetscCall(ISDestroy(&cperm));
66
67 /* create Cperm = rperm*C*icperm */
68 PetscCall(PetscOptionsGetBool(NULL, NULL, "-testmyordering", &TestMyorder, NULL));
69 if (TestMyorder) {
70 PetscCall(MatGetOrdering_myordering(C, MATORDERINGQMD, &rperm, &cperm));
71 printf("myordering's rperm:\n");
72 PetscCall(ISView(rperm, PETSC_VIEWER_STDOUT_SELF));
73 PetscCall(ISInvertPermutation(cperm, PETSC_DECIDE, &icperm));
74 PetscCall(ISGetIndices(rperm, &rperm_ptr));
75 PetscCall(ISGetIndices(icperm, &cperm_ptr));
76 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m * n, m * n, 5, NULL, &Cperm));
77 for (i = 0; i < m * n; i++) {
78 PetscCall(MatGetRow(C, rperm_ptr[i], &ncols, &cols, &vals));
79 for (j = 0; j < ncols; j++) {
80 /* printf(" (%d %d %g)\n",i,cperm_ptr[cols[j]],vals[j]); */
81 PetscCall(MatSetValues(Cperm, 1, &i, 1, &cperm_ptr[cols[j]], &vals[j], INSERT_VALUES));
82 }
83 }
84 PetscCall(MatAssemblyBegin(Cperm, MAT_FINAL_ASSEMBLY));
85 PetscCall(MatAssemblyEnd(Cperm, MAT_FINAL_ASSEMBLY));
86 PetscCall(ISRestoreIndices(rperm, &rperm_ptr));
87 PetscCall(ISRestoreIndices(icperm, &cperm_ptr));
88
89 PetscCall(ISDestroy(&rperm));
90 PetscCall(ISDestroy(&cperm));
91 PetscCall(ISDestroy(&icperm));
92 PetscCall(MatDestroy(&Cperm));
93 }
94
95 PetscCall(MatDestroy(&C));
96 PetscCall(PetscFinalize());
97 return 0;
98 }
99
100 #include <petsc/private/matimpl.h>
101 /* This is modified from MatGetOrdering_Natural() */
MatGetOrdering_myordering(Mat mat,MatOrderingType type,IS * irow,IS * icol)102 PetscErrorCode MatGetOrdering_myordering(Mat mat, MatOrderingType type, IS *irow, IS *icol)
103 {
104 PetscInt n, i, *ii;
105 PetscBool done;
106 MPI_Comm comm;
107
108 PetscFunctionBegin;
109 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
110 PetscCall(MatGetRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, &n, NULL, NULL, &done));
111 PetscCall(MatRestoreRowIJ(mat, 0, PETSC_FALSE, PETSC_TRUE, NULL, NULL, NULL, &done));
112 PetscCheck(done, PETSC_COMM_WORLD, PETSC_ERR_SUP, "MatRestoreRowIJ fails!");
113 /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
114 PetscCall(PetscMalloc1(n, &ii));
115 for (i = 0; i < n; i++) ii[i] = n - i - 1; /* replace your index here */
116 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_COPY_VALUES, irow));
117 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, ii, PETSC_OWN_POINTER, icol));
118 PetscCall(ISSetPermutation(*irow));
119 PetscCall(ISSetPermutation(*icol));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
123 /*TEST
124
125 test:
126
127 TEST*/
128