How do I generate the upper triangular indices from a list?

The solution is straightforward: Subsets, specifically

Subsets[{1,2,3}, {2}]

gives

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

To generate the lower indices, just Reverse them

Reverse /@ Subsets[{1,2,3}, {2}]

which gives

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

Just for the sake of functional programming here is a solution which is neither in speed nor memory comparable to Subsets but which only uses basic list operations and pattern matching:

makePairs[{e_, es__}] := Join[{e, #} & /@ {es}, makePairs[{es}]];
makePairs[{e_}] := {};

It grabs the first element of the list and makes pairs with the other elements. Then it calls itself with the list where the first element is removed and joins all results when the recursion finishes.

In[37]:= makePairs[{1, 2, 3}]

(*
  Out[37]= {{1, 2}, {1, 3}, {2, 3}}
*)

Your bonus question can be implemented by simply reversing {e,#} and changing the order of the arguments to Join.

Combinations in combinatorics have a special meaning. Accordingly, you can calculate the combinations for more dimensions. Note that these always contain the indices of diagonal elements. See threads on StackOverflow here (combinations) and here (combinations with repetition).

This solution is the function version of the answer provided by Yaroslav Bulatov (link above):

Combinations[elem_List] := Combinations[elem, All];
Combinations[elem_List, All | Full] := Combinations[elem, Length@Union@elem];
Combinations[elem_List, n_Integer] := Module[{coef2vars, coefs, set},
   set = Union@elem;
   coef2vars[list_] := Join @@ (MapIndexed[Table[set[[First@#2]], {#1}] &, list]);
   coefs = Flatten[Permutations /@ IntegerPartitions[n, {Length@set}, Range[0, n]], 1];
   Sort@(coef2vars /@ coefs)
   ];

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

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