Actual source code: wbm.c

  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 include superlu_dist includes? */
 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 *, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscScalar *, PetscInt *, PetscInt *);
 26: #endif

 28: /*
 29:   MatGetOrdering_WBM - Find the nonsymmetric reordering of the graph which maximizes the product of diagonal entries,
 30:     using weighted bipartite graph matching. This is MC64 in the Harwell-Boeing library.
 31: */
 32: PETSC_INTERN PetscErrorCode MatGetOrdering_WBM(Mat mat, MatOrderingType type, IS *row, IS *col)
 33: {
 34:   PetscScalar    *a, *dw;
 35:   const PetscInt *ia, *ja;
 36:   PetscInt        job = 5;
 37:   PetscInt       *perm, nrow, ncol, nnz, liw, *iw, ldw;
 38:   PetscBool       done;
 39: #if defined(PETSC_HAVE_SUPERLU_DIST)
 40:   PetscInt num, info[10], icntl[10], i;
 41: #endif

 43:   MatGetRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done);
 44:   ncol = nrow;
 45:   nnz  = ia[nrow];
 47:   MatSeqAIJGetArray(mat, &a);
 48:   switch (job) {
 49:   case 1:
 50:     liw = 4 * nrow + ncol;
 51:     ldw = 0;
 52:     break;
 53:   case 2:
 54:     liw = 2 * nrow + 2 * ncol;
 55:     ldw = ncol;
 56:     break;
 57:   case 3:
 58:     liw = 8 * nrow + 2 * ncol + nnz;
 59:     ldw = nnz;
 60:     break;
 61:   case 4:
 62:     liw = 3 * nrow + 2 * ncol;
 63:     ldw = 2 * ncol + nnz;
 64:     break;
 65:   case 5:
 66:     liw = 3 * nrow + 2 * ncol;
 67:     ldw = nrow + 2 * ncol + nnz;
 68:     break;
 69:   }

 71:   PetscMalloc3(liw, &iw, ldw, &dw, nrow, &perm);
 72: #if defined(PETSC_HAVE_SUPERLU_DIST)
 73:   PetscCallExternal(mc64id_dist, icntl);
 74:   icntl[0] = 0;  /* allow printing error messages (f2c'd code uses if non-negative, ignores value otherwise) */
 75:   icntl[1] = -1; /* suppress warnings */
 76:   icntl[2] = -1; /* ignore diagnostic output [default] */
 77:   icntl[3] = 0;  /* perform consistency checks [default] */
 78:   PetscCallExternal(mc64ad_dist, &job, &nrow, &nnz, ia, ja, a, &num, perm, &liw, iw, &ldw, dw, icntl, info);
 79:   MatRestoreRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done);
 80:   for (i = 0; i < nrow; ++i) perm[i]--;
 81:   /* If job == 5, dw[0..ncols] contains the column scaling and dw[ncols..ncols+nrows] contains the row scaling */
 82:   ISCreateStride(PETSC_COMM_SELF, nrow, 0, 1, row);
 83:   ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, col);
 84:   PetscFree3(iw, dw, perm);
 85:   return 0;
 86: #else
 87:   SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "WBM using MC64 does not support complex numbers");
 88: #endif
 89: }