nvgraphPagerankでPagerankを求める

CodeZineにPagerankアルゴリズムの解説記事を書きました。
推移確率行列:H に対し

　λI = HI

を満たす固有ベクトル(eigenvector):I がすなわち Pagerank となります。
※ Iは確率ベクトルなので総和は常に1、なので固有値(eigenvalue):λ は1.0 です。

この問題を解いてくれるのが nvGRAPH の nvgraphPagerank、CodeZineで使ったサンプルデータで nvgraphPagerank を使ってみましたよ。
※ codeの半分は COO-format から CSC-format への変換でできてます。

#include <cuda_runtime.h>
#include <nvgraph.h>

#include <iostream>
#include <algorithm>

using namespace std;

void pagerank(nvgraphHandle_t handle, nvgraphCSCTopology32I_t topology, float* weights, float*bookmark, float* probs, int iter);

template<typename T>
void dump(T* data, int n) {
  T* items = new T[n];
  cudaMemcpy(items, data, n*sizeof(T), cudaMemcpyDefault);
  for_each( items, items+n, [](T item) { cerr << item << " "; });
  cerr << endl;
  delete[] items;
}

int main() {
  const float r0 = 0.0f;
  const float r1 = 1.0f;
  const float r2 = 1.0f / 2.0f;
  const float r3 = 1.0f / 3.0f;
  const float r8 = 1.0f / 8.0f;

  // 元ネタ: COO-format
  const int nvertices =  8; // ページ数
  const int nedges    = 17; // リンク数
  int source_indices[nedges] = // リンク元
   {  0,  0,  1,  2,  2,  3,  3,  3,  4,  4,  4,  5,  6,  6,  6,  7,  7 }; 
  int destination_indices[nedges] = // リンク先
   {  1,  2,  3,  1,  4,  1,  4,  5,  5,  6,  7,  7,  0,  4,  7,  5,  6 };
  float weights[nedges] = // 各リンクの重み
   { r2, r2, r1, r2, r2, r3, r3, r3, r3, r3, r3, r1, r3, r3, r3, r2, r2 };

  // リンク先が少なくともひとつあるページは 1.0f、
  // リンク先のない行き止まりページは 0.0f
  float bookmark[nedges];
  fill_n(bookmark, nedges, 1.0f);
  for ( int dest : destination_indices ) {
    bookmark[dest] = 0.0f;
  }

  // 上記の元ネタを coo_topology に転写
  nvgraphCOOTopology32I_st coo_topology;
  float* coo_weights;
  cudaMalloc(&coo_topology.source_indices     , nedges*sizeof(int));
  cudaMalloc(&coo_topology.destination_indices, nedges*sizeof(int));
  cudaMalloc(&coo_weights, nedges*sizeof(float));

  coo_topology.nvertices = nvertices;
  coo_topology.nedges    = nedges;
  cudaMemcpy(coo_topology.source_indices     , source_indices     , nedges*sizeof(int), cudaMemcpyDefault);
  cudaMemcpy(coo_topology.destination_indices, destination_indices, nedges*sizeof(int), cudaMemcpyDefault);
  coo_topology.tag = NVGRAPH_UNSORTED;
  cudaMemcpy(coo_weights, weights, nedges*sizeof(float), cudaMemcpyDefault);

  // nvGRAPH の出番はここから
  nvgraphHandle_t handle;
  nvgraphCreate(&handle);
  cudaDataType_t data_type = CUDA_R_32F; // 32bit-real (=float)

  {
  // Pagerankを求めるには CSR-formatに変換せにゃならん
  nvgraphCSCTopology32I_st csc_topology;
  float* csc_weights;
  cudaMalloc(&csc_topology.source_indices, nedges*sizeof(int));
  cudaMalloc(&csc_topology.destination_offsets, (nvertices+1)*sizeof(int));
  cudaMalloc(&csc_weights, nedges*sizeof(float));

  nvgraphConvertTopology(handle, 
                         NVGRAPH_COO_32, &coo_topology, coo_weights,  // in  type, topology, weigts
                         &data_type,                                  // type of weights
                         NVGRAPH_CSC_32, &csc_topology, csc_weights); // out type, topology, weigts

  cerr << "CSC : " <<endl;
  cerr << "nvertices           = " << csc_topology.nvertices << endl;
  cerr << "nedges              = " << csc_topology.nedges    << endl;
  cerr << "source_indices      = "; dump(csc_topology.source_indices     , csc_topology.nedges);
  cerr << "destination_offsets = "; dump(csc_topology.destination_offsets, csc_topology.nvertices+1);
  cerr << "weights             = "; dump(csc_weights, csc_topology.nedges);

  // Pagerankを求めて出力
  float probs[nvertices] = { r1, r0, r0, r0, r0, r0, r0, r0 };
  pagerank(handle, &csc_topology, csc_weights, bookmark, probs, 100);
  cout << "\nPagerank: \n";
  for ( int i = 0; i < nvertices; ++i ) {
    cout << "page[" << i << "] : " << probs[i] << endl;
  }

  cudaFree(csc_topology.source_indices);
  cudaFree(csc_topology.destination_offsets);
  cudaFree(csc_weights);
  }

  cudaFree(coo_topology.source_indices);
  cudaFree(coo_topology.destination_indices);
  cudaFree(coo_weights);

  nvgraphDestroy(handle);
  cudaDeviceReset();
}

void pagerank(nvgraphHandle_t handle, nvgraphCSCTopology32I_t topology, float* weights, float* bookmark, float* probs, int iter) {
  nvgraphGraphDescr_t descr;
  // setup
  nvgraphCreateGraphDescr(handle, &descr);

  // topology
  nvgraphSetGraphStructure(handle, descr, topology, NVGRAPH_CSC_32);
  cudaDataType_t dimT[2] = { CUDA_R_32F, CUDA_R_32F };
  nvgraphAllocateVertexData(handle, descr, 2, dimT);
  nvgraphAllocateEdgeData(handle, descr, 1, dimT);
  // weight
  nvgraphSetEdgeData(handle, descr, weights, 0);
  // bookmark
  nvgraphSetVertexData(handle, descr, bookmark, 0);
  // initial guess
  nvgraphSetVertexData(handle, descr, probs, 1);

  // pagerank!!
  float damping_factor = 0.85f; // 85%はリンクを辿り,残る15%は気ままにジャンプ
  nvgraphPagerank(handle, descr,
    0,       // weight    : EdgeData[0] 
    &damping_factor, 
    0,       // bookmark  : VertexData[0]
    1,       // has guess
    1,       // pagerank  : VertexData[1]
    1.0e-6f, // tolerance 誤差 : これより小さくなったら計算終了
    iter     // 最大繰り返し回数
  );
  // get result
  nvgraphGetVertexData(handle, descr, probs, 1);

  // teardown
  nvgraphDestroyGraphDescr(handle, descr);
}