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:   PetscErrorCode    ierr;
 12:   PetscScalar       v;
 13:   PetscMPIInt       size;
 14:   IS                rperm,cperm,icperm;
 15:   const PetscInt    *rperm_ptr,*cperm_ptr,*cols;
 16:   const PetscScalar *vals;
 17:   PetscBool         TestMyorder=PETSC_FALSE;

 19:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 20:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 21:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

 23:   /* create the matrix for the five point stencil, YET AGAIN */
 24:   MatCreateSeqAIJ(PETSC_COMM_SELF,m*n,m*n,5,NULL,&C);
 25:   MatSetUp(C);
 26:   for (i=0; i<m; i++) {
 27:     for (j=0; j<n; j++) {
 28:       v = -1.0;  Ii = j + n*i;
 29:       if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 30:       if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 31:       if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 32:       if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
 33:       v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 34:     }
 35:   }
 36:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 37:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 39:   MatGetOrdering(C,MATORDERINGND,&rperm,&cperm);
 40:   ISView(rperm,PETSC_VIEWER_STDOUT_SELF);
 41:   ISDestroy(&rperm);
 42:   ISDestroy(&cperm);

 44:   MatGetOrdering(C,MATORDERINGRCM,&rperm,&cperm);
 45:   ISView(rperm,PETSC_VIEWER_STDOUT_SELF);
 46:   ISDestroy(&rperm);
 47:   ISDestroy(&cperm);

 49:   MatGetOrdering(C,MATORDERINGQMD,&rperm,&cperm);
 50:   ISView(rperm,PETSC_VIEWER_STDOUT_SELF);
 51:   ISDestroy(&rperm);
 52:   ISDestroy(&cperm);

 54:   /* create Cperm = rperm*C*icperm */
 55:   PetscOptionsGetBool(NULL,NULL,"-testmyordering",&TestMyorder,NULL);
 56:   if (TestMyorder) {
 57:     MatGetOrdering_myordering(C,MATORDERINGQMD,&rperm,&cperm);
 58:     printf("myordering's rperm:\n");
 59:     ISView(rperm,PETSC_VIEWER_STDOUT_SELF);
 60:     ISInvertPermutation(cperm,PETSC_DECIDE,&icperm);
 61:     ISGetIndices(rperm,&rperm_ptr);
 62:     ISGetIndices(icperm,&cperm_ptr);
 63:     MatCreateSeqAIJ(PETSC_COMM_SELF,m*n,m*n,5,NULL,&Cperm);
 64:     for (i=0; i<m*n; i++) {
 65:       MatGetRow(C,rperm_ptr[i],&ncols,&cols,&vals);
 66:       for (j=0; j<ncols; j++) {
 67:         /* printf(" (%d %d %g)\n",i,cperm_ptr[cols[j]],vals[j]); */
 68:         MatSetValues(Cperm,1,&i,1,&cperm_ptr[cols[j]],&vals[j],INSERT_VALUES);
 69:       }
 70:     }
 71:     MatAssemblyBegin(Cperm,MAT_FINAL_ASSEMBLY);
 72:     MatAssemblyEnd(Cperm,MAT_FINAL_ASSEMBLY);
 73:     ISRestoreIndices(rperm,&rperm_ptr);
 74:     ISRestoreIndices(icperm,&cperm_ptr);

 76:     ISDestroy(&rperm);
 77:     ISDestroy(&cperm);
 78:     ISDestroy(&icperm);
 79:     MatDestroy(&Cperm);
 80:   }

 82:   MatDestroy(&C);
 83:   PetscFinalize();
 84:   return ierr;
 85: }

 87: #include <petsc/private/matimpl.h>
 88: /* This is modified from MatGetOrdering_Natural() */
 89: PetscErrorCode MatGetOrdering_myordering(Mat mat,MatOrderingType type,IS *irow,IS *icol)
 90: {
 92:   PetscInt       n,i,*ii;
 93:   PetscBool      done;
 94:   MPI_Comm       comm;

 97:   PetscObjectGetComm((PetscObject)mat,&comm);
 98:   MatGetRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,&n,NULL,NULL,&done);
 99:   MatRestoreRowIJ(mat,0,PETSC_FALSE,PETSC_TRUE,NULL,NULL,NULL,&done);
100:   if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
101:     PetscMalloc1(n,&ii);
102:     for (i=0; i<n; i++) ii[i] = n-i-1; /* replace your index here */
103:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_COPY_VALUES,irow);
104:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,PETSC_OWN_POINTER,icol);
105:   } else SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"MatRestoreRowIJ fails!");
106:   ISSetPermutation(*irow);
107:   ISSetPermutation(*icol);
108:   return(0);
109: }




114: /*TEST

116:    test:

118: TEST*/