Most elegant way to pass and return 2D Array (Matrix) typesafe with special constraints (NxN) in C++

Question 1

Your template idea is surely doable. The interesting part is the squareMatrixT class declaration with the typedef. The defined type is a reference to a one-dimensional array of the specified size. That's what the index operator returns (a plain old array, well, as ref). The rest of the code is just proof of concept clutter (objects of this type can be passed as reference or as value, all data is contained in the object).

#include<iostream>
using namespace std;

const size_t MSIZE = 3;

template<typename T, size_t S> class squareMatrixT
{
    T matrix[S][S];
    typedef T (&oneDimRefT)[S];   // the interesting part

    public:

    const size_t size = S;
    oneDimRefT operator[](size_t i) { return matrix[i]; }
};

typedef squareMatrixT<int,MSIZE> sqm3T;   // convenience and brevity

void init(sqm3T &sqm3)
{
    for(int i=0; i<sqm3.size; i++)
    {
        for(int j=0; j<sqm3.size; j++)
        {
            sqm3[i][j] = i*sqm3.size + j;
        }
    }   
}

void  print(sqm3T sqm3)
{
    for(int i=0; i<sqm3.size; i++)
    {
        for(int j=0; j<sqm3.size; j++)
        {
            cout << sqm3[i][j] << ' ';
        }
        cout << endl;
    }   
}

int main(void)
{
    sqm3T sqm33main;
    init(sqm33main);   // by ref
    print(sqm33main);  // by val
    return 0;
}

Question 2

MSDN suggests that std::tr1::array is available. This eases usage of C-style arrays by wrapping them inside a class. This way not only simplifies the syntax (at least to me), but most importantly allows assigning, copying and therefore also returning them from functions:

template<typename T, size_t S>
Someclass<T, S>::Someclass(std::tr1::array<S, std::tr1::array<S, T> > const&);

template<typename T, size_t S>
std::tr1::array<S, std::tr1::array<S, T> > Someclass<T, S>:doSomething();

This is very close to what you wished, but differs from the C# solution in that it requires the size of the matrices to be known at compile time.

If you wish to allow for matrices whose sizes are known only at runtime, you will need to provide your own container (or find a 3rd party library providing a class for dynamic square matrices), possibly based on std::vector along the lines of the following:

template<typename T>
class matrix {
    std::vector<T> data;
    size_t size;

public:
    matrix(size_t size) : data(size * size), size(size) { }
    T& operator()(size_t i, size_t j) { return data.at(i * size + j); }
    T const& operator() const(size_t i, size_t j) { return data.at(i * size + j); }
    T& operator[](std::pair<size_t, size_t> p) { return data.at(p.first * size + p.second); }
    T const& operator[](std::pair<size_t, size_t> p) const { return data.at(p.first * size + p.second); }
};

Note, that overloading the array subscript operator to take more than one argument is sadly not allowed in c++.