Actual source code: wbm.c
petsc-3.6.1 2015-08-06
1: #include <petscmat.h>
2: #include <petsc/private/matorderimpl.h>
4: #if defined(PETSC_HAVE_SUPERLU_DIST)
6: /* SuperLU_DIST bundles f2ced mc64ad_() from HSL */
8: /*
9: SuperLU_dist uses a common flag for both Fortran mangling and BLAS/LAPACK mangling which
10: corresponds to the PETSc BLAS/LAPACK mangling flag (we pass this flag to configure SuperLU_dist)
11: */
13: /* Why not incude superlu_dist inludes? */
14: # if defined(PETSC_BLASLAPACK_CAPS)
15: # define mc64id_dist MC64ID_DIST
16: # define mc64ad_dist MC64AD_DIST
18: # elif !defined(PETSC_BLASLAPACK_UNDERSCORE)
19: # define mc64id_dist mc64id_dist
20: # define mc64ad_dist mc64ad_dist
22: # endif
24: PETSC_EXTERN PetscInt mc64id_dist(PetscInt*);
25: PETSC_EXTERN PetscInt mc64ad_dist(const PetscInt*, PetscInt*, PetscInt*, const PetscInt*, const PetscInt*n, PetscScalar*, PetscInt*,
26: PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscScalar*, PetscInt*, PetscInt*);
27: #endif
29: /*
30: MatGetOrdering_WBM - Find the nonsymmetric reordering of the graph which maximizes the product of diagonal entries,
31: using weighted bipartite graph matching. This is MC64 in the Harwell-Boeing library.
32: */
35: PETSC_EXTERN PetscErrorCode MatGetOrdering_WBM(Mat mat, MatOrderingType type, IS *row, IS *col)
36: {
37: PetscScalar *a, *dw;
38: const PetscInt *ia, *ja;
39: const PetscInt job = 5;
40: PetscInt *perm, nrow, ncol, nnz, liw, *iw, ldw, i;
41: PetscBool done;
42: PetscErrorCode ierr;
45: MatGetRowIJ(mat,1,PETSC_TRUE,PETSC_TRUE,&nrow,&ia,&ja,&done);
46: ncol = nrow;
47: nnz = ia[nrow];
48: if (!done) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get rows for matrix");
49: MatSeqAIJGetArray(mat, &a);
50: switch (job) {
51: case 1: liw = 4*nrow + ncol; ldw = 0;break;
52: case 2: liw = 2*nrow + 2*ncol; ldw = ncol;break;
53: case 3: liw = 8*nrow + 2*ncol + nnz; ldw = nnz;break;
54: case 4: liw = 3*nrow + 2*ncol; ldw = 2*ncol + nnz;break;
55: case 5: liw = 3*nrow + 2*ncol; ldw = nrow + 2*ncol + nnz;break;
56: }
58: PetscMalloc3(liw,&iw,ldw,&dw,nrow,&perm);
59: #if defined(PETSC_HAVE_SUPERLU_DIST)
60: {
61: PetscInt num, info[10], icntl[10];
63: mc64id_dist(icntl);
64: if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"HSL mc64id_dist returned %d\n",ierr);
65: icntl[0] = 0; /* allow printing error messages (f2c'd code uses if non-negative, ignores value otherwise) */
66: icntl[1] = -1; /* suppress warnings */
67: icntl[2] = -1; /* ignore diagnostic output [default] */
68: icntl[3] = 0; /* perform consistency checks [default] */
69: mc64ad_dist(&job, &nrow, &nnz, ia, ja, a, &num, perm, &liw, iw, &ldw, dw, icntl, info);
70: if (ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"HSL mc64ad_dist returned %d\n",ierr);
71: }
72: #else
73: SETERRQ(PetscObjectComm((PetscObject) mat), PETSC_ERR_SUP, "WBM using MC64 does not support complex numbers");
74: #endif
75: MatRestoreRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done);
76: for (i = 0; i < nrow; ++i) perm[i]--;
77: /* If job == 5, dw[0..ncols] contains the column scaling and dw[ncols..ncols+nrows] contains the row scaling */
78: ISCreateStride(PETSC_COMM_SELF, nrow, 0, 1, row);
79: ISCreateGeneral(PETSC_COMM_SELF,nrow,perm,PETSC_COPY_VALUES,col);
80: PetscFree3(iw,dw,perm);
81: return(0);
82: }