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