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 Parameters:
13: . L - The graph Laplacian matrix
15: Level: intermediate
17: .seealso: `MatChop()`, `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;
28: MatCreate(PetscObjectComm((PetscObject)A), L);
29: MatGetSize(A, &M, &N);
30: MatGetLocalSize(A, &m, &n);
31: MatSetSizes(*L, m, n, M, N);
32: MatGetOwnershipRange(A, &rStart, &rEnd);
33: PetscMalloc2(m, &dnnz, m, &onnz);
34: for (r = rStart; r < rEnd; ++r) {
35: const PetscScalar *vals;
36: const PetscInt *cols;
37: PetscInt ncols, newcols, c;
38: PetscBool hasdiag = PETSC_FALSE;
40: dnnz[r - rStart] = onnz[r - rStart] = 0;
41: MatGetRow(A, r, &ncols, &cols, &vals);
42: for (c = 0, newcols = 0; c < ncols; ++c) {
43: if (cols[c] == r) {
44: ++newcols;
45: hasdiag = PETSC_TRUE;
46: ++dnnz[r - rStart];
47: } else if (PetscAbsScalar(vals[c]) >= tol) {
48: if ((cols[c] >= rStart) && (cols[c] < rEnd)) ++dnnz[r - rStart];
49: else ++onnz[r - rStart];
50: ++newcols;
51: }
52: }
53: if (!hasdiag) {
54: ++newcols;
55: ++dnnz[r - rStart];
56: }
57: colMax = PetscMax(colMax, newcols);
58: MatRestoreRow(A, r, &ncols, &cols, &vals);
59: }
60: MatSetFromOptions(*L);
61: MatXAIJSetPreallocation(*L, 1, dnnz, onnz, NULL, NULL);
62: MatSetUp(*L);
63: PetscMalloc2(colMax, &newCols, colMax, &newVals);
64: for (r = rStart; r < rEnd; ++r) {
65: const PetscScalar *vals;
66: const PetscInt *cols;
67: PetscInt ncols, newcols, c;
68: PetscBool hasdiag = PETSC_FALSE;
70: MatGetRow(A, r, &ncols, &cols, &vals);
71: for (c = 0, newcols = 0; c < ncols; ++c) {
72: if (cols[c] == r) {
73: newCols[newcols] = cols[c];
74: newVals[newcols] = dnnz[r - rStart] + onnz[r - rStart] - 1;
75: ++newcols;
76: hasdiag = PETSC_TRUE;
77: } else if (PetscAbsScalar(vals[c]) >= tol) {
78: newCols[newcols] = cols[c];
79: newVals[newcols] = -1.0;
80: ++newcols;
81: }
83: }
84: if (!hasdiag) {
85: newCols[newcols] = r;
86: newVals[newcols] = dnnz[r - rStart] + onnz[r - rStart] - 1;
87: ++newcols;
88: }
89: MatRestoreRow(A, r, &ncols, &cols, &vals);
90: MatSetValues(*L, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
91: }
92: PetscFree2(dnnz, onnz);
93: MatAssemblyBegin(*L, MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(*L, MAT_FINAL_ASSEMBLY);
95: PetscFree2(newCols, newVals);
96: return 0;
97: }
99: /*
100: MatGetOrdering_Spectral - Find the symmetric reordering of the graph by .
101: */
102: PETSC_INTERN PetscErrorCode MatGetOrdering_Spectral(Mat A, MatOrderingType type, IS *row, IS *col)
103: {
104: Mat L;
105: const PetscReal eps = 1.0e-12;
107: MatCreateLaplacian(A, eps, PETSC_FALSE, &L);
108: {
109: /* Check Laplacian */
110: PetscReal norm;
111: Vec x, y;
113: MatCreateVecs(L, &x, NULL);
114: VecDuplicate(x, &y);
115: VecSet(x, 1.0);
116: MatMult(L, x, y);
117: VecNorm(y, NORM_INFINITY, &norm);
119: VecDestroy(&x);
120: VecDestroy(&y);
121: }
122: /* Compute Fiedler vector (right now, all eigenvectors) */
123: #if defined(PETSC_USE_COMPLEX)
124: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Spectral partitioning does not support complex numbers");
125: #else
126: {
127: Mat LD;
128: PetscScalar *a;
129: PetscReal *realpart, *imagpart, *eigvec, *work;
130: PetscReal sdummy;
131: PetscBLASInt bn, bN, lwork = 0, lierr, idummy;
132: PetscInt n, i, evInd, *perm, tmp;
134: MatConvert(L, MATDENSE, MAT_INITIAL_MATRIX, &LD);
135: MatGetLocalSize(LD, &n, NULL);
136: MatDenseGetArray(LD, &a);
137: PetscBLASIntCast(n, &bn);
138: PetscBLASIntCast(n, &bN);
139: PetscBLASIntCast(5 * n, &lwork);
140: PetscBLASIntCast(1, &idummy);
141: PetscMalloc4(n, &realpart, n, &imagpart, n * n, &eigvec, lwork, &work);
142: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
143: PetscCallBLAS("LAPACKgeev", LAPACKgeev_("N", "V", &bn, a, &bN, realpart, imagpart, &sdummy, &idummy, eigvec, &bN, work, &lwork, &lierr));
145: PetscFPTrapPop();
146: MatDenseRestoreArray(LD, &a);
147: MatDestroy(&LD);
148: /* Check lowest eigenvalue and eigenvector */
149: PetscMalloc1(n, &perm);
150: for (i = 0; i < n; ++i) perm[i] = i;
151: PetscSortRealWithPermutation(n, realpart, perm);
152: evInd = perm[0];
154: evInd = perm[1];
156: evInd = perm[0];
157: for (i = 0; i < n; ++i) {
159: }
160: /* Construct Fiedler partition */
161: evInd = perm[1];
162: for (i = 0; i < n; ++i) perm[i] = i;
163: PetscSortRealWithPermutation(n, &eigvec[evInd * n], perm);
164: for (i = 0; i < n / 2; ++i) {
165: tmp = perm[n - 1 - i];
166: perm[n - 1 - i] = perm[i];
167: perm[i] = tmp;
168: }
169: ISCreateGeneral(PETSC_COMM_SELF, n, perm, PETSC_OWN_POINTER, row);
170: PetscObjectReference((PetscObject)*row);
171: *col = *row;
173: PetscFree4(realpart, imagpart, eigvec, work);
174: MatDestroy(&L);
175: return 0;
176: }
177: #endif
178: }