1
2 #include <petscmat.h>
3 #include <petsc/private/matorderimpl.h>
4
5 /*
6 MatGetOrdering_QMD - Find the Quotient Minimum Degree ordering of a given matrix.
7 */
MatGetOrdering_QMD(Mat mat,MatOrderingType type,IS * row,IS * col)8 PETSC_INTERN PetscErrorCode MatGetOrdering_QMD(Mat mat,MatOrderingType type,IS *row,IS *col)
9 {
10 PetscInt i, *deg,*marker,*rchset,*nbrhd,*qsize,*qlink,nofsub,*iperm,nrow,*perm;
11 PetscErrorCode ierr;
12 const PetscInt *ia,*ja;
13 PetscBool done;
14
15 PetscFunctionBegin;
16 ierr = MatGetRowIJ(mat,1,PETSC_TRUE,PETSC_TRUE,&nrow,&ia,&ja,&done);CHKERRQ(ierr);
17 if (!done) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get rows for matrix");
18
19 ierr = PetscMalloc1(nrow,&perm);CHKERRQ(ierr);
20 ierr = PetscMalloc5(nrow,&iperm,nrow,°,nrow,&marker,nrow,&rchset,nrow,&nbrhd);CHKERRQ(ierr);
21 ierr = PetscMalloc2(nrow,&qsize,nrow,&qlink);CHKERRQ(ierr);
22 /* WARNING - genqmd trashes ja */
23 SPARSEPACKgenqmd(&nrow,ia,ja,perm,iperm,deg,marker,rchset,nbrhd,qsize,qlink,&nofsub);
24 ierr = MatRestoreRowIJ(mat,1,PETSC_TRUE,PETSC_TRUE,NULL,&ia,&ja,&done);CHKERRQ(ierr);
25
26 ierr = PetscFree2(qsize,qlink);CHKERRQ(ierr);
27 ierr = PetscFree5(iperm,deg,marker,rchset,nbrhd);CHKERRQ(ierr);
28 for (i=0; i<nrow; i++) perm[i]--;
29 ierr = ISCreateGeneral(PETSC_COMM_SELF,nrow,perm,PETSC_COPY_VALUES,row);CHKERRQ(ierr);
30 ierr = ISCreateGeneral(PETSC_COMM_SELF,nrow,perm,PETSC_OWN_POINTER,col);CHKERRQ(ierr);
31 PetscFunctionReturn(0);
32 }
33