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