Can I get the geodist matrix of all pairs(N*N) directly? Thanks!

[whitesay]

Instead of using for loop to get the geodist, i.e. : solver.compute_distance(i), return N \times 1.
Can directly return the matrix of N \times N.

[nmwsharp]

There is no built-in for this currently.

Interestingly, most algorithms for computing geodesic distance are not any faster at computing all-pairs distance than just computing each distance vector yourself and stacking them. So simply calling the existing routines in a for-loop with an np.concatenate() at the end is a pretty good solution for this.