haskell-igraph-0.8.0: igraph/src/prpack_solver.cpp
#include "prpack_solver.h"
#include "prpack_utils.h"
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <algorithm>
using namespace prpack;
using namespace std;
void prpack_solver::initialize() {
geg = NULL;
gsg = NULL;
sg = NULL;
sccg = NULL;
owns_bg = true;
}
prpack_solver::prpack_solver(const prpack_csc* g) {
initialize();
TIME(read_time, bg = new prpack_base_graph(g));
}
prpack_solver::prpack_solver(const prpack_int64_csc* g) {
initialize();
TIME(read_time, bg = new prpack_base_graph(g));
}
prpack_solver::prpack_solver(const prpack_csr* g) {
initialize();
TIME(read_time, bg = new prpack_base_graph(g));
}
prpack_solver::prpack_solver(const prpack_edge_list* g) {
initialize();
TIME(read_time, bg = new prpack_base_graph(g));
}
prpack_solver::prpack_solver(prpack_base_graph* g, bool owns_bg) {
initialize();
this->owns_bg = owns_bg;
TIME(read_time, bg = g);
}
prpack_solver::prpack_solver(const char* filename, const char* format, const bool weighted) {
initialize();
TIME(read_time, bg = new prpack_base_graph(filename, format, weighted));
}
prpack_solver::~prpack_solver() {
if (owns_bg) {
delete bg;
}
delete geg;
delete gsg;
delete sg;
delete sccg;
}
int prpack_solver::get_num_vs() {
return bg->num_vs;
}
prpack_result* prpack_solver::solve(const double alpha, const double tol, const char* method) {
return solve(alpha, tol, NULL, NULL, method);
}
prpack_result* prpack_solver::solve(
const double alpha,
const double tol,
const double* u,
const double* v,
const char* method) {
double preprocess_time = 0;
double compute_time = 0;
prpack_result* ret = NULL;
// decide which method to run
string m;
if (strcmp(method, "") != 0)
m = string(method);
else {
if (bg->num_vs < 128)
m = "ge";
else if (sccg != NULL)
m = "sccgs";
else if (sg != NULL)
m = "sg";
else
m = "sccgs";
if (u != v)
m += "_uv";
}
// run the appropriate method
if (m == "ge") {
if (geg == NULL) {
TIME(preprocess_time, geg = new prpack_preprocessed_ge_graph(bg));
}
TIME(compute_time, ret = solve_via_ge(
alpha,
tol,
geg->num_vs,
geg->matrix,
u));
} else if (m == "ge_uv") {
if (geg == NULL) {
TIME(preprocess_time, geg = new prpack_preprocessed_ge_graph(bg));
}
TIME(compute_time, ret = solve_via_ge_uv(
alpha,
tol,
geg->num_vs,
geg->matrix,
geg->d,
u,
v));
} else if (m == "gs") {
if (gsg == NULL) {
TIME(preprocess_time, gsg = new prpack_preprocessed_gs_graph(bg));
}
TIME(compute_time, ret = solve_via_gs(
alpha,
tol,
gsg->num_vs,
gsg->num_es,
gsg->heads,
gsg->tails,
gsg->vals,
gsg->ii,
gsg->d,
gsg->num_outlinks,
u,
v));
} else if (m == "gserr") {
if (gsg == NULL) {
TIME(preprocess_time, gsg = new prpack_preprocessed_gs_graph(bg));
}
TIME(compute_time, ret = solve_via_gs_err(
alpha,
tol,
gsg->num_vs,
gsg->num_es,
gsg->heads,
gsg->tails,
gsg->ii,
gsg->num_outlinks,
u,
v));
} else if (m == "sgs") {
if (sg == NULL) {
TIME(preprocess_time, sg = new prpack_preprocessed_schur_graph(bg));
}
TIME(compute_time, ret = solve_via_schur_gs(
alpha,
tol,
sg->num_vs,
sg->num_no_in_vs,
sg->num_no_out_vs,
sg->num_es,
sg->heads,
sg->tails,
sg->vals,
sg->ii,
sg->d,
sg->num_outlinks,
u,
sg->encoding,
sg->decoding));
} else if (m == "sgs_uv") {
if (sg == NULL) {
TIME(preprocess_time, sg = new prpack_preprocessed_schur_graph(bg));
}
TIME(compute_time, ret = solve_via_schur_gs_uv(
alpha,
tol,
sg->num_vs,
sg->num_no_in_vs,
sg->num_no_out_vs,
sg->num_es,
sg->heads,
sg->tails,
sg->vals,
sg->ii,
sg->d,
sg->num_outlinks,
u,
v,
sg->encoding,
sg->decoding));
} else if (m == "sccgs") {
if (sccg == NULL) {
TIME(preprocess_time, sccg = new prpack_preprocessed_scc_graph(bg));
}
TIME(compute_time, ret = solve_via_scc_gs(
alpha,
tol,
sccg->num_vs,
sccg->num_es_inside,
sccg->heads_inside,
sccg->tails_inside,
sccg->vals_inside,
sccg->num_es_outside,
sccg->heads_outside,
sccg->tails_outside,
sccg->vals_outside,
sccg->ii,
sccg->d,
sccg->num_outlinks,
u,
sccg->num_comps,
sccg->divisions,
sccg->encoding,
sccg->decoding));
} else if (m == "sccgs_uv") {
if (sccg == NULL) {
TIME(preprocess_time, sccg = new prpack_preprocessed_scc_graph(bg));
}
TIME(compute_time, ret = solve_via_scc_gs_uv(
alpha,
tol,
sccg->num_vs,
sccg->num_es_inside,
sccg->heads_inside,
sccg->tails_inside,
sccg->vals_inside,
sccg->num_es_outside,
sccg->heads_outside,
sccg->tails_outside,
sccg->vals_outside,
sccg->ii,
sccg->d,
sccg->num_outlinks,
u,
v,
sccg->num_comps,
sccg->divisions,
sccg->encoding,
sccg->decoding));
} else {
// TODO: throw exception
}
ret->method = m.c_str();
ret->read_time = read_time;
ret->preprocess_time = preprocess_time;
ret->compute_time = compute_time;
ret->num_vs = bg->num_vs;
ret->num_es = bg->num_es;
return ret;
}
// VARIOUS SOLVING METHODS ////////////////////////////////////////////////////////////////////////
prpack_result* prpack_solver::solve_via_ge(
const double alpha,
const double tol,
const int num_vs,
const double* matrix,
const double* uv) {
prpack_result* ret = new prpack_result();
// initialize uv values
const double uv_const = 1.0/num_vs;
const int uv_exists = (uv) ? 1 : 0;
uv = (uv) ? uv : &uv_const;
// create matrix A
double* A = new double[num_vs*num_vs];
for (int i = 0; i < num_vs*num_vs; ++i)
A[i] = -alpha*matrix[i];
for (int i = 0; i < num_vs*num_vs; i += num_vs + 1)
++A[i];
// create vector b
double* b = new double[num_vs];
for (int i = 0; i < num_vs; ++i)
b[i] = uv[uv_exists*i];
// solve and normalize
ge(num_vs, A, b);
normalize(num_vs, b);
// clean up and return
delete[] A;
ret->num_es_touched = -1;
ret->x = b;
return ret;
}
prpack_result* prpack_solver::solve_via_ge_uv(
const double alpha,
const double tol,
const int num_vs,
const double* matrix,
const double* d,
const double* u,
const double* v) {
prpack_result* ret = new prpack_result();
// initialize u and v values
const double u_const = 1.0/num_vs;
const double v_const = 1.0/num_vs;
const int u_exists = (u) ? 1 : 0;
const int v_exists = (v) ? 1 : 0;
u = (u) ? u : &u_const;
v = (v) ? v : &v_const;
// create matrix A
double* A = new double[num_vs*num_vs];
for (int i = 0; i < num_vs*num_vs; ++i)
A[i] = -alpha*matrix[i];
for (int i = 0, inum_vs = 0; i < num_vs; ++i, inum_vs += num_vs)
for (int j = 0; j < num_vs; ++j)
A[inum_vs + j] -= alpha*u[u_exists*i]*d[j];
for (int i = 0; i < num_vs*num_vs; i += num_vs + 1)
++A[i];
// create vector b
double* b = new double[num_vs];
for (int i = 0; i < num_vs; ++i)
b[i] = (1 - alpha)*v[v_exists*i];
// solve
ge(num_vs, A, b);
// clean up and return
delete[] A;
ret->num_es_touched = -1;
ret->x = b;
return ret;
}
// Vanilla Gauss-Seidel.
prpack_result* prpack_solver::solve_via_gs(
const double alpha,
const double tol,
const int num_vs,
const int num_es,
const int* heads,
const int* tails,
const double* vals,
const double* ii,
const double* d,
const double* num_outlinks,
const double* u,
const double* v) {
prpack_result* ret = new prpack_result();
const bool weighted = vals != NULL;
// initialize u and v values
const double u_const = 1.0/num_vs;
const double v_const = 1.0/num_vs;
const int u_exists = (u) ? 1 : 0;
const int v_exists = (v) ? 1 : 0;
u = (u) ? u : &u_const;
v = (v) ? v : &v_const;
// initialize the eigenvector (and use personalization vector)
double* x = new double[num_vs];
for (int i = 0; i < num_vs; ++i)
x[i] = 0;
// initialize delta
double delta = 0;
// run Gauss-Seidel
ret->num_es_touched = 0;
double err = 1, c = 0;
do {
if (weighted) {
for (int i = 0; i < num_vs; ++i) {
double new_val = 0;
const int start_j = tails[i];
const int end_j = (i + 1 != num_vs) ? tails[i + 1] : num_es;
for (int j = start_j; j < end_j; ++j)
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads[j]]*vals[j];
new_val = alpha*new_val + (1 - alpha)*v[v_exists*i];
delta -= alpha*x[i]*d[i];
new_val += delta*u[u_exists*i];
new_val /= 1 - alpha*(d[i]*u[u_exists*i] + (1 - d[i])*ii[i]);
delta += alpha*new_val*d[i];
COMPENSATED_SUM(err, x[i] - new_val, c);
x[i] = new_val;
}
} else {
for (int i = 0; i < num_vs; ++i) {
const double old_val = x[i]*num_outlinks[i];
double new_val = 0;
const int start_j = tails[i];
const int end_j = (i + 1 != num_vs) ? tails[i + 1] : num_es;
for (int j = start_j; j < end_j; ++j)
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads[j]];
new_val = alpha*new_val + (1 - alpha)*v[v_exists*i];
if (num_outlinks[i] < 0) {
delta -= alpha*old_val;
new_val += delta*u[u_exists*i];
new_val /= 1 - alpha*u[u_exists*i];
delta += alpha*new_val;
} else {
new_val += delta*u[u_exists*i];
new_val /= 1 - alpha*ii[i];
}
COMPENSATED_SUM(err, old_val - new_val, c);
x[i] = new_val/num_outlinks[i];
}
}
// update iteration index
ret->num_es_touched += num_es;
} while (err >= tol);
// undo num_outlinks transformation
if (!weighted)
for (int i = 0; i < num_vs; ++i)
x[i] *= num_outlinks[i];
// return results
ret->x = x;
return ret;
}
// Implement a gauss-seidel-like process with a strict error bound
// we return a solution with 1-norm error less than tol.
prpack_result* prpack_solver::solve_via_gs_err(
const double alpha,
const double tol,
const int num_vs,
const int num_es,
const int* heads,
const int* tails,
const double* ii,
const double* num_outlinks,
const double* u,
const double* v) {
prpack_result* ret = new prpack_result();
// initialize u and v values
const double u_const = 1.0/num_vs;
const double v_const = 1.0/num_vs;
const int u_exists = (u) ? 1 : 0;
const int v_exists = (v) ? 1 : 0;
u = (u) ? u : &u_const;
v = (v) ? v : &v_const;
// Note to Dave, we can't rescale v because we could be running this
// same routine from multiple threads.
// initialize the eigenvector (and use personalization vector)
double* x = new double[num_vs];
for (int i = 0; i < num_vs; ++i) {
x[i] = 0.;
}
// initialize delta
double delta = 0.;
// run Gauss-Seidel, note that we store x/deg[i] throughout this
// iteration.
int64_t maxedges = (int64_t)((double)num_es*std::min(
log(tol)/log(alpha),
(double)PRPACK_SOLVER_MAX_ITERS));
ret->num_es_touched = 0;
double err=1., c = 0.;
do {
// iterate through vertices
for (int i = 0; i < num_vs; ++i) {
double old_val = x[i]*num_outlinks[i]; // adjust back to the "true" value.
double new_val = 0.;
int start_j = tails[i], end_j = (i + 1 != num_vs) ? tails[i + 1] : num_es;
for (int j = start_j; j < end_j; ++j) {
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads[j]];
}
new_val = alpha*new_val + alpha*ii[i]*old_val + (1.0-alpha)*v[v_exists*i];
new_val += delta*u[u_exists*i]; // add the dangling node adjustment
if (num_outlinks[i] < 0) {
delta += alpha*(new_val - old_val);
}
// note that new_val > old_val, but the fabs is just for
COMPENSATED_SUM(err, -(new_val - old_val), c);
x[i] = new_val/num_outlinks[i];
}
// update iteration index
ret->num_es_touched += num_es;
} while (err >= tol && ret->num_es_touched < maxedges);
if (err >= tol) {
ret->converged = 0;
} else {
ret->converged = 1;
}
// undo num_outlinks transformation
for (int i = 0; i < num_vs; ++i)
x[i] *= num_outlinks[i];
// return results
ret->x = x;
return ret;
}
// Gauss-Seidel using the Schur complement to separate dangling nodes.
prpack_result* prpack_solver::solve_via_schur_gs(
const double alpha,
const double tol,
const int num_vs,
const int num_no_in_vs,
const int num_no_out_vs,
const int num_es,
const int* heads,
const int* tails,
const double* vals,
const double* ii,
const double* d,
const double* num_outlinks,
const double* uv,
const int* encoding,
const int* decoding,
const bool should_normalize) {
prpack_result* ret = new prpack_result();
const bool weighted = vals != NULL;
// initialize uv values
const double uv_const = 1.0/num_vs;
const int uv_exists = (uv) ? 1 : 0;
uv = (uv) ? prpack_utils::permute(num_vs, uv, encoding) : &uv_const;
// initialize the eigenvector (and use personalization vector)
double* x = new double[num_vs];
for (int i = 0; i < num_vs - num_no_out_vs; ++i)
x[i] = uv[uv_exists*i]/(1 - alpha*ii[i])/((weighted) ? 1 : num_outlinks[i]);
// run Gauss-Seidel for the top left part of (I - alpha*P)*x = uv
ret->num_es_touched = 0;
double err, c;
do {
// iterate through vertices
int num_es_touched = 0;
err = c = 0;
#pragma omp parallel for firstprivate(c) reduction(+:err, num_es_touched) schedule(dynamic, 64)
for (int i = num_no_in_vs; i < num_vs - num_no_out_vs; ++i) {
double new_val = 0;
const int start_j = tails[i];
const int end_j = (i + 1 != num_vs) ? tails[i + 1] : num_es;
if (weighted) {
for (int j = start_j; j < end_j; ++j)
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads[j]]*vals[j];
COMPENSATED_SUM(err, fabs(uv[uv_exists*i] + alpha*new_val - (1 - alpha*ii[i])*x[i]), c);
new_val = (alpha*new_val + uv[uv_exists*i])/(1 - alpha*ii[i]);
x[i] = new_val;
} else {
for (int j = start_j; j < end_j; ++j)
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads[j]];
COMPENSATED_SUM(err, fabs(uv[uv_exists*i] + alpha*new_val - (1 - alpha*ii[i])*x[i]*num_outlinks[i]), c);
new_val = (alpha*new_val + uv[uv_exists*i])/(1 - alpha*ii[i]);
x[i] = new_val/num_outlinks[i];
}
num_es_touched += end_j - start_j;
}
// update iteration index
ret->num_es_touched += num_es_touched;
} while (err/(1 - alpha) >= tol);
// solve for the dangling nodes
int num_es_touched = 0;
#pragma omp parallel for reduction(+:num_es_touched) schedule(dynamic, 64)
for (int i = num_vs - num_no_out_vs; i < num_vs; ++i) {
x[i] = 0;
const int start_j = tails[i];
const int end_j = (i + 1 != num_vs) ? tails[i + 1] : num_es;
for (int j = start_j; j < end_j; ++j)
x[i] += x[heads[j]]*((weighted) ? vals[j] : 1);
x[i] = (alpha*x[i] + uv[uv_exists*i])/(1 - alpha*ii[i]);
num_es_touched += end_j - start_j;
}
ret->num_es_touched += num_es_touched;
// undo num_outlinks transformation
if (!weighted)
for (int i = 0; i < num_vs - num_no_out_vs; ++i)
x[i] *= num_outlinks[i];
// normalize x to get the solution for: (I - alpha*P - alpha*u*d')*x = (1 - alpha)*v
if (should_normalize)
normalize(num_vs, x);
// return results
ret->x = prpack_utils::permute(num_vs, x, decoding);
delete[] x;
if (uv_exists)
delete[] uv;
return ret;
}
prpack_result* prpack_solver::solve_via_schur_gs_uv(
const double alpha,
const double tol,
const int num_vs,
const int num_no_in_vs,
const int num_no_out_vs,
const int num_es,
const int* heads,
const int* tails,
const double* vals,
const double* ii,
const double* d,
const double* num_outlinks,
const double* u,
const double* v,
const int* encoding,
const int* decoding) {
// solve uv = u
prpack_result* ret_u = solve_via_schur_gs(
alpha,
tol,
num_vs,
num_no_in_vs,
num_no_out_vs,
num_es,
heads,
tails,
vals,
ii,
d,
num_outlinks,
u,
encoding,
decoding,
false);
// solve uv = v
prpack_result* ret_v = solve_via_schur_gs(
alpha,
tol,
num_vs,
num_no_in_vs,
num_no_out_vs,
num_es,
heads,
tails,
vals,
ii,
d,
num_outlinks,
v,
encoding,
decoding,
false);
// combine the u and v cases
return combine_uv(num_vs, d, num_outlinks, encoding, alpha, ret_u, ret_v);
}
/** Gauss-Seidel using strongly connected components.
* Notes:
* If not weighted, then we store x[i] = "x[i]/outdegree" to
* avoid additional arithmetic. We don't do this for the weighted
* case because the adjustment may not be constant.
*/
prpack_result* prpack_solver::solve_via_scc_gs(
const double alpha,
const double tol,
const int num_vs,
const int num_es_inside,
const int* heads_inside,
const int* tails_inside,
const double* vals_inside,
const int num_es_outside,
const int* heads_outside,
const int* tails_outside,
const double* vals_outside,
const double* ii,
const double* d,
const double* num_outlinks,
const double* uv,
const int num_comps,
const int* divisions,
const int* encoding,
const int* decoding,
const bool should_normalize) {
prpack_result* ret = new prpack_result();
const bool weighted = vals_inside != NULL;
// initialize uv values
const double uv_const = 1.0/num_vs;
const int uv_exists = (uv) ? 1 : 0;
uv = (uv) ? prpack_utils::permute(num_vs, uv, encoding) : &uv_const;
// CHECK initialize the solution with one iteration of GS from x=0.
double* x = new double[num_vs];
for (int i = 0; i < num_vs; ++i)
x[i] = uv[uv_exists*i]/(1 - alpha*ii[i])/((weighted) ? 1 : num_outlinks[i]);
// create x_outside
double* x_outside = new double[num_vs];
// run Gauss-Seidel for (I - alpha*P)*x = uv
ret->num_es_touched = 0;
for (int comp_i = 0; comp_i < num_comps; ++comp_i) {
const int start_comp = divisions[comp_i];
const int end_comp = (comp_i + 1 != num_comps) ? divisions[comp_i + 1] : num_vs;
const bool parallelize = end_comp - start_comp > 512;
// initialize relevant x_outside values
for (int i = start_comp; i < end_comp; ++i) {
x_outside[i] = 0;
const int start_j = tails_outside[i];
const int end_j = (i + 1 != num_vs) ? tails_outside[i + 1] : num_es_outside;
for (int j = start_j; j < end_j; ++j)
x_outside[i] += x[heads_outside[j]]*((weighted) ? vals_outside[j] : 1.);
ret->num_es_touched += end_j - start_j;
}
double err, c;
do {
int num_es_touched = 0;
err = c = 0;
if (parallelize) {
// iterate through vertices
#pragma omp parallel for firstprivate(c) reduction(+:err, num_es_touched) schedule(dynamic, 64)
for (int i = start_comp; i < end_comp; ++i) {
double new_val = x_outside[i];
const int start_j = tails_inside[i];
const int end_j = (i + 1 != num_vs) ? tails_inside[i + 1] : num_es_inside;
if (weighted) {
for (int j = start_j; j < end_j; ++j) {
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads_inside[j]]*vals_inside[j];
}
COMPENSATED_SUM(err, fabs(uv[uv_exists*i] + alpha*new_val - (1 - alpha*ii[i])*x[i]), c);
x[i] = (alpha*new_val + uv[uv_exists*i])/(1 - alpha*ii[i]);
} else {
for (int j = start_j; j < end_j; ++j) {
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads_inside[j]];
}
COMPENSATED_SUM(err, fabs(uv[uv_exists*i] + alpha*new_val - (1 - alpha*ii[i])*x[i]*num_outlinks[i]), c);
x[i] = (alpha*new_val + uv[uv_exists*i])/(1 - alpha*ii[i])/num_outlinks[i];
}
num_es_touched += end_j - start_j;
}
} else {
for (int i = start_comp; i < end_comp; ++i) {
double new_val = x_outside[i];
const int start_j = tails_inside[i];
const int end_j = (i + 1 != num_vs) ? tails_inside[i + 1] : num_es_inside;
if (weighted) {
for (int j = start_j; j < end_j; ++j) {
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads_inside[j]]*vals_inside[j];
}
COMPENSATED_SUM(err, fabs(uv[uv_exists*i] + alpha*new_val - (1 - alpha*ii[i])*x[i]), c);
x[i] = (alpha*new_val + uv[uv_exists*i])/(1 - alpha*ii[i]);
} else {
for (int j = start_j; j < end_j; ++j) {
// TODO: might want to use compensation summation for large: end_j - start_j
new_val += x[heads_inside[j]];
}
COMPENSATED_SUM(err, fabs(uv[uv_exists*i] + alpha*new_val - (1 - alpha*ii[i])*x[i]*num_outlinks[i]), c);
x[i] = (alpha*new_val + uv[uv_exists*i])/(1 - alpha*ii[i])/num_outlinks[i];
}
num_es_touched += end_j - start_j;
}
}
// update iteration index
ret->num_es_touched += num_es_touched;
} while (err/(1 - alpha) >= tol*(end_comp - start_comp)/num_vs);
}
// undo num_outlinks transformation
if (!weighted)
for (int i = 0; i < num_vs; ++i)
x[i] *= num_outlinks[i];
// normalize x to get the solution for: (I - alpha*P - alpha*u*d')*x = (1 - alpha)*v
if (should_normalize)
normalize(num_vs, x);
// return results
ret->x = prpack_utils::permute(num_vs, x, decoding);
delete[] x;
delete[] x_outside;
if (uv_exists)
delete[] uv;
return ret;
}
prpack_result* prpack_solver::solve_via_scc_gs_uv(
const double alpha,
const double tol,
const int num_vs,
const int num_es_inside,
const int* heads_inside,
const int* tails_inside,
const double* vals_inside,
const int num_es_outside,
const int* heads_outside,
const int* tails_outside,
const double* vals_outside,
const double* ii,
const double* d,
const double* num_outlinks,
const double* u,
const double* v,
const int num_comps,
const int* divisions,
const int* encoding,
const int* decoding) {
// solve uv = u
prpack_result* ret_u = solve_via_scc_gs(
alpha,
tol,
num_vs,
num_es_inside,
heads_inside,
tails_inside,
vals_inside,
num_es_outside,
heads_outside,
tails_outside,
vals_outside,
ii,
d,
num_outlinks,
u,
num_comps,
divisions,
encoding,
decoding,
false);
// solve uv = v
prpack_result* ret_v = solve_via_scc_gs(
alpha,
tol,
num_vs,
num_es_inside,
heads_inside,
tails_inside,
vals_inside,
num_es_outside,
heads_outside,
tails_outside,
vals_outside,
ii,
d,
num_outlinks,
v,
num_comps,
divisions,
encoding,
decoding,
false);
// combine u and v
return combine_uv(num_vs, d, num_outlinks, encoding, alpha, ret_u, ret_v);
}
// VARIOUS HELPER METHODS /////////////////////////////////////////////////////////////////////////
// Run Gaussian-Elimination (note: this changes A and returns the solution in b)
void prpack_solver::ge(const int sz, double* A, double* b) {
// put into triangular form
for (int i = 0, isz = 0; i < sz; ++i, isz += sz)
for (int k = 0, ksz = 0; k < i; ++k, ksz += sz)
if (A[isz + k] != 0) {
const double coeff = A[isz + k]/A[ksz + k];
A[isz + k] = 0;
for (int j = k + 1; j < sz; ++j)
A[isz + j] -= coeff*A[ksz + j];
b[i] -= coeff*b[k];
}
// backwards substitution
for (int i = sz - 1, isz = (sz - 1)*sz; i >= 0; --i, isz -= sz) {
for (int j = i + 1; j < sz; ++j)
b[i] -= A[isz + j]*b[j];
b[i] /= A[isz + i];
}
}
// Normalize a vector to sum to 1.
void prpack_solver::normalize(const int length, double* x) {
double norm = 0, c = 0;
for (int i = 0; i < length; ++i) {
COMPENSATED_SUM(norm, x[i], c);
}
norm = 1/norm;
for (int i = 0; i < length; ++i)
x[i] *= norm;
}
// Combine u and v results.
prpack_result* prpack_solver::combine_uv(
const int num_vs,
const double* d,
const double* num_outlinks,
const int* encoding,
const double alpha,
const prpack_result* ret_u,
const prpack_result* ret_v) {
prpack_result* ret = new prpack_result();
const bool weighted = d != NULL;
double delta_u = 0;
double delta_v = 0;
for (int i = 0; i < num_vs; ++i) {
if ((weighted) ? (d[encoding[i]] == 1) : (num_outlinks[encoding[i]] < 0)) {
delta_u += ret_u->x[i];
delta_v += ret_v->x[i];
}
}
const double s = ((1 - alpha)*alpha*delta_v)/(1 - alpha*delta_u);
const double t = 1 - alpha;
ret->x = new double[num_vs];
for (int i = 0; i < num_vs; ++i)
ret->x[i] = s*ret_u->x[i] + t*ret_v->x[i];
ret->num_es_touched = ret_u->num_es_touched + ret_v->num_es_touched;
// clean up and return
delete ret_u;
delete ret_v;
return ret;
}