Array input arguments in nonlinear expressions

Context: Yesterday I had a chat to @rluce about nonlinear expressions, particularly as they relate to Gurobi's upcoming nonlinear interface. We broadly agree on scalar nonlinear functions, and he had some preliminary ideas for vector-inputs.

cc @ccoffrin and @chriscoey for thoughts.

Our ultimate goal is to support examples like the following:

```julia
using JuMP, Ipopt
model = Model(Ipopt.Optimizer)
@variable(model, x[1:2, 1:2])
@objective(model, Max, log(det(X)))
# or perhaps the easier to manage;
@objective(model, Max, log_det(X))
```
```julia
using JuMP, Ipopt
model = Model(Ipopt.Optimizer)
@variable(model, x[1:3] >= 1)
p = 2
@objective(model, Min, norm(x, p))
```

Another key consumer of this would be https://github.com/jump-dev/MiniZinc.jl.

## Changes required in JuMP

* We'd need to support `Array` in `GenericNonlinearExpr` and their mapping to `moi_function`. This seems pretty easy.

## Changes required in MOI

1. We'd need to support `Array` in `MOI.(Scalar,Vector)NonlinearFunction`. This seems pretty easy.
2. We'd need to support `Array` in `MOI.Nonlinear` and be able to compute derivatives, etc.

The tricky part is all in 2.

We'd likely need some sort of `Node(NODE_VECTOR, parent, n)`.

But matrices are a bit more complicated. We'd need to encode `(rows, cols)`. One option would be to store the size as a packed `value::Int64`. That'd mean that we couldn't have matrices with side dimension greater than `typemax(Int32)`... but that seems okay.

```julia
encode(m::Int64, n::Int64) = encode(Int32(m), Int32(n))
encode(m::Int32, n::Int32) = reinterpret(Int64, (m, n))
decode(x::Int64) = Int64.(reinterpret(Tuple{Int32,Int32}, x))
```

`norm(x)` would look something like:
```julia
expr = Expression(
[
Node(NODE_CALL_MULTIVARIATE, -1, OP_NORM ),
Node(NODE_ARRAY, 1, 3 #= encode(3, 0) =#),
Node(NODE_VARIABLE, 2, 1 #= x1 =# ),
Node(NODE_VARIABLE, 2, 2 #= x2 =# ),
Node(NODE_VARIABLE, 2, 3 #= x3 =# ),
],
[],
);
```

`log_det(X)` would then look something like:
```julia
expr = Expression(
[
Node(NODE_CALL_MULTIVARIATE, -1, OP_LOGDET ),
Node(NODE_ARRAY, 1, 8589934594 #= encode(2, 2) =#),
Node(NODE_VARIABLE, 2, 1 #= x11 =# ),
Node(NODE_VARIABLE, 2, 2 #= x21 =# ),
Node(NODE_VARIABLE, 2, 2 #= x12 =# ),
Node(NODE_VARIABLE, 2, 3 #= x22 =# ),
],
[],
);
```

Once you have the data structure, it seems to me that the AD should follow fairly well.

## Output arguments

Still absolutely no idea.

Start with GenericNonlinearExpr, MOI.(Scalar,Vector)NonlinearFunction, MOI.Nonlinear, and the proposed Node and Expression structures described in the issue. Trace how scalar nonlinear expressions and derivatives are currently represented, then determine how array and matrix inputs should be represented and differentiated. Done means the design and implementation support examples such as norm(x, p) and log_det(X), with tests covering their derivatives.

Written by the indexing model from the issue text.