Nekomata + -e, 7 bytes

ᴶ∑Ťᴶ∑←ž

Attempt This Online!

ᴶ∑Ťᴶ∑←ž
ᴶ∑          Non-deterministically split the input into parts and sum each part
  Ť         Transpose
   ᴶ∑       Again non-deterministically split and sum
     ←      Decrement each element by 1
      ž     Check that all elements are zero

-e checks if there is a valid solution.

Python3, 559 bytes

P=lambda a,V,l:[[[i.pop(0)for _ in range(t)]for t in a]for i in V[:l]]
W=lambda a,B,l,F=0:all(sum(map(sum,j))==1 if F else sum(map(sum,j))<=1 for j in zip(*P(a,[[*i]for i in B],l)))
def f(b):
 q,s=[[1]],[]
 for a in q:
  if sum(a)==len(b[0]):
   Q=[[1]]
   for A in Q:
    X=[*b];C=[[X.pop(0)for _ in range(y)]for y in A]
    if sum(A)==len(b):
     if all(W(a,B,len(B),1)for B in C):return 1
     continue
    if all(W(a,B,len(B))for B in C):Q+=[A[:-1]+[A[-1]+1],A+[1]]
   continue
  if W(a,b,1):q+=(o:=[p for p in[a[:-1]+[a[-1]+1],a+[1]]if p not in s]);s+=o

Try it online!

Jelly, 20 bytes

ŒṖS€€ỊȦ$ÐfZ€µ⁺€Ȧ€€FS

This is still a brute-force solution, but it's quite a bit faster than my other answer – which cannot cope with the last two test cases on TIO – and handles all test cases in ~4 seconds.

Try it online!

How it works

ŒṖS€€ỊȦ$ÐfZ€µ⁺€Ȧ€€FS  Main link. Argument: M (matrix, array of rows)

ŒṖ                    Compute all partitions, i.e., all groupings of M's rows.
  S€€                 Map sum over all individual groupings, collapsing the grouped
                      rows into a single row.
        Ðf            Filter; keep only those partially collapsed matrices for
                      which the link to the left returns a truthy value.
       $                Group the two links to the left into a monadic chain.
     Ị                    Insignificant; map 0 and 1 to 1, greater integers to 0.
      Ȧ                   All; return 1 iff the matrix contains no zeroes.
          Z€          Zip/transpose all kept matrices,
            µ         Combine all links to the left into a monadic chain.
             ⁺€       Duplicate the chain and map it over the individual results
                      from the first call. We now have all possible combinations
                      of row and column groupings (represented by the corresponding
                      matrices of collapsed rows and columns) that do not have a
                      2 anywhere. However, they still may contain zeroes.
               Ȧ€€    Map the all atom over the matrices, returning 1 only for
                      matrices that consist entirely of ones.
                  FS  Flatten and sum, counting the number of valid divisions.

Haskell, 116 bytes

import Data.List
m(a:b)=[a:e|e<-m b]++[zipWith(+)a d:e|d:e<-m b];m e=[e]
d=(any$any$all$all(==1)).map(m.transpose).m

Try it online!

Jelly, 22 bytes

ŒṖḅ1S€€=1Ȧ€
SŒṖ⁸ṁþÇ€FS

Try it online!

Pyth, 30 29 26 24 23 bytes

sm.Asmmq1ssbCkds./MC./M

Try it online.

I'm sure this will get shorter. This is O(2^mn), where m and n are the width and height of the array, but completes the last two test cases in 45 seconds on my laptop on battery (i5-5200U with limited performance).

Outputs the number of solutions.

Explanation

Five-dimensional arrays are really fun to work with._^</sarcasm> You are not supposed to understand how this works even with the explanation.

                    ./M    Find all partitions of each row. Now we have a list of rows,
                           each containing the ways to split the row, each containing
                           the parts of the split (3D).
                   C       Transpose. Now we have a list of ways to split the columns,
                           each containing the rows, each containing the parts of the
                           row (3D).
                ./M        Find all partitions of each row list. Now we have a list of
                           ways to split the columns, each containing the ways to split
                           the rows, each containing the bunch of rows, each containing 
                           the rows in the bunch, each containing the parts of the row
                           (6D).
               s           Combine the ways to split rows & columns into one array (5D).
 m            d            Do the following for each way to split rows & columns (4D):
     m       k                 Do the following for each bunch of rows (3D):
            C                      Transpose the array. We now have a list of column
                                   groups, each containing the row parts (3D).
      m    b                       Do the following for each column group (2D):
          s                            Combine the row parts in the column group. We now
                                       have the list of cells in this row/column group
                                       (1D).
         s                             Sum the cells.
       q1                              Check if the sum is one.
                                   We now have the list of booleans that tell if each
                                   row/column group is valid (1D).
                               We now have the 2D list of booleans that tell if each
                               row/column group in each bunch of rows is valid. (2D)
    s                          Combine the 2D list of booleans to 1D.
  .A                           Check if all values are truthy; if the split is valid.
                           We now have the validity of each split.
s                          Sum the list to get the number of valid solutions.