Given a number $n \in \Bbb{N}$. In how many ways can $n$ be written as $\prod_{i=1}^{k}n_i$ such that $n1|n2|\ldots|n_k|n$? I came to this problem in a question a few days ago. I have not found any duplicates so I assume there are none. In fact, the proposed problem is equivalent to "How many abelian groups are there with order $n$?". First of all I tried to solve the problem where $n=p^a$ with $p$ prime. Looking at some examples like $p=p^{i}p^{i.j}p^{a-i-i.j}$ I came to the conclusion that each solution corresponds to the sizes of a partition of a set with $a$ elements (in the example the partition sizes are $(i,i.j,a-i-i.j)$, and vice versa. Then I tried some solutions where $n=p^aq^b$ where $p$ and $q$ prime and coprime. For every solution we write the exponents of $p$ and $q$ in the factors of the decomposition in a matrix where the rows are these exponents per factor and the columns correspond to the exponents of $a$ resp. $b$. An example Let $n=72$ then a solution is $72=12.6$ (I reversed the order on purpose). This would give rise to a matrix of the form $\bigl( \begin{smallmatrix} 2 & 1 \\ 1 & 1 \end{smallmatrix} \bigr)$, because $12=2^2.3^1$ and $6=2^1.3^1$. In the same way the decomposition $72=6.6.2$ corresponds to the matrix $\bigl( \begin{smallmatrix} 1 & 1 \\ 1 & 1 \\ 1 & 0 \end{smallmatrix}\bigr)$. The only conditions on these matrices is that the sum of the first column is $3$, because $2$ occurs with exponent $3$ in $72$, that the sum of the elements of the second colum is $2$ and that the elements of each column are non increasing. I suspect that, if I order these partitions by decreasing size I can fill in for the first column any list of sizes of a partition of $3$ elements, and any list of sizes of a partition of $2$ elements in the second column. This would imply that there are $3.2=6$ abelian groups of order $72$. If my suspicion is right then for any number $n=p_1^{e_1}p_2^{e_2}.\ldots.p_k^{e_k}$ there are $\prod_{i=1}^{k}P(e^i)$ abelian groups of order $n$ with $P(a)$ the number of partitions of $a$ elements. 
 A: If you write "as the product of k numbers", and you found the answer for $p^a$, then I'd say the answer for $p^a q^b r^c$ would be the product of the answers for $p^a$, $q^b$ and $r^c$. 
If you exclude trivial answers where $n_1 = 1$, then the largest k would be the largest of the prime exponents in the factorisation. So all you need is f (a, k) = "number of ways to write a as the sum of k non-decreasing non-negative integers"; you factor $n>1$ into the product of powers of primes $p_i^{a_i}$, let K = largest of the $a_i$, then for each $1 ≤ k ≤ K$ calculate the product of $f (a_i, k)$ over all i in the factorisation, and add the products. 
A: The following remarks  establish an exact formula for  these chains in
terms of the cycle index of  the symmetric group on $v$ elements, with
$v$ being the  maximal exponent of a prime  in the prime factorization
of $n.$
Suppose $$n  = \prod_{q=1}^m p_q^{v_q}$$ is the  prime factoriation of
$n.$ Now if we consider the  exponents $w_1, w_2, w_3, \ldots, w_k$ of
$p_q$ dividing $n_1, n_2, n_3, \ldots, n_k$ then we must have
$$w_1\le w_2\le w_3\le\cdots\le w_k
\quad\text{and}\quad \sum_{j=1}^k w_k = v_q.$$
The  number  of such  $k$-tuples  can  be  calculated with  the  Polya
Enumeration Theorem (PET) and has the species equation
$$\mathfrak{M}_{=k}
(\epsilon + \mathcal{Z} + \mathcal{Z}^2 + \mathcal{Z}^3 + \cdots).$$
This translates to the generating function equation
$$[z^{v_q}] Z(S_k)\left(\frac{1}{1-z}\right).$$
It follows that the count $P_{n,k}$ of such tuples/chains is given by
(this is the promised exact formula)
$$P_{n,k} = 
\prod_{q=1}^m [z^{v_q}] Z(S_k)\left(\frac{1}{1-z}\right).$$
This includes  chains that have a  prefix of a string  of one factors,
which we  are apparently not counting here.  These inadmissible chains
among  the chains  counted by  $P_{n,k}$ have  the property  that they
correspond  bijectively  to the  chains  from  $P_{n,k-1}.$ Hence  the
desired value $Q_{n, k}$ is given by
$$Q_{n,k} = P_{n,k} - P_{n, k-1}.$$
Here are some examples of the cycle indices that are used:
$$Z(S_3) = 1/6\,{a_{{1}}}^{3}+1/2\,a_{{2}}a_{{1}}+1/3\,a_{{3}},$$ 
$$Z(S_4) = 1/24\,{a_{{1}}}^{4}+1/4\,a_{{2}}{a_{{1}}}^{2}
\\+1/3\,a_{{3}}a_{{1}}+1/8\,{a_{{2}}}^{2}+1/4\,a_{{4}}$$ 
and
$$Z(S_5) = {\frac {{a_{{1}}}^{5}}{120}}+1/12\,a_{{2}}{a_{{1}}}^{3}\\
+1/6\,a_{{3}}{a_{{1}}}^{2}+1/8\,a_{{1}}{a_{{2}}}^{2}\\
+1/4\,a_{{4}}a_{{1}}+1/6\,a_{{2}}a_{{3}}+1/5\,a_{{5}}.$$ 
It  is  not  difficult to  see  that  the  maximum  $k$ such  that  an
admissible chain for $n$ exists is the maximal exponent call it $v$ of
a prime $p$  in  the  prime  factorization  of  $n.$  (This  was  also 
observed by the  first responder  above.)  The  chain  is  obtained by
distributing $v_q$ copies of  $p_q$ into the $v$ slots  starting  from
the right, with divisibility being determined from the left. The value
in a slot is the product of the primes that have been distributed into
it.

It follows that  the total number $T_n$ of  admissible chains is given
by
$$T_n = \sum_{k=1}^v Q_{n,k}
= \sum_{k=1}^v (P_{n,k} - P_{n,{k-1}})
= P_{n, v}
= \prod_{q=1}^m [z^{v_q}] Z(S_v)\left(\frac{1}{1-z}\right).$$
This gives the sequence
$$1, 1, 1, 2, 1, 1, 1, 3, 2, 1, 1, 2, 1, 1, 1, 5, 1, 2, 1, 2, 1,\\
1, 1, 3, 2, 1, 3, 2, 1, 1, 1, 7, 1, 1, 1, 4, 1, 1, 1, 3, 1, 1, 1,\\
2, 2, 1, 1, 5, 2, 2, 1, 2, 1, 3, 1, 3, 1, 1, 1, 2, 1, 1, 2, 11,\ldots$$
which is OEIS A000688
where we find that indeed it  counts the number of Abelian groups as
surmised by the OP.
The Maple  code for  these including a  verification of the  values of
$Q_{n,k}$ by brute force (routine Q_ex) is as follows:

with(numtheory);
with(combinat);

pet_cycleind_symm :=
proc(n)
local p, s;
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;


pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;

    res := ind;

    polyvars := indets(poly);
    indvars := indets(ind);

    for v in indvars do
        pot := op(1, v);

        subs1 :=
        [seq(polyvars[k]=polyvars[k]^pot,
             k=1..nops(polyvars))];

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;

P :=
proc(n, k)
    option remember;
    local pv, p, ind, gf;

    ind := pet_cycleind_symm(k);
    gf := pet_varinto_cind(1/(1-z), ind);

    p := 1;

    for pv in ifactors(n)[2] do
        p := p*coeftayl(gf, z=0, pv[2]);
    od;

    p;
end;

Q :=
proc(n, k)
    option remember;
    P(n,k) - P(n, k-1);
end;

T :=
proc(n)
    option remember;
    local pv, mx;

    mx := -1;
    for pv in ifactors(n)[2] do
        if pv[2] > mx then
            mx := pv[2];
        fi;
    od;

    P(n, mx);
end;

chooserep :=
proc(data, sofar, k, res)
    if k=0 then
        res[sofar] := 1;
        return;
    fi;

    if nops(data) > 0 then
        if nops(sofar) = 0 or data[1] mod sofar[-1] = 0 then
            chooserep(data, [op(sofar), data[1]], k-1, res);
        fi;
        chooserep([op(2..nops(data), data)], sofar, k, res);
    fi;
end;

Q_ex :=
proc(n, k)
    option remember;
    local res, tupl, divchoose;

    res :=0;

    divchoose := table();
    chooserep(sort(convert(divisors(n) minus {1}, list)),
              [], k, divchoose);

    for tupl in [indices(divchoose, 'nolist')] do
        if convert(tupl, `*`) = n then
            res := res+1;
        fi;
    od;

    res;
end;

Addendum. One possible improvement of Q_ex would be not to collect all chains and then check them for admissibility and instead check admissibility as they are generated. The code looks the way it does because it started with the Maple command choose as its paradigm.
