Actual source code: spectral.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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_INTERN 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: }