Actual source code: spectral.c
1: #include <petscmat.h>
2: #include <petscblaslapack.h>
4: /*@
5: MatCreateLaplacian - Create the matrix Laplacian, with all values in the matrix less than the tolerance set to zero
7: Input Parameters:
8: + A - The matrix
9: . tol - The zero tolerance
10: - weighted - Flag for using edge weights
12: Output Parameter:
13: . L - The graph Laplacian matrix
15: Level: intermediate
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;
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;
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;
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: }
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;
108: PetscFunctionBegin;
109: PetscCall(MatCreateLaplacian(A, eps, PETSC_FALSE, &L));
110: {
111: /* Check Laplacian */
112: PetscReal norm;
113: Vec x, y;
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;
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;
175: PetscCall(PetscFree4(realpart, imagpart, eigvec, work));
176: PetscCall(MatDestroy(&L));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
179: #endif
180: }