how to efficiently apply function to all pairs of a set (and collect the results)

There are two built-in functions to generate pairs, either with (Tuples) or without (Subsets) duplication. Since your question states the number of iterations as $n*(n-1)/2$ I believe you want the latter:

set = {1, 2, 3, 4};
Subsets[set, {2}]

{{1, 2}, {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4}}

The short notation for Apply at level 1 is @@@, so this gives f for each pair:

f @@@ Subsets[set, {2}]

{f[1, 2], f[1, 3], f[1, 4], f[2, 3], f[2, 4], f[3, 4]}

This is in my opinion the most elegant code to produce this specific result, and it is quite fast, but it is not memory efficient if you only need to collect a result for a low percentage of the pairs. Let's define f as follows (I use {} here generically):

f[a_, b_] /; b~Divisible~a := {a, b}

f[___] = Sequence[];

If we now compute all the pairs for Range[5000] it takes one gigabyte of memory:

pairs = Range@5000 ~Subsets~ {2};
pairs // ByteCount

1099780032

And applying f we see that of the nearly 12.5 million pairs we need a return for only 38376 of them which takes only 3MB of storage:

r1 = f @@@ pairs;

ByteCount[r1]

Length[r1]

3377120

38376

Yet, the maximum memory used is 1.6GB:

 MaxMemoryUsed[]

1629793536

A simple method to reduce memory consumption is to process the subsets in blocks, rather than all at once, as follows:

set = Range@5000;
n = Length@set;
max = n (n - 1)/2;
block = 10000;

Timing[
 r2 =
   Join @@ Array[
     f @@@ Subsets[set, {2}, block {# - 1, #} + {1, 0}] &,
     ⌈max/block⌉
   ];
]

Length[r2]

MaxMemoryUsed[]

{8.299, Null}

38376

19769800

This only uses a maximum of ~20MB of memory, only a few MB over the baseline on my system.
(It issues a Subsets::take message but there is no error.)

My preferred method

Another method, and the one that I prefer, is to compute the pairs more manually allowing f to be embedded in the process so as to not generate all pairs beforehand. This uses Outer to effect the pair generation for each element separately (that is, all pairs starting with a certain element).

pairMap[f_, s_] := Module[{ss = s}, 
  Flatten[Outer[f, {#}, ss = Rest@ss, 1] & /@ Most@s, 2] ]

pairMap[f, Range@5000] // Length // Timing

MaxMemoryUsed[]

{7.816, 38376}

19430512

This also uses only a small amount of memory, and testing should bear out that it is faster as well. A variation of this method that may be even faster is to not build an output expression at all, relying instead on Sow and Reap to gather your results:

g[a_, b_] /; b ~Divisible~ a := Sow[ {a, b} ]

pairScan[f_, s_] := Module[{ss = s}, Outer[f, {#}, ss = Rest@ss, 1] & ~Scan~ Most@s ]

Reap[ pairScan[g, Range@5000] ][[2, 1]] // Length // Timing

MaxMemoryUsed[]

{6.583, 38376}

18757552

An argument for Do loops

While Outer is somewhat faster, after further consideration and conferring with jVincent I think perhaps after all a Do loop is as good as anything. One could write pairScan in this way:

pairScan2[f_, s_] := Module[{ss = s}, Do[f[i, j], {i, ss}, {j, ss = Rest@ss}] ]

Reap[ pairScan2[g, Range@5000] ][[2, 1]] // Length // Timing

MaxMemoryUsed[]

{7.613, 38376}

18711080

vertices = Range[10];
pairs = Tuples[vertices, 2];
func[x__] := First@x <= Last@x;
edges = Pick[pairs, func[#] & /@ pairs];
(* or *)
edges = Pick[pairs, Boole[func[#]] & /@ pairs, 1];
Graph[DirectedEdge @@@ edges, VertexLabels -> "Name", ImagePadding -> 20]

 Graph[DirectedEdge @@@ edges, VertexLabels -> "Name",  ImagePadding -> 20]

vertices = Range[100];
pairs = Tuples[vertices, 2];
func2[x__] := PrimeQ[(First@x)^(Last@x) + 3];
edges = Pick[pairs, func2[#] & /@ pairs];
Graph[DirectedEdge @@@ edges, ImagePadding -> 20,
  VertexLabels -> "Name",
  VertexLabelStyle -> Red,
  VertexStyle -> White,
  VertexShapeFunction -> "Point",
  EdgeStyle -> Directive[Opacity[.5], Hue[.15, .5, .8]],
  Background -> Black,
  EdgeShapeFunction -> ({Arrowheads[.01], Arrow[#1, 0.01]} &), 
  ImageSize -> 800]

Update: Using SparseArray with AdjacencyGraph and Graph:

func2[x__] := PrimeQ[(First@x)^(Last@x) + 3]; (* or your choice of test function *)
adjmtrx = SparseArray[{{i_, j_} /; func2[{i, j}] -> 1}, {20, 20}];
styling = {ImagePadding -> 20,
     VertexLabels -> "Name",
     VertexLabelStyle -> Red,
     VertexStyle -> White,
     VertexShapeFunction -> "Point",
     EdgeStyle -> Directive[Opacity[.5], Hue[.15, .5, .8]],
     Background -> Black,
     EdgeShapeFunction -> ({Arrowheads[.02], Arrow[#1, 0.01]} &), 
     ImageSize -> 400};

Row[{AdjacencyGraph[adjmtrx, styling], Spacer[5], 
      Graph[Range@First@Dimensions[adjmtrx], 
           DirectedEdge @@@ adjmtrx["NonzeroPositions"], styling]}]

Yes, there are definitely shorter, non-loop ways to do this. Define your link-defining function like this (with whatever test you need as the first argument to If:

islinked[a_Integer?Positive, b_Integer?Positive] := 
 If[Mod[a, 3] == 0 && Mod[b, 2] == 1, a -> b, {}]

You can then Apply this at the level of each row using the @@@ shorthand.
I would use Tuples rather than Subsets to allow for both 1->2 and 2->1. Eliminate the empty cases with a replacement rule as shown.

gg = islinked @@@ Tuples[Range[10], 2] /. {} -> Sequence[]

(* {3 -> 1, 3 -> 3, 3 -> 5, 3 -> 7, 3 -> 9, 6 -> 1, 6 -> 3, 6 -> 5, 
    6 -> 7, 6 -> 9, 9 -> 1, 9 -> 3, 9 -> 5, 9 -> 7, 9 -> 9} *)

Graph[gg]