The given condition can be consider the requirement $AX+s=b$ of an ILP in standard form. However, as noted in the question, there is no function to maximize. Well, that is not big obstacle at all. We can specify any function to maximize. For example, we can maximize the amount of $m_1$. Or maximize the amount of $m_2$.
There are various algorithms that solve ILP. However, the problem of integer linear programming is NP-hard. It is not clear which algorithm might be efficient for the current problem.
Another approach with a linear algorithm
Let us use math a bit. What we need is the method of elimination of variables and the extended Euclidean algorithm that solves linear Diophantine equation of two variables. Throughout the explanation below, I will skip various easier cases so that we can focus on the main idea. All variables will be integers.
Let the amount of $m_i$ be $x_i\in\Bbb N_{>0}$. The requirements are
$$ a_1x_1+a_2x_2+a_3x_3=A\\
c_1x_1+c_2x_2+c_3x_3=C$$
Consider $x_1$ a constant, we can solve $x_2, x_3$.
$$rx_3=p_2x_1-q_2\\
rx_2=p_3x_1-q_3$$
where $p_3=c_1a_3-c_3a_1$, $\ q_3=Ca_3-Ac_3$, $\ p_2=-(c_1a_2-c_2a_1)$, $\ q_2=-(Ca_2-Ac_2)$, $\ r=a_2c_3-c_2a_3$.
I will skip the easy case when $r=0$. Assume $r\not=0$.
The original problem is equivalent to finding a positive integer $x_1$ such that both $p_2x_1-q_2$ and $p_3x_1-q_3$ are positive multiples of $r$. Thanks to the elimination of variables $x_2$ and $x_3$, we will deal with a single integral variable $x_1$ only.
Fact. Let $a, b, n$ are integers, $a\not=0$, $b\not=0$. We can use the extended Euclidean algorithm or here to find $d=\gcd(a,b)$ as well as an integral solution $(x,y)=(x_0,y_0)$ of the equation $$ax+by=d.$$
- If $n$ is a multiple of $d$, all integral solutions to $ax+by=n$ are given by $$(x_0\frac nd+\frac bdk, y_0\frac nd-\frac adk),$$ where $k\in\Bbb Z$.
- If $n$ is not a multiple of $d$, there is no integral solution to $ax+by=n$.
Consider the condition that $p_2x_1-q_2$ is a positive multiple of $r$. I will skip the easy case when $p_2=0$. This condition asks us to solve the Diophantine equation $p_2x_1+ry_1=q_2$ with positive $x_1$ and $y_1$. I will skip the easy case when $q_2$ is not a multiple of $\gcd(p_2, r)$. So all integral solution to $p_2x_1+ry_1=q_2$ are $(x_1, y_1)=(s_0+s_1k_1, s_2-s_3k_1)$, where $k_1\in\Bbb Z$ for some constant $s_0, s_1, s_2, s_3\in\Bbb Z$ that can be computed according to the fact above. The condition that $x_1$ and $y_1$ are positive will mean $k_1$ must be in some interval (possible empty or infinite). In summary, $x_1=s_0+s_1k_1$ for some integer constant $s_0, s_1$ where integer $k_1$ is in some known interval.
Similarly, with some easier cases skipped, the condition that $p_3x_1-q_3$ is a positive multiple of $r$ will mean $x_1=t_0+t_1k_2$ for some integer constant $t_0, t_1$ where integer $k_2$ is in some known interval.
Combining both conditions, we find the problem is reduced to finding integers $k_1, k_2$ in some know intervals respectively such that $s_0+s_1k_1=t_0+t_1k_2$. The last equation is again a linear equation of two variables $(k_1, k_2)$, which can be solved by the fact above.
If we implement the approach above properly, the algorithm will have $O(\log(\max(c_1, c_2, c_3, a_1, a_2, a_3, A, C)))$ time-complexity and $O(1)$ space-complexity.
