**Built-in option**

This sidesteps most of your code, so it might not be what you are looking for, but I believe your goal can be achieved with _Mathematica_'s built-in image processing capability, specifically: `MorphologicalComponents`! 

Define a new `clustering` function

    clustering1[config_] := Module[{output, csizes, cindices},
       output =  MorphologicalComponents[Image@Abs@config, CornerNeighbors -> False];
       csizes = Rest@Sort@Tally@Flatten@output;
       cindices = Module[
          {sa = SparseArray[output], xx, yy, sa1, sa2},
          sa1 = sa["NonzeroValues"];
          xx = GatherBy[Range@Length@sa1, sa1[[#]] &];
          sa2 = sa["NonzeroPositions"];
          yy = sa2[[#]] & /@ xx;
          Transpose[{sa1[[xx[[All, 1]]]], yy}]];
         ];
       {csizes, cindices, output}
      ];

and apply:

    inputConfig = {{-1, -1, -1, 0}, {-1, 0, 0, -1}, {-1, -1, 0, -1}, {0, -1, 0, 0}};
    clustering1@inputConfig

>     {{{1, 7}, {2, 2}}, 
     {{1, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}}, {2, {{2, 4}, {3, 4}}}},
     {{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}

The first item in the output is the list of `{clusterNumber, clusterSize}`, the second the list of `{clusterNumber, clusterIndices}`, and the third is the input array with cluster numbers replacing filled sites.

**EDIT**

It appears that finding the cluster indices with `Position` is extremely slow with large datasets, although the rest of the code is relatively fast. With thanks to [ciao][1] (see below), I've replaced it with a much faster construction.

**Mimicking the recursive algorithm**

Here is my take on what the OP was actually asking for: a _Mathematic_ approach to a recursive algorithm. The actual algorithm is not much changed - the use of `Sow` and `Reap` in the OP's code is a good idea, and I agree that the recursion needed is not amenable to `Nest` or `Fold` - all I can do is clean it up a little and exchange the `Do` for a `Scan` over the populated sites:

    clustering2[config_] := Module[
       {output = config,
        cnum = 0, length = Length@config,
        csearch, clusters, c},
       
       csearch[{i_, j_}, cnum_] := If[
         output[[i, j]] == -1,
         output[[i, j]] = cnum;
         Sow[{i, j}];
         csearch[#, cnum] & /@ Select[
           {i, j} + # & /@ {{0, 1}, {0, -1}, {1, 0}, {-1, 0}},
           (1 <= First@# <= length && 1 <= Last@# <= length) &]
        ];
       
       clusters = Reap[
          Scan[
           If[output[[Sequence @@ #]] == -1,
              cnum++;
              c = Reap[csearch[#, cnum]][[2, 1]];
              Sow[{cnum, Length@c, c}]] &,
           SparseArray[config]["NonzeroPositions"]]
        ][[2, 1]];
       
       {clusters, output}
      ];

Applying this to `inputConfig` gives

>     {{{1, 7, {{1, 1}, {1, 2}, {1, 3}, {2, 1}, {3, 1}, {3, 2}, {4, 2}}}, 
      {2, 2, {{2, 4}, {3, 4}}}}, 
     {{1, 1, 1, 0}, {1, 0, 0, 2}, {1, 1, 0, 2}, {0, 1, 0, 0}}}

The first item in the output is the list of `{clusterNumber, clusterSize, clusterIndices}` and the second is the input array with cluster numbers replacing filled sites. 

I've run this on a 1000x1000 array with filling factor ~0.5, and although it is about 3 times slower than the first method, it does work. 

**Comments**

I'd be interested to know if anybody can think of a more clever way of implementing the algorithm. _Mathematica_ has the `SparseArray` object, which is a good fit for the input and output configurations we are dealing with here, so one might consider using those from the beginning for the input. Thanks to [ciao][1], who pointed out that one can extract the positions of non-zero elements easily from a `SparseArray` with `["NonzeroElements"]`, we can restrict our attention to those only.


  [1]: http://mathematica.stackexchange.com/users/11467/ciao