Actual source code: spectral.c

petsc-3.10.5 2019-03-28
Report Typos and Errors
  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: }