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