Simple matrix, vector, permutation utilities

I am going to be rolling out the rest of my branch-and-bound algorithm in the next few posts.  To make that easier, in this post I introduce some common matrix, vector, and permutation methods.  It turns out that for technical computing applications, C#'s extension methods (introduced in 3.0) are awesome.  With vectors it's great because you retain the control of having direct array access, but you get the nice object-oriented notation.

I've left out comments and error-checking, and it's not super-optimized, but you get the idea.  None of these methods end up being the bottleneck.  Most of the methods are obvious, but I do want to comment on the permutation methods.  The goal of QAP is to find an optimal assignment of facilities to locations, represented as a permutation.  In the course of our branch-and-bound algorithm, we'll be making partial assignments - that is, only some of the entries in the permutation will be filled in.  By convention, p == -1 will mean that facility i is unassigned.  An index where p < 0 is called "unused".  When we branch, we'll want to pick an unused facility (or location), and try assigning all unused locations (or facilities) to it.   The use of C#'s iterator and yield concepts really help here.

That said, here's the code.

public static class MatrixUtilities {
  #region Vector
  public static void ConstantFill(this T[] data, T val)
  {
    for (int i = 0; i < data.Length; i++) {
      data = val;
    }
  }

  public static int ArgMin(this double[] data, out double best) {
    if (data.Length == 0) {
      best = Double.MinValue;
      return -1;
    }
    int iBest = 0;
    best = data;
    for (int i = 1; i < data.Length; i++) {
      if (best > data) {
        best = data;
        iBest = i;
      }
    }
    return iBest;
  }

  public static int ArgMax(this double[] data, out double best) {
    if (data.Length == 0) { best = Double.MinValue; return -1; }
    int iBest = 0;
    best = data;
    for (int i = 1; i < data.Length; i++) {
      if (best < data) {
        best = data;
        iBest = i;
      }
    }
    return iBest;
  }
#endregion

#region Permutation
public static void Swap(this int[] p, int i, int j) {
  int temp = p;
  p = p;
  p = temp;
}

public static int FindUnused(this int[] data, int index) {
  foreach (int unused in data.UnusedIndices()) {
    if (index-- <= 0) {
      return unused;
    }
  }
  return -1;
}

public static IEnumerable UnusedIndices(this int[] data) {
  for (int i = 0; i < data.Length; i++) {
    if (data < 0) {
      yield return i;
    }
  }
}

public static string PermutationToString(this int[] p, bool oneBased) {
  StringBuilder build = new StringBuilder(p.Length * 4);
  int width = ((int)Math.Log10(p.Length)) + 2;
  build.Append("[");
  for (int i = 0; i < p.Length; i++) {
    int index = oneBased ? p + 1 : p;
    build.Append(index.ToString().PadLeft(width));
  }
  build.Append("]");
  return build.ToString();
}

#endregion

#region Matrix
public static string MatrixToString(this double[][] A) {
  if (A != null) {
    StringBuilder build = new StringBuilder(A.Length * A.Length * 3);
    for (int i = 0; i < A.Length; i++) {
      if (A != null) {
        for (int j = 0; j < A.Length; j++) {
          build.AppendFormat("{0,4}", A);
          build.Append(" ");
        }
        build.AppendLine();
      }
    }
    return build.ToString();
  }
  return null;
}

public static double[][] NewMatrix(int m, int n) {
  double[][] M = new double[];
  for (int i = 0; i < M.Length; i++) {
    M = new double;
  }
  return M;
}

public static double[][] NewMatrix(int n) {
  return NewMatrix(n, n);
}
#endregion
}