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()
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) {++newcols; ++dnnz[r-rStart];}
54: colMax = PetscMax(colMax, newcols);
55: MatRestoreRow(A, r, &ncols, &cols, &vals);
56: }
57: MatSetFromOptions(*L);
58: MatXAIJSetPreallocation(*L, 1, dnnz, onnz, NULL, NULL);
59: MatSetUp(*L);
60: PetscMalloc2(colMax,&newCols,colMax,&newVals);
61: for (r = rStart; r < rEnd; ++r) {
62: const PetscScalar *vals;
63: const PetscInt *cols;
64: PetscInt ncols, newcols, c;
65: PetscBool hasdiag = PETSC_FALSE;
67: MatGetRow(A, r, &ncols, &cols, &vals);
68: for (c = 0, newcols = 0; c < ncols; ++c) {
69: if (cols[c] == r) {
70: newCols[newcols] = cols[c];
71: newVals[newcols] = dnnz[r-rStart]+onnz[r-rStart]-1;
72: ++newcols;
73: hasdiag = PETSC_TRUE;
74: } else if (PetscAbsScalar(vals[c]) >= tol) {
75: newCols[newcols] = cols[c];
76: newVals[newcols] = -1.0;
77: ++newcols;
78: }
80: }
81: if (!hasdiag) {
82: newCols[newcols] = r;
83: newVals[newcols] = dnnz[r-rStart]+onnz[r-rStart]-1;
84: ++newcols;
85: }
86: MatRestoreRow(A, r, &ncols, &cols, &vals);
87: MatSetValues(*L, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
88: }
89: PetscFree2(dnnz,onnz);
90: MatAssemblyBegin(*L, MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(*L, MAT_FINAL_ASSEMBLY);
92: PetscFree2(newCols,newVals);
93: return 0;
94: }
96: /*
97: MatGetOrdering_Spectral - Find the symmetric reordering of the graph by .
98: */
99: PETSC_INTERN PetscErrorCode MatGetOrdering_Spectral(Mat A, MatOrderingType type, IS *row, IS *col)
100: {
101: Mat L;
102: const PetscReal eps = 1.0e-12;
104: MatCreateLaplacian(A, eps, PETSC_FALSE, &L);
105: {
106: /* Check Laplacian */
107: PetscReal norm;
108: Vec x, y;
110: MatCreateVecs(L, &x, NULL);
111: VecDuplicate(x, &y);
112: VecSet(x, 1.0);
113: MatMult(L, x, y);
114: VecNorm(y, NORM_INFINITY, &norm);
116: VecDestroy(&x);
117: VecDestroy(&y);
118: }
119: /* Compute Fiedler vector (right now, all eigenvectors) */
120: #if defined(PETSC_USE_COMPLEX)
121: SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "Spectral partitioning does not support complex numbers");
122: #else
123: {
124: Mat LD;
125: PetscScalar *a;
126: PetscReal *realpart, *imagpart, *eigvec, *work;
127: PetscReal sdummy;
128: PetscBLASInt bn, bN, lwork = 0, lierr, idummy;
129: PetscInt n, i, evInd, *perm, tmp;
131: MatConvert(L, MATDENSE, MAT_INITIAL_MATRIX, &LD);
132: MatGetLocalSize(LD, &n, NULL);
133: MatDenseGetArray(LD, &a);
134: PetscBLASIntCast(n, &bn);
135: PetscBLASIntCast(n, &bN);
136: PetscBLASIntCast(5*n,&lwork);
137: PetscBLASIntCast(1,&idummy);
138: PetscMalloc4(n,&realpart,n,&imagpart,n*n,&eigvec,lwork,&work);
139: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
140: PetscStackCallBLAS("LAPACKgeev", LAPACKgeev_("N","V",&bn,a,&bN,realpart,imagpart,&sdummy,&idummy,eigvec,&bN,work,&lwork,&lierr));
142: PetscFPTrapPop();
143: MatDenseRestoreArray(LD,&a);
144: MatDestroy(&LD);
145: /* Check lowest eigenvalue and eigenvector */
146: PetscMalloc1(n, &perm);
147: for (i = 0; i < n; ++i) perm[i] = i;
148: PetscSortRealWithPermutation(n,realpart,perm);
149: evInd = perm[0];
151: evInd = perm[1];
153: evInd = perm[0];
154: for (i = 0; i < n; ++i) {
156: }
157: /* Construct Fiedler partition */
158: evInd = perm[1];
159: for (i = 0; i < n; ++i) perm[i] = i;
160: PetscSortRealWithPermutation(n, &eigvec[evInd*n], perm);
161: for (i = 0; i < n/2; ++i) {
162: tmp = perm[n-1-i];
163: perm[n-1-i] = perm[i];
164: perm[i] = tmp;
165: }
166: ISCreateGeneral(PETSC_COMM_SELF, n, perm, PETSC_OWN_POINTER, row);
167: PetscObjectReference((PetscObject) *row);
168: *col = *row;
170: PetscFree4(realpart,imagpart,eigvec,work);
171: MatDestroy(&L);
172: return 0;
173: }
174: #endif
175: }