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