Actual source code: spectral.c
petsc-3.6.1 2015-08-06
1: #include <petscmat.h> /*I <petscmat.h> I*/
2: #include <petscblaslapack.h>
6: /*@
7: MatCreateLaplacian - Create the matrix Laplacian, with all values in the matrix less than the tolerance set to zero
9: Input Parameters:
10: + A - The matrix
11: . tol - The zero tolerance
12: - weighted - Flag for using edge weights
14: Output Parameters:
15: . L - The graph Laplacian matrix
17: Level: intermediate
19: .seealso: MatChop()
20: @*/
21: PetscErrorCode MatCreateLaplacian(Mat A, PetscReal tol, PetscBool weighted, Mat *L)
22: {
23: PetscScalar *newVals;
24: PetscInt *newCols;
25: PetscInt rStart, rEnd, r, colMax = 0;
26: PetscInt *dnnz, *onnz;
27: PetscInt m, n, M, N;
31: if (weighted) SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "Will get to this soon");
32: MatCreate(PetscObjectComm((PetscObject) A), L);
33: MatGetSize(A, &M, &N);
34: MatGetLocalSize(A, &m, &n);
35: MatSetSizes(*L, m, n, M, N);
36: MatGetOwnershipRange(A, &rStart, &rEnd);
37: PetscMalloc2(m,&dnnz,m,&onnz);
38: for (r = rStart; r < rEnd; ++r) {
39: const PetscScalar *vals;
40: const PetscInt *cols;
41: PetscInt ncols, newcols, c;
42: PetscBool hasdiag = PETSC_FALSE;
44: dnnz[r-rStart] = onnz[r-rStart] = 0;
45: MatGetRow(A, r, &ncols, &cols, &vals);
46: for (c = 0, newcols = 0; c < ncols; ++c) {
47: if (cols[c] == r) {
48: ++newcols;
49: hasdiag = PETSC_TRUE;
50: ++dnnz[r-rStart];
51: } else if (PetscAbsScalar(vals[c]) >= tol) {
52: if ((cols[c] >= rStart) && (cols[c] < rEnd)) ++dnnz[r-rStart];
53: else ++onnz[r-rStart];
54: ++newcols;
55: }
56: }
57: if (!hasdiag) {++newcols; ++dnnz[r-rStart];}
58: colMax = PetscMax(colMax, newcols);
59: MatRestoreRow(A, r, &ncols, &cols, &vals);
60: }
61: MatSetFromOptions(*L);
62: MatXAIJSetPreallocation(*L, 1, dnnz, onnz, NULL, NULL);
63: MatSetUp(*L);
64: 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: 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: if (newcols > colMax) SETERRQ(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: MatRestoreRow(A, r, &ncols, &cols, &vals);
91: MatSetValues(*L, 1, &r, newcols, newCols, newVals, INSERT_VALUES);
92: }
93: PetscFree2(dnnz,onnz);
94: MatAssemblyBegin(*L, MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(*L, MAT_FINAL_ASSEMBLY);
96: PetscFree2(newCols,newVals);
97: return(0);
98: }
100: /*
101: MatGetOrdering_Spectral - Find the symmetric reordering of the graph by .
102: */
105: PETSC_EXTERN PetscErrorCode MatGetOrdering_Spectral(Mat A, MatOrderingType type, IS *row, IS *col)
106: {
107: Mat L;
108: PetscInt *perm, tmp;
109: const PetscReal eps = 1.0e-12;
110: PetscErrorCode ierr;
113: MatCreateLaplacian(A, eps, PETSC_FALSE, &L);
114: {
115: /* Check Laplacian */
116: PetscReal norm;
117: Vec x, y;
119: MatCreateVecs(L, &x, NULL);
120: VecDuplicate(x, &y);
121: VecSet(x, 1.0);
122: MatMult(L, x, y);
123: VecNorm(y, NORM_INFINITY, &norm);
124: if (norm > 1.0e-10) SETERRQ(PetscObjectComm((PetscObject) y), PETSC_ERR_PLIB, "Invalid graph Laplacian");
125: VecDestroy(&x);
126: VecDestroy(&y);
127: }
128: /* Compute Fiedler vector (right now, all eigenvectors) */
129: {
130: Mat LD;
131: PetscScalar *a;
132: PetscReal *realpart, *imagpart, *eigvec, *work;
133: #ifndef PETSC_USE_COMPLEX
134: PetscReal sdummy;
135: #endif
136: PetscBLASInt bn, bN, lwork, lierr, idummy;
137: PetscInt n, i, evInd;
139: MatConvert(L, MATDENSE, MAT_INITIAL_MATRIX, &LD);
140: MatGetLocalSize(LD, &n, NULL);
141: MatDenseGetArray(LD, &a);
142: PetscBLASIntCast(n, &bn);
143: PetscBLASIntCast(n, &bN);
144: PetscBLASIntCast(5*n,&lwork);
145: PetscBLASIntCast(1,&idummy);
146: PetscMalloc4(n,&realpart,n,&imagpart,n*n,&eigvec,lwork,&work);
147: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
148: #ifdef PETSC_USE_COMPLEX
149: SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "Spectral partitioning does not support complex numbers");
150: #elif defined(PETSC_HAVE_ESSL)
151: SETERRQ(PetscObjectComm((PetscObject) A), PETSC_ERR_SUP, "Spectral partitioning does not support ESSL Lapack Routines");
152: #else
153: PetscStackCall("LAPACKgeev", LAPACKgeev_("N","V",&bn,a,&bN,realpart,imagpart,&sdummy,&idummy,eigvec,&bN,work,&lwork,&lierr));
154: #endif
155: if (lierr) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in LAPACK routine %d", (int) lierr);
156: PetscFPTrapPop();
157: MatDenseRestoreArray(LD,&a);
158: MatDestroy(&LD);
159: /* Check lowest eigenvalue and eigenvector */
160: PetscMalloc1(n, &perm);
161: for (i = 0; i < n; ++i) perm[i] = i;
162: PetscSortRealWithPermutation(n,realpart,perm);
163: evInd = perm[0];
164: 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");
165: evInd = perm[1];
166: 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");
167: evInd = perm[0];
168: for (i = 0; i < n; ++i) {
169: if (fabs(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]);
170: }
171: /* Construct Fiedler partition */
172: evInd = perm[1];
173: for (i = 0; i < n; ++i) perm[i] = i;
174: PetscSortRealWithPermutation(n, &eigvec[evInd*n], perm);
175: for (i = 0; i < n/2; ++i) {
176: tmp = perm[n-1-i];
177: perm[n-1-i] = perm[i];
178: perm[i] = tmp;
179: }
180: ISCreateGeneral(PETSC_COMM_SELF, n, perm, PETSC_OWN_POINTER, row);
181: PetscObjectReference((PetscObject) *row);
182: *col = *row;
184: PetscFree4(realpart,imagpart,eigvec,work);
185: MatDestroy(&L);
186: }
187: return(0);
188: }