Python3, 659 bytes

Long, but fast.

def R(n,k):
 q=[(n,[])]
 for n,p in q:
  if len(p)>k:continue
  if len(p)==k and not n:yield tuple(sorted(p));continue
  for j in range(1,n+1):q+=[(n-j,p+[j])]
def M(e,t,I,l):
 q=[(l,{*I},[])]
 for l,I,m in q:
  if[]==l:yield m;continue
  q+=[(l[1:],I-{i},m+[(i,l[0])])for i in I if e[i]+l[0]<=t[i]]
def f(t):
 q,r=[([0 for _ in t],[])],0
 for e,v in q:
  if e==t:r+=1;continue
  if(l:=[i for i,a in enumerate(e)if a<t[i]]):
   i,*I=l
   for j in range(1,len(I)+1):
    for l in{*R(t[i]-e[i],j)}:
     for U in map(eval,{str(sorted(K))for K in M(e,t,I,[*l])}):
      E=[*e];E[i]=t[i]
      for x,y in U:E[x]+=y
      q+=[(E,v+[(i,x,y)for x,y in U])]
 return r

Try it online!

Haskell, 57 bytes

f(h:t)=sum[f l|l<-mapM(\n->[0..n])t,sum t-sum l==h]
f _=1

Try it online!

Haskell, 75 bytes

f(h:t)=sum$f<$>h%t
f[]=1
r%(h:t)=[h-i:x|i<-[0..r],x<-(r-i)%t]
r%_=[[]|r==0]

Try it online!

It's faster (and easier to understand) when i<-[0..r] is replaced by i<-[0..min r h]. Also, it's faster to give the degrees (whose order obviously doesn't matter) in increasing order.

Haskell, 127 bytes

e=[]:e
h a b=mapM id$b<$a
f a=length[0|x<-h a$h a[0..maximum a],x==foldr(zipWith(:))e x,all(<1)$zipWith(!!)x[0..],map sum x==a]

Try it online!

Very very inefficient (superexponential), but I believe it should work in theory. Exceeds Tio's 60 second limit with [3, 3, 1, 1] as input, but works with [1, 1, 1, 1].

Considers the adjacency matrices of the multigraphs.

h a$h a[0..maximum a] is a list of all square matrices with entries between 0 and maximum a, with size length a.

x==foldr(zipWith(:))e x checks if a matrix is symmetric.

all(<1)$zipWith(!!)x[0..] checks that the diagonal entries are zero, since there are no loops.

map sum x==a checks that each vertex has the proper degree.