calculus and analysis - How to deal with complicated Gaussian integrals in Mathematica?

As we know, for most Gaussian integrals, we can get the analytical result. Now I have many Gaussian integrals to treat, which have the following general form,

Integrate[
   x1^n1 x2^n2 x3^n3 ... xd^
     nd Exp[-{x1, x2, x3 ... xd}.matrix.{x1, x2, x3 ... xd} + c1 x1 + 
        c2 x2 ... cd xd], {x1, -Infinity, Infinity},
     {x2, -Infinity, Infinity} ... {xd, -Infinity, Infinity}]. 

The matrix is a symmetrical positive definite matrix.

We can get the analytical integral result for this general form. However, it is a difficult task for Mathematica to deal with the similiar expressions, which takes so much time. I have written a package for treating this problem.

But, it is still not sufficient because I need to deal with lots of these kinds of integral expressions. So the computation time is important.

So have you ever seen the package for doing this kind of Gaussian integral? Or can you give me some suggestions about this problem? Thank you very much!

Thanks all of you. For the concept MultinormalDistribution, it is very helpful for me.

By the way, How can I transform the following integral expression to the standard form of MultinormalDistribution in a simple way?

vec1 := {r1x, r1y, r1z};
vec2 := {r2x, r2y, r2z};
vect := {r1x, r1y, r1z, r2x, r2y, r2z};
vcof := {c1, c2, c3, c4, c5, c6};
Integrate[Exp[-vect.mat.vect + vcof.vect] (vec1.vec2)^n (vec1.vec1)^
l (vec2.vec2)^m, {r1x, -Infinity, Infinity}, {r1y, -Infinity, 

Infinity}, {r1z, -Infinity, Infinity}, {r2x, -Infinity, 
Infinity}, {r2y, -Infinity, Infinity}, {r2z, -Infinity, Infinity}];

The mat is a symmetrical positive definite sparse matrix. And n, l, m are integers.
Thanks!