return D;
}
-
-void SparseMatrix_page_rank(SparseMatrix A, real teleport_probablity, int weighted, real epsilon, real **page_rank){
- /* A(i,j)/Sum_k A(i,k) gives the probablity of the random surfer walking from i to j
- A: n x n square matrix
- weighted: whether to use the wedge weights (matrix entries)
- page_rank: array of length n. If *page_rank was null on entry, will be assigned.
-
- */
- int n = A->n;
- int i, j;
- int *ia = A->ia, *ja = A->ja;
- real *x, *y, *diag, res;
- real *a = NULL;
- int iter = 0;
-
- assert(A->m == n);
- assert(teleport_probablity >= 0);
-
- if (weighted){
- switch (A->type){
- case MATRIX_TYPE_REAL:
- a = (real*) A->a;
- break;
- case MATRIX_TYPE_COMPLEX:/* take real part */
- a = MALLOC(sizeof(real)*n);
- for (i = 0; i < n; i++) a[i] = ((real*) A->a)[2*i];
- break;
- case MATRIX_TYPE_INTEGER:
- a = MALLOC(sizeof(real)*n);
- for (i = 0; i < n; i++) a[i] = ((int*) A->a)[i];
- break;
- case MATRIX_TYPE_PATTERN:
- case MATRIX_TYPE_UNKNOWN:
- default:
- weighted = FALSE;
- break;
- }
- }
-
-
- if (!(*page_rank)) *page_rank = MALLOC(sizeof(real)*n);
- x = *page_rank;
-
- diag = MALLOC(sizeof(real)*n);
- for (i = 0; i < n; i++) diag[i] = 0;
- y = MALLOC(sizeof(real)*n);
-
- for (i = 0; i < n; i++) x[i] = 1./n;
-
- /* find the column sum */
- if (weighted){
- for (i = 0; i < n; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (ja[j] == i) continue;
- diag[i] += fabs(a[j]);
- }
- }
- } else {
- for (i = 0; i < n; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (ja[j] == i) continue;
- diag[i]++;
- }
- }
- }
- for (i = 0; i < n; i++) diag[i] = 1./MAX(diag[i], MACHINEACC);
-
- /* iterate */
- do {
- iter++;
- for (i = 0; i < n; i++) y[i] = 0;
- if (weighted){
- for (i = 0; i < n; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (ja[j] == i) continue;
- y[ja[j]] += a[j]*x[i]*diag[i];
- }
- }
- } else {
- for (i = 0; i < n; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (ja[j] == i) continue;
- y[ja[j]] += x[i]*diag[i];
- }
- }
- }
- for (i = 0; i < n; i++){
- y[i] = (1-teleport_probablity)*y[i] + teleport_probablity/n;
- }
-
- /*
- fprintf(stderr,"\n============\nx=");
- for (i = 0; i < n; i++) fprintf(stderr,"%f,",x[i]);
- fprintf(stderr,"\nx=");
- for (i = 0; i < n; i++) fprintf(stderr,"%f,",y[i]);
- fprintf(stderr,"\n");
- */
-
- res = 0;
- for (i = 0; i < n; i++) res += fabs(x[i] - y[i]);
- if (Verbose) fprintf(stderr,"page rank iter -- %d, res = %f\n",iter, res);
- memcpy(x, y, sizeof(real)*n);
- } while (res > epsilon);
-
- FREE(y);
- FREE(diag);
- if (a && a != A->a) FREE(a);
-}