1 #include <petscmat.h>
2 #include <petsc/private/matorderimpl.h>
3
4 #if defined(PETSC_HAVE_SUPERLU_DIST)
5
6 /* SuperLU_DIST bundles f2ced mc64ad_() from HSL */
7
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 */
12
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
17
18 #elif !defined(PETSC_BLASLAPACK_UNDERSCORE)
19 #define mc64id_dist mc64id_dist
20 #define mc64ad_dist mc64ad_dist
21
22 #endif
23
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
27
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 */
MatGetOrdering_WBM(Mat mat,MatOrderingType type,IS * row,IS * col)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
42
43 PetscFunctionBegin;
44 PetscCall(MatGetRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, &nrow, &ia, &ja, &done));
45 ncol = nrow;
46 nnz = ia[nrow];
47 PetscCheck(done, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot get rows for matrix");
48 PetscCall(MatSeqAIJGetArray(mat, &a));
49 switch (job) {
50 case 1:
51 liw = 4 * nrow + ncol;
52 ldw = 0;
53 break;
54 case 2:
55 liw = 2 * nrow + 2 * ncol;
56 ldw = ncol;
57 break;
58 case 3:
59 liw = 8 * nrow + 2 * ncol + nnz;
60 ldw = nnz;
61 break;
62 case 4:
63 liw = 3 * nrow + 2 * ncol;
64 ldw = 2 * ncol + nnz;
65 break;
66 case 5:
67 liw = 3 * nrow + 2 * ncol;
68 ldw = nrow + 2 * ncol + nnz;
69 break;
70 }
71
72 PetscCall(PetscMalloc3(liw, &iw, ldw, &dw, nrow, &perm));
73 #if defined(PETSC_HAVE_SUPERLU_DIST)
74 PetscCallExternal(mc64id_dist, icntl);
75 icntl[0] = 0; /* allow printing error messages (f2c'd code uses if non-negative, ignores value otherwise) */
76 icntl[1] = -1; /* suppress warnings */
77 icntl[2] = -1; /* ignore diagnostic output [default] */
78 icntl[3] = 0; /* perform consistency checks [default] */
79 PetscCallExternal(mc64ad_dist, &job, &nrow, &nnz, ia, ja, a, &num, perm, &liw, iw, &ldw, dw, icntl, info);
80 PetscCall(MatRestoreRowIJ(mat, 1, PETSC_TRUE, PETSC_TRUE, NULL, &ia, &ja, &done));
81 for (i = 0; i < nrow; ++i) perm[i]--;
82 /* If job == 5, dw[0..ncols] contains the column scaling and dw[ncols..ncols+nrows] contains the row scaling */
83 PetscCall(ISCreateStride(PETSC_COMM_SELF, nrow, 0, 1, row));
84 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nrow, perm, PETSC_COPY_VALUES, col));
85 PetscCall(PetscFree3(iw, dw, perm));
86 PetscFunctionReturn(PETSC_SUCCESS);
87 #else
88 SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "WBM using MC64 does not support complex numbers");
89 #endif
90 }
91