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