1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: static char help[] = "Computes a GSVD associated with a complex network analysis problem.\n"
 12:   "The command line options are:\n"
 13:   "  -file <filename>, text file in TSV format containing undirected graph.\n"
 14:   "  -beta <beta>, optional decay parameter.\n\n";

 16: /*
 17:     Computes the Katz high-order proximity embedding of a graph via a partial GSVD.

 19:     [1] M. Ou et al. "Asymmetric transitivity preserving graph embedding" Proc. of
 20:         ACM SIGKDD 2016 - https://doi.org/10.1145/2939672.2939751
 21: */

 23: #include <slepcsvd.h>
 24: #include "network.h"

 26: PetscErrorCode SpectralRadius(Mat A,PetscReal *rho)
 27: {
 28:   EPS         eps;
 29:   PetscInt    nconv;
 30:   PetscScalar kr,ki;

 32:   PetscFunctionBeginUser;
 33:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
 34:   PetscCall(EPSSetOperators(eps,A,NULL));
 35:   PetscCall(EPSSetProblemType(eps,EPS_NHEP));
 36:   PetscCall(EPSSetDimensions(eps,1,PETSC_DEFAULT,PETSC_DEFAULT));
 37:   PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_MAGNITUDE));
 38:   PetscCall(EPSSetFromOptions(eps));
 39:   PetscCall(EPSSolve(eps));
 40:   PetscCall(EPSGetConverged(eps,&nconv));
 41:   *rho = 0.0;
 42:   if (nconv) {
 43:     PetscCall(EPSGetEigenvalue(eps,0,&kr,&ki));
 44:     *rho = SlepcAbsEigenvalue(kr,ki);
 45:   }
 46:   PetscCall(EPSDestroy(&eps));
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: int main(int argc,char **argv)
 51: {
 52:   Graph          graph;           /* directed graph */
 53:   Mat            A,B;             /* matrices */
 54:   SVD            svd;             /* singular value problem solver context */
 55:   PetscReal      beta=0.0,rho=0.0;
 56:   char           filename[PETSC_MAX_PATH_LEN];
 57:   PetscBool      flg,terse;

 59:   PetscFunctionBeginUser;
 60:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

 62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63:         Load graph and build adjacency matrix
 64:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 66:   PetscCall(PetscOptionsGetString(NULL,NULL,"-graph",filename,sizeof(filename),&flg));
 67:   PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must use -graph <name> to indicate a file containing the graph in TSV format");

 69:   PetscCall(GraphCreate(PETSC_COMM_WORLD,&graph));
 70:   PetscCall(GraphPreload(graph,filename));
 71:   PetscCheck(graph->type==GRAPH_DIRECTED,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"A directed graph is required");
 72:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGraph with %" PetscInt_FMT " vertices and %" PetscInt_FMT " edges.\n\n",graph->nvertices,graph->nedges));
 73:   if (graph->weight!=GRAPH_WEIGHT_UNWEIGHTED) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"WARNING: ignoring weights in input graph\n"));

 75:   PetscCall(GraphPreallocate(graph,filename));
 76:   PetscCall(GraphLoadUnweighted(graph,filename));

 78:   PetscCall(MatViewFromOptions(graph->adjacency,NULL,"-adjacency"));

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:       Build high-order proximity matrices (Katz): A=I-beta*G, B=beta*G
 82:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 84:   PetscCall(PetscOptionsGetReal(NULL,NULL,"-beta",&beta,NULL));
 85:   if (beta==0.0) {  /* compute decay parameter beta */
 86:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Computing spectral radius...\n"));
 87:     PetscCall(SpectralRadius(graph->adjacency,&rho));
 88:     if (rho) {
 89:       beta = 0.8*rho;
 90:       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Done, rho=%g, setting decay parameter beta=%g\n\n",(double)rho,(double)beta));
 91:     } else {
 92:       beta = 0.5;
 93:       PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Failed, using default decay parameter beta=%g\n\n",(double)beta));
 94:     }
 95:   }
 96:   PetscCall(MatDuplicate(graph->adjacency,MAT_COPY_VALUES,&A));
 97:   PetscCall(MatScale(A,-beta));
 98:   PetscCall(MatShift(A,1.0));
 99:   PetscCall(MatDuplicate(graph->adjacency,MAT_COPY_VALUES,&B));
100:   PetscCall(MatScale(B,beta));

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:            Create the singular value solver and set various options
104:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   PetscCall(SVDCreate(PETSC_COMM_WORLD,&svd));
107:   PetscCall(SVDSetOperators(svd,A,B));
108:   PetscCall(SVDSetProblemType(svd,SVD_GENERALIZED));
109:   PetscCall(SVDSetType(svd,SVDTRLANCZOS));
110:   PetscCall(SVDSetFromOptions(svd));

112:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113:                       Solve the problem and print solution
114:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

116:   PetscCall(SVDSolve(svd));

118:   /* show detailed info unless -terse option is given by user */
119:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
120:   if (terse) PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,NULL));
121:   else {
122:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
123:     PetscCall(SVDConvergedReasonView(svd,PETSC_VIEWER_STDOUT_WORLD));
124:     PetscCall(SVDErrorView(svd,SVD_ERROR_NORM,PETSC_VIEWER_STDOUT_WORLD));
125:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
126:   }
127:   PetscCall(SVDDestroy(&svd));
128:   PetscCall(GraphDestroy(&graph));
129:   PetscCall(MatDestroy(&A));
130:   PetscCall(MatDestroy(&B));
131:   PetscCall(SlepcFinalize());
132:   return 0;
133: }
134: /*TEST

136:    build:
137:       depends: network.o

139:    test:
140:       args: -graph ${SLEPC_DIR}/share/slepc/datafiles/graphs/out.moreno_taro_taro -svd_nsv 4 -terse
141:       filter: sed -e 's/4.38031/4.38032/' | sed -e 's/3.75089/3.7509/' | sed -e 's/3.00071/3.00072/'

143: TEST*/