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 */ 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