1 #include <petscmat.h> /*I <petscmat.h> I*/ 2 #include <petscblaslapack.h> 3 4 /*@ 5 MatCreateLaplacian - Create the matrix Laplacian, with all values in the matrix less than the tolerance set to zero 6 7 Input Parameters: 8 + A - The matrix 9 . tol - The zero tolerance 10 - weighted - Flag for using edge weights 11 12 Output Parameter: 13 . L - The graph Laplacian matrix 14 15 Level: intermediate 16 17 .seealso: `MatFilter()`, `MatGetGraph()` 18 @*/ 19 PetscErrorCode MatCreateLaplacian(Mat A, PetscReal tol, PetscBool weighted, Mat *L) 20 { 21 PetscScalar *newVals; 22 PetscInt *newCols; 23 PetscInt rStart, rEnd, r, colMax = 0; 24 PetscInt *dnnz, *onnz; 25 PetscInt m, n, M, N; 26 27 PetscFunctionBegin; 28 PetscCheck(!weighted, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Will get to this soon"); 29 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), L)); 30 PetscCall(MatGetSize(A, &M, &N)); 31 PetscCall(MatGetLocalSize(A, &m, &n)); 32 PetscCall(MatSetSizes(*L, m, n, M, N)); 33 PetscCall(MatGetOwnershipRange(A, &rStart, &rEnd)); 34 PetscCall(PetscMalloc2(m, &dnnz, m, &onnz)); 35 for (r = rStart; r < rEnd; ++r) { 36 const PetscScalar *vals; 37 const PetscInt *cols; 38 PetscInt ncols, newcols, c; 39 PetscBool hasdiag = PETSC_FALSE; 40 41 dnnz[r - rStart] = onnz[r - rStart] = 0; 42 PetscCall(MatGetRow(A, r, &ncols, &cols, &vals)); 43 for (c = 0, newcols = 0; c < ncols; ++c) { 44 if (cols[c] == r) { 45 ++newcols; 46 hasdiag = PETSC_TRUE; 47 ++dnnz[r - rStart]; 48 } else if (PetscAbsScalar(vals[c]) >= tol) { 49 if ((cols[c] >= rStart) && (cols[c] < rEnd)) ++dnnz[r - rStart]; 50 else ++onnz[r - rStart]; 51 ++newcols; 52 } 53 } 54 if (!hasdiag) { 55 ++newcols; 56 ++dnnz[r - rStart]; 57 } 58 colMax = PetscMax(colMax, newcols); 59 PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals)); 60 } 61 PetscCall(MatSetFromOptions(*L)); 62 PetscCall(MatXAIJSetPreallocation(*L, 1, dnnz, onnz, NULL, NULL)); 63 PetscCall(MatSetUp(*L)); 64 PetscCall(PetscMalloc2(colMax, &newCols, colMax, &newVals)); 65 for (r = rStart; r < rEnd; ++r) { 66 const PetscScalar *vals; 67 const PetscInt *cols; 68 PetscInt ncols, newcols, c; 69 PetscBool hasdiag = PETSC_FALSE; 70 71 PetscCall(MatGetRow(A, r, &ncols, &cols, &vals)); 72 for (c = 0, newcols = 0; c < ncols; ++c) { 73 if (cols[c] == r) { 74 newCols[newcols] = cols[c]; 75 newVals[newcols] = dnnz[r - rStart] + onnz[r - rStart] - 1; 76 ++newcols; 77 hasdiag = PETSC_TRUE; 78 } else if (PetscAbsScalar(vals[c]) >= tol) { 79 newCols[newcols] = cols[c]; 80 newVals[newcols] = -1.0; 81 ++newcols; 82 } 83 PetscCheck(newcols <= colMax, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Overran work space"); 84 } 85 if (!hasdiag) { 86 newCols[newcols] = r; 87 newVals[newcols] = dnnz[r - rStart] + onnz[r - rStart] - 1; 88 ++newcols; 89 } 90 PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals)); 91 PetscCall(MatSetValues(*L, 1, &r, newcols, newCols, newVals, INSERT_VALUES)); 92 } 93 PetscCall(PetscFree2(dnnz, onnz)); 94 PetscCall(MatAssemblyBegin(*L, MAT_FINAL_ASSEMBLY)); 95 PetscCall(MatAssemblyEnd(*L, MAT_FINAL_ASSEMBLY)); 96 PetscCall(PetscFree2(newCols, newVals)); 97 PetscFunctionReturn(PETSC_SUCCESS); 98 } 99 100 /* 101 MatGetOrdering_Spectral - Find the symmetric reordering of the graph by . 102 */ 103 PETSC_INTERN PetscErrorCode MatGetOrdering_Spectral(Mat A, MatOrderingType type, IS *row, IS *col) 104 { 105 Mat L; 106 const PetscReal eps = 1.0e-12; 107 108 PetscFunctionBegin; 109 PetscCall(MatCreateLaplacian(A, eps, PETSC_FALSE, &L)); 110 { 111 /* Check Laplacian */ 112 PetscReal norm; 113 Vec x, y; 114 115 PetscCall(MatCreateVecs(L, &x, NULL)); 116 PetscCall(VecDuplicate(x, &y)); 117 PetscCall(VecSet(x, 1.0)); 118 PetscCall(MatMult(L, x, y)); 119 PetscCall(VecNorm(y, NORM_INFINITY, &norm)); 120 PetscCheck(norm <= 1.0e-10, PetscObjectComm((PetscObject)y), PETSC_ERR_PLIB, "Invalid graph Laplacian"); 121 PetscCall(VecDestroy(&x)); 122 PetscCall(VecDestroy(&y)); 123 } 124 /* Compute Fiedler vector (right now, all eigenvectors) */ 125 #if defined(PETSC_USE_COMPLEX) 126 SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Spectral partitioning does not support complex numbers"); 127 #else 128 { 129 Mat LD; 130 PetscScalar *a; 131 PetscReal *realpart, *imagpart, *eigvec, *work; 132 PetscReal sdummy; 133 PetscBLASInt bn, bN, lwork = 0, lierr, idummy; 134 PetscInt n, i, evInd, *perm, tmp; 135 136 PetscCall(MatConvert(L, MATDENSE, MAT_INITIAL_MATRIX, &LD)); 137 PetscCall(MatGetLocalSize(LD, &n, NULL)); 138 PetscCall(MatDenseGetArray(LD, &a)); 139 PetscCall(PetscBLASIntCast(n, &bn)); 140 PetscCall(PetscBLASIntCast(n, &bN)); 141 PetscCall(PetscBLASIntCast(5 * n, &lwork)); 142 PetscCall(PetscBLASIntCast(1, &idummy)); 143 PetscCall(PetscMalloc4(n, &realpart, n, &imagpart, n * n, &eigvec, lwork, &work)); 144 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 145 PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "V", &bn, a, &bN, realpart, imagpart, &sdummy, &idummy, eigvec, &bN, work, &lwork, &lierr)); 146 PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int)lierr); 147 PetscCall(PetscFPTrapPop()); 148 PetscCall(MatDenseRestoreArray(LD, &a)); 149 PetscCall(MatDestroy(&LD)); 150 /* Check lowest eigenvalue and eigenvector */ 151 PetscCall(PetscMalloc1(n, &perm)); 152 for (i = 0; i < n; ++i) perm[i] = i; 153 PetscCall(PetscSortRealWithPermutation(n, realpart, perm)); 154 evInd = perm[0]; 155 PetscCheck(realpart[evInd] <= 1.0e-12 && imagpart[evInd] <= 1.0e-12, PetscObjectComm((PetscObject)L), PETSC_ERR_PLIB, "Graph Laplacian must have lowest eigenvalue 0"); 156 evInd = perm[1]; 157 PetscCheck(realpart[evInd] >= 1.0e-12 || imagpart[evInd] >= 1.0e-12, PetscObjectComm((PetscObject)L), PETSC_ERR_PLIB, "Graph Laplacian must have only one zero eigenvalue"); 158 evInd = perm[0]; 159 for (i = 0; i < n; ++i) { 160 PetscCheck(PetscAbsReal(eigvec[evInd * n + i] - eigvec[evInd * n + 0]) <= 1.0e-10, PetscObjectComm((PetscObject)L), PETSC_ERR_PLIB, "Graph Laplacian must have constant lowest eigenvector ev_%" PetscInt_FMT " %g != ev_0 %g", i, (double)eigvec[evInd * n + i], (double)eigvec[evInd * n + 0]); 161 } 162 /* Construct Fiedler partition */ 163 evInd = perm[1]; 164 for (i = 0; i < n; ++i) perm[i] = i; 165 PetscCall(PetscSortRealWithPermutation(n, &eigvec[evInd * n], perm)); 166 for (i = 0; i < n / 2; ++i) { 167 tmp = perm[n - 1 - i]; 168 perm[n - 1 - i] = perm[i]; 169 perm[i] = tmp; 170 } 171 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, perm, PETSC_OWN_POINTER, row)); 172 PetscCall(PetscObjectReference((PetscObject)*row)); 173 *col = *row; 174 175 PetscCall(PetscFree4(realpart, imagpart, eigvec, work)); 176 PetscCall(MatDestroy(&L)); 177 PetscFunctionReturn(PETSC_SUCCESS); 178 } 179 #endif 180 } 181