FR: Matrix inputs to multi_normal_*

It would be really nice to have signatures for the multi_normal_* family that accept matrix inputs for at least the data and means arguments, possibly with a new argument that designates whether they should be treated as row-wise or column-wise manner. For example, see below for a centered-parameterized hierarchical model (of the same core structure as SUG1.13), where I currently have to do two format conversions, one to convert a matrix to an array of row-vectors for input to multi_normal_cholesky, then another conversion of the first argument to a matrix for use with rows_dot_product. The latter is made easier by #931, but it would be nice to skip it altogether.

data{

	// nXg: number of cols in the group-level predictor matrix
	int<lower=2> nXg ;

	// rXg: number of rows in the group-level predictor matrix
	int<lower=nXg> rXg ;

	// Xg: group-level predictor matrix
	matrix[rXg,nXg] Xg ;

	// nI: number of individuals
	int<lower=nXg> nI ;

	// iXg: which group each individual is associated with
	array[nI] int<lower=1,upper=rXg> iXg ;

	// nXq: number of cols in the condition-level predictor matrix
	int<lower=2> nXq ;

	// rXq: number of rows in the condition-level predictor matrix
	int<lower=nXq> rXq ;

	// Xq: condition-level predictor matrix
	matrix[rXq,nXq] Xq ;

	// iXq: which individual is associated with each row in Xq
	array[rXq] int<lower=1,upper=nI> iXq ;

	// nY: num entries in the observation vectors
	int<lower=nXg*nXq*nI> nY ;

	// Y: observations
	vector[nY] Y ;

	// yXq: which row in Xq is associated with each observation in Y
	array[nY] int<lower=1,upper=rXq> yXq ;

}
parameters{

	// Y_sd: magnitude of observation-level variability
	real<lower=0> Y_sd ;

	// Z: coefficients associated with predictors (group-level, condition-level, & interactions)
	matrix[nXg,nXq] Z ;

	// iZq_sd: magnitude of variability among individuals within a group
	vector<lower=0>[nXq] iZq_sd ;

	// iZq_cholcorr: cholesky-factor of correlation structure associated with variability among individuals on influence of within-individual predictors
	cholesky_factor_corr[nXq] iZq_cholcorr ;

	// iZq: by-individual coefficients (centered parameterization)
	array[nI] row_vector[nXq] iZq ;

}
model{

	////
	// group-level structure
	////

	// standard-normal priors on all group-level coefficients
	to_vector(Z) ~ std_normal() ;

	// using the group predictors and coefficients, compute condition coefficients for each group
	// NOT SURE THIS IS THE MOST EFFICIENT WAY TO DO THIS (ESP. GIVEN NECESSARY LATER CONVERSION)
	matrix[rXg,nXq] gZq ;
	for(this_nXq in 1:nXq){
		gZq[,this_nXq]= rows_dot_product(
			rep_matrix(to_row_vector(Z[,this_nXq]),rXg)
			, Xg
		) ;
	}

	//convert gZq from matrix to array of row-vectors
	array[rXg] row_vector[nXq] gZq_arr ;
	for(this_rXg in 1:rXg){
		gZq_arr[this_rXg] = gZq[this_rXg] ;
	}


	////
	// individual-level structure
	////

	// positive-standard-normal priors on magnitude of variability among individuals within a group
	iZq_sd ~ std_normal() ;

	// flat prior on correlations
	iZq_cholcorr ~ lkj_corr_cholesky(1) ;

	// multi-normal structure for iZq
	iZq ~ multi_normal_cholesky(
		gZq_arr[iXg]
		, diag_pre_multiply(iZq_sd, iZq_cholcorr)
	) ;

	//convert iZq from array of row-vectors to matrix
	// (hopefully replaceable by `to_matrix(iZ)` soon,
	// see: https://github.com/stan-dev/cmdstan/issues/1015 )
	matrix[nI,nXq] iZq_mat ;
	for(this_nI in 1:nI){
		iZq_mat[this_nI] = iZq[this_nI] ;
	}

	// using the indivividual condition coefficients and predictors, compute
	// values for each individual
	vector[nXq] iZq_dot_Xq = rows_dot_product( iZq_mat[iXq] , Xq ) ;

	////
	// observation-level structure
	////

	// prior peaked around .8 for magnitude of observation-level variability
	Y_sd ~ weibull(2,1) ;

	Y ~ normal(
		iZq_dot_Xq[yXq]
		, Y_sd
	) ;

}

Start with the multi_normal_* signatures and the matrix/array conversions in the example; the issue does not name implementation files or tests. Define the supported matrix inputs and row-wise versus column-wise behavior, then verify that the shown hierarchical model no longer needs those conversions.