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: }