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: }