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 @*/
MatCreateLaplacian(Mat A,PetscReal tol,PetscBool weighted,Mat * L)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 */
MatGetOrdering_Spectral(Mat A,MatOrderingType type,IS * row,IS * col)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 %" PetscBLASInt_FMT, 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