Outline
A. I will first clarify the question.
B. Then I discuss an approach which unfortunately does not work, but which might nevertheless give some insight.
C. Finally I discuss a some non-UFD GCD domain $S_i$ and sketch proofs that they still fulfill the property that the OP is interested in, i.e. are not "counterexamples" in the sense of the OP. This is however not for the "reason of being an UFD".
My answer was informed by an answer by Bill Dubuque to another question.
$
\require{begingroup}\begingroup
\DeclareMathOperator{\Frac}{Frac}
\DeclareMathOperator{\lcm}{lcm}
\DeclareMathOperator{\mult}{mult}
$
General notation
$R$ will always be an integral domain, i.e. a commutative ring with $1$ that has no zero divisors.
For $a,b\in R$ which have a gcd we know: $\gcd(a,b)\mid a$ (the "greatest common divisor" is a "divisor", as it should be) and thus there is $a'\in R$ such that $a'\gcd(a,b)=a$. We will abbreviate this element $a'=\frac{a}{\gcd(a,b)}$, even though the element $\frac{1}{\gcd{a,b}}\in\Frac(R)$ in the fraction field might not be in $R$.
Question from the title:
Given $a,b\in R$ a GCD domain it is always possible to find $d_1,d_2\in R$ which have $d_1d_2=\lcm(a,b)$ and that are divisors in the sense $d_a:=d_1\mid a$ and $d_b:=d_2\mid b$. Set e.g. $d_a=a, d_b=\frac{b}{\gcd(a,b)}$ or symmetrically $d_a=\frac{a}{\gcd(a,b)}, d_b=b$.
See also this related question and its answer.
Question from the body:
Given $a,b\in R$ a GCD domain, is it always possible to find $d_a:=d_1, d_b:=d_2\in R$ which satisfy $d_ad_b=\lcm(a,b),\quad$ $d_a\mid a,\quad$ $d_b\mid b$ and $\gcd(d_a,d_b)=1$?
Currently I do not have an answer for that, but I will share some deliberations. Notation: We say the pair $(d_a,d_b)$ are gl-divisors of the pair $(a,b)$ where "gl-divisor" abbreviates "gcd-less lcm-full divisors". A GCD domain that admits gl-divisors for every pair $(a,b)\in R^2$ is called a gl GCD domain.
One could try the ansatz $d_a=a, d_b=\frac{b}{\gcd(a,b)}$ again, but that would even fail for the nicest of domains. A counterexample for $\mathbb{Z}$: $a=2^2\cdot 3=12, b=2\cdot 3^2=18$ then $\gcd\left(a,\frac{b}{\gcd(a,b)}\right)=\gcd(12,3)=3\neq 1$.
In a UFD one could use the prime factorization to remedy this as pointed out by the OP. Here is another take on the situation which unfortunately will boil down to the same result ("it works for UFDs") in the end, but still might give some inspiration on how to go about a potential counterexample.
Let again $R$ be a general GCD domain, $a,b\in R$ arbitrary and $d_a=a, d_b=\frac{b}{\gcd(a,b)}$. We can factor $d_a=a=\gcd(a,b)\cdot\frac{a}{\gcd(a,b)}$. The problem above was $\gcd\left(a,\frac{b}{\gcd(a,b)}\right)\neq 1$. We can pin down the problem further. Set: $A:=\frac{a}{\gcd(a,b)}, B:=\frac{b}{\gcd(a,b)}$ and let C be a common divisor of $A$ and $B$, i.e. there are $C_A,C_B\in R$ such that $A=C_A\cdot C, B=C_B\cdot C$. Then we have
$$
a=A\cdot\gcd(a,b)=C_A\cdot C\cdot\gcd(a,b)\quad\text{ and }\quad
b=B\cdot\gcd(a,b)=C_B\cdot C\cdot\gcd(a,b).
$$
This implies $C\cdot\gcd(a,b)\mid\gcd(a,b)$ and therefore $C\in R^{\times}$ must be a unit (recall that gcds are only unique up to units). But then
$$
\gcd\left(\frac{a}{\gcd(a,b)},\frac{b}{\gcd(a,b)}\right)=1.
$$
So the failure of our ansatz cannot come from the factor $\frac{a}{\gcd(a,b)}$ of $a$, but must come from the factor $\gcd(a,b)$ of $a$.
This suggests the following procedure: Set $d_{a,1}=a, d_{b,1}=\frac{b}{\gcd(a,b)}$. If $\gcd(d_{a,1},d_{b,1})=1$ we are done. Otherwise set:
$$
d_{a,2}=\frac{d_{a,1}}{\gcd(d_{a,1},d_{b,1})}\quad\text{ and }\quad d_{b,2}=d_{b,1}\cdot\gcd(d_{a,1},d_{b,1}).
$$
Unfortunately $\gcd\left(d_{a,2},d_{b,2}\right)\neq 1$ is still possible. E.g. for $R=\mathbb{Z}$ and $a=2^3\cdot 3^2=72,b=2^2\cdot 3^3=108$ we have $d_{a,1}=72,d_{b,1}=3$ and $d_{a,2}=24,d_{b,2}=9$, but still $\gcd\left(d_{a,2},d_{b,2}\right)=3\neq 1$.
Next we iterate this process by setting
$$
d_{a,n+1}=\frac{d_{a,n}}{\gcd(d_{a,n},d_{b,n})}\quad\text{ and }\quad d_{b,n+1}=d_{b,n}\cdot\gcd(d_{a,n},d_{b,n}).
$$
This defines a sequence of ascending ideal on the "$a$-side":
$$
(d_{a,1})\subset(d_{a,2})\subset\ldots\subset(d_{a,n})\subset\ldots
$$
If $R$ was Noetherian, i.e. the ACC (ascending chain condition) holds for all ideals, then such a sequence had to stabilize, i.e. $\gcd(d_{a,m},d_{b,m})=1$ for some $m\in\mathbb{N}_{>0}$ and $R$ was a gl GCD domain. Even the weaker condition of ACC for principal ideals would suffice. Unfortunately GCD domains that satisfy ACC for principal ideals are exactly UFDs, cf. wikipedia.
But that does of course not mean that non-UFD GCD rings necessarily cannot be gl GCD domains. It just means that the method sketched above may fail - it might even work for some non-UFD GCD rings, although I do not have an example.
An example: The ring $S_1$ of entire functions (on $\mathbb{C}$)
Any UFD is a Bézout domain. Any integral domain is a Bézout domain if and only if it is Prüfer domain and a GCD domain at the same time. An example of a Bézout domain is the ring $S$ of holomorphic functions on $\mathbb{C}$.
By the Weierstraß factorization theorem we can construct a holomorphic functions $f,g$ using elementary factors $E_{k}$ as follows:
$$
f(z)=\prod\limits_{i=1}^\infty E_{p_n}\left(\frac{z}{a_n}\right),\quad
g(z)=\prod\limits_{i=1}^\infty E_{p_n}\left(\frac{z}{b_n}\right),
$$
where $(a_n)=(1,2,3,\ldots)$ is the sequence of positive integers, $(b_n)=(1,2,2,3,3,3,4,\ldots,k-1,\stackrel{k-\text{fold}}{\overbrace{k,\ldots,k}},k+1,\ldots)$ and $(p_n)=(2,2,2,\ldots)$ the constant sequence $2$. The idea is that $f$ has zeros exactly at the positive integers $z_i=i\in\mathbb{Z}_{>0}$ with constant multiplicity $1$ and $g$ has the same zeros but with increasing multiplicities $i$. Set
$$
a=a(z)=g\left(\frac{1}{2}(z+1)\right)\cdot g\left(\frac{1}{2}z\right)\cdot f\left(\frac{1}{2}(z+1)\right),
b=b(z)=g\left(\frac{1}{2}(z+1)\right)\cdot g\left(\frac{1}{2}z\right)\cdot f\left(\frac{1}{2}z\right)
$$
The zeros of $a,b$ are still at the integers, for an odd positive integer $z_o$ the multiplicity of $a$ is $1$ more than for $b$ - as formula $\mult_{z_o}(a)=\mult_{z_o}(b)+1$ - and vice versa for even positive integers $z_e$: $\mult_{z_e}(a)+1=\mult_{z_e}(b)$.
The procedure above produces $d_{a,n},d_{b,n}$ in step $n$ such that for an even positive integer $z_k=2k$ we have $\mult_{z_k}\gcd(d_{a,n},d_{b,n})>0$ if $k\geq 2^{n-2}$. Since we can take $k$ arbitrarily large we see that this approach will never terminate globally. However viewing each "finite place" $z_k$ isolated it does terminate in the sense, that $\gcd(d_{a,n},d_{b,n})$ will not vanish at that $z_k$ for all sufficiently large $n$.
However if we set
$$
d_a=g\left(\frac{1}{2}(z+1)\right)\cdot f\left(\frac{1}{2}(z+1)\right),
d_b=g\left(\frac{1}{2}z\right)\cdot f\left(\frac{1}{2}z\right)
$$
then $(d_a,d_b)$ are gl-divisors for $(a,b)$: let $c=\gcd(d_a,d_b)$. If $c$ was not a unit, then $c$ had a zero at some $z_0\in\mathbb{C}$. However all zeros of $d_a$ are at the odd positive integers and all zeros of $d_b$ are at the even positive integers, so there cannot be a common zero $z_0$. Thus $c$ is a unit and can be chosen to be $c=1$. The other properties are straightforward or in the case of the lcm-condition involve similarly checking that functions vanish to the right multiplicity at the right places.
The last paragraph motivates a proof that $S_1$, the ring of entire functions, is actually a gl GCD domain, which we quickly sketch. The zeros of entire functions are isolated thus countable. Given $a,b\in S_1$ let $\mathcal{Z}=\{z_1,z_2,z_3,\ldots\}$ be the indexed set of all zeros of $a$ or $b$. Partition $\mathcal{Z}=\mathcal{Z}_a\dot\cup\mathcal{Z}_b$ as follows: if $\mult_{z_i}(a)\geq\mult_{z_i}(b)$ then $z_i\in\mathcal{Z}_a$, otherwise $z_i\in\mathcal{Z}_b$. By the Weierstraß factorization theorem we can construct a holomorphic function $d_a$ that vanishes exactly on $\mathcal{Z}_a$ with the same multiplicities as $a$ (a convergence condition for the existence of $d_a$ has to be checked, but this follows quickly from the existence of a Weierstraß factorization for $a$). Set $d_b=\frac{\lcm(a,b)}{d_a}$ (this is possible by similar reasons as for the existence of $\frac{a}{\gcd(a,b)}$). $d_b$ vanishes exactly on $\mathcal{Z}_b$ with the same multiplicities as $b$. Thus any $(a,b)$ admits a pair of gl-divisors $(d_a,d_b)$, the properties are streight forward to check. In summary $S_1$ is a gl GCD domain.
Here are some further examples of interesting rings in this context:
the ring $\bar{\mathbb{Z}}$ of algebraic integers is a non-UFD Bézout domain, cf. wikipedia.
for any Bézout domain $R$ that is not a field, i.e. is different from its fraction field $F=\Frac(R)$, the ring $S=R+xF[x]$ (here $x$ is a variable independent of $R$) is a non-UFD Bézout domain, cf. wikipedia. In particular one can take for a UFD such as $R=\mathbb{Z}$.
for a field $k$ the monoid rings $k[\mathbb{Q}^+_0]$ and $k[\mathbb{R}^+_0]$ are non-UFD Bézout domain, cf. an answer to this question or wikipedia for a more general formulation.
if $R$ is a non-atomic GCD domain then $R[X_1,\ldots,X_m]$ is non-Bézout GCD domain for any $m\in\mathbb{N}_{>0}$, cf. wikipedia. E.g. the ring $R$ of holomorphic functions on $\mathbb{C}$ is a non-atomic GCD domain by the same source.
the book Hutchins, "Examples of Commutative Rings", gives examples of GCD-domains (called pdeudo-Bézout there): examples 33 (p. 67) is non-Bézout, for example 60 (p. 84) it is not stated if it is Bézout; it contains several more examples of Bézout domains some of them non-UFD.
the answers to this question also give examples worth considering, although the focus is an non-UFD domains. Some of these examples are GCD domains but it is not guaranteed that all of them are.
Remark: On the positive side the iterative procedure described above gives in principle an algorithm to compute gl-divisors in any UFD. The UFD has to be algorithmic, i.e. its elements can effectively stored and ring operations effectively computed. This algorithm might however not be very efficient.
Update: I have thought about the rings $S_2=\mathbb{Z}+x\mathbb{Q}[x]$, $S_3=\bar{\mathbb{Z}}$ and $S_4=k[\mathbb{Q}^+_0]$. I will sketch proofs that they all are gl GCD domains. The ideas will be fairly different and relay on the explicit ring given - there does not seem to be a unifying idea.
Example: The ring $S_2=\mathbb{Z}+x\mathbb{Q}[x]$
The irreducible elements $S_2$ are prime (that's true in any GCD domain, cf. Smith, "Fermat’s last theorem and Bezout’s theorem in GCD domains") and they are classified (up to units) as follows (cf. this question and the link to the duplicate therein): they are either
- a prime $p\in\mathbb{Z}$ or
- an irreducible polynomial $f\in\mathbb{Q}[x]$ with constant
coefficient $1$.
This gives rise to four different types of prime ideals:
- height 2 prime ideals of the form $(p)$,
- height 1 prime ideals of the form $(f)$,
- height 1 prime ideal $I_x=(x,\frac{x}{2},\frac{x}{3},\frac{x}{4},\ldots)$,
- the height 0 prime ideal $(0)$.
Note that the prime ideal $I_x$ in 3. is not finitely generated reflecting that $S_2$ is non-Noetherian. The ideals 1. give examples that Krull's principal ideal theorem does not need to hold in non-Noetherian rings. We also have $I_x\subsetneq(p)$ for every prime number.
From this one sees that we have "almost a unique prime factorization". If $a\in S_2$ has non-vanishing constant coefficient it can be written as a product of powers of irreducibles: $a=\prod_ip_i^{r_i}\cdot\prod_jf_j^{s_j}$. If its constant coefficient is $0$ we have $a=\frac{m}{n}x\cdot x^{u_1}\cdot\prod_jf_j^{s_j}$.
If we take $a,b\in R$ and want to find $\gcd(a,b)$ and $d_a,d_b$ we can use theses "factorizations". If $a,b$ are of the "same type" the gcd can be taken to have the factorization with the minimum of the exponents for $a$ resp. $b$ for the same factor, $d_a$ takes the factors where the exponent is at least as large as the matching exponent for $b$ and similarly for $d_b$. The interesting case is when the factorizations are of different type.
Let $a=\prod_ip_i^{r_i}\cdot\prod_jf_j^{s_j}$ and $b=\frac{m}{n}x\cdot x^{u_1}\cdot\prod_jf_j^{t_j}$. We deal with the $f_j$ factors as before, so we may wlog. assume $a=\prod_ip_i^{r_i}$ and $b=\frac{m}{n}x\cdot x^{u_1}$. In this case $\gcd(a,b)=a, \lcm(a,b)=b, d_a=a, d_b=\frac{m}{na}x\cdot x^{u_1}$.
Thus by explicit study of the multiplicative structure of $S_2$ we see that it is a gl GCD domain.
Example: The ring $S_3=\bar{\mathbb{Z}}$ of algebraic integers
This question asks about $S_3$ being Bézout. Although a proof is not explicitly given one of the answers illustrates the strategy with examples. Given explicit $a,b\in S_3$ the $\gcd(a,b)$ might not exist in the ring of integers $\mathfrak{o}_K$ of $K=\mathbb{Q}(a,b)$. But there is always an explicitly constructible finite extension field $L/K$ such that $\gcd(a,b)$ exists in $\mathfrak{o}_L$ - this does not mean that $\mathfrak{o}_L$ is itself a GCD domain or even Bézout, but only that enough prime ideals of $\mathfrak{o}_K$ become products of principle ideals of $\mathfrak{o}_L$ such that $\gcd(a,b)$ exists and $a,b$ admit suitable factorizations (not necessarily by irreducibles, e.g. if the extension $L$ is larger than necessary) facilitating the computation of $\gcd(a,b)$. Using this factorizations we should also be able to find suitable $d_a, d_b$ - I am sketchy here and did not work out the details, i.e. I did not prove this. If it does not work out let me know and let me know why and how.
In this example the ring $S_3$ has the structure of a union: $\bar{\mathbb{Z}}=\bigcup_i\mathfrak{o}_i$, where one even can assume that it is an ascending chain, i.e. $\mathfrak{o}_i\subset\mathfrak{o}_j$ for $i<j$ (assuming that the index set is suitably ordered). The idea is that for each pair $(a,b)$ there is some $\mathfrak{o}_i$ where $a,b$ factor nicely enough (although possibly not into irreducibles) such that $\gcd(a,b), d_a, d_b$ exist in $\mathfrak{o}_i$ and such that $\gcd(a,b), d_a, d_b$ keep their properties in extension $\mathfrak{o}_j$ all the way up to $\bar{\mathbb{Z}}$.
Does this strategy generally work out without refinements? Maybe. Let's look back at $S_2$. define $e_n=\prod_{i=1}^np_{i}^{n+1-i}$ where $p_i$ are ordered prime numbers. E.g. $e_1=2^1, e_2=2^2\cdot 3^1, e_3=2^3\cdot 3^2\cdot 5^1$ etc. Look at the following rings:
$$
S_{2,n}=\mathbb{Z}+\frac{1}{e_n}x\mathbb{Q}\left[\frac{1}{e_n}x\right].
$$
E.g. $35+\frac{1}{6}x+\frac{1}{216}x^2\in S_{2,3}\setminus S_{2,2}$: it is in $S_{2,3}$, because $6\mid 360$ and $216\mid 360^2$; it is not in $S_{2,2}$, because $216\nmid 12^2$.
We have $S_2=\bigcup_{n}S_{2,n}$. For $a=3,b=x$ we have $\gcd_{S_{2,1}}(3,x)=1$, but for every $n>1$ we have $\gcd_{S_{2,n}}(3,x)=3$ which matches $\gcd_{S_2}(3,x)=3$. Also $d_a=3,d_b=\frac{1}{3}x$ work out in $S_{2,n}$ for $n>1$ (but not for $n=1$).
If we wanted to generalize this tower strategy for general GCD domains $R$, one would have to show first that $R$ admits suitable towers. Currently I do not see how this should be possible without using anything about the specifics of the ring in question. It might also be necessary to make the tower depend on the given $a,b$. Anyway that tower strategy might open up an avenue for a proof that an GCD domain is gl or aid in the construction of a counterexample.
Note that Hutchins, "Examples of Commutative Rings", mentions in example 31 (p. 66) that any countable tower of PIDs is Bézout. This example deals with $S_4$ given in another notation.
Example: The ring $S_4=k[\mathbb{Q}_0^+]$ for a field $k$
According to an answer to this question any element of $S_4$ can be written as $cx^{q}$ for some $c\in S_4^{\times}$ and $q\in\mathbb{Q}_0^+$. For $a=c_ax^{q_a},b=c_bx^{q_b}\in S_4\setminus\{0\}$ assume wlog. $q_a\geq q_b$. We see that $S_4$ is a gl GCD domain by observing $\gcd(a,b)=b\lcm(a,b)=a$ and $d_a=a,d_b=1$.
Similarly one sees that $k[\mathbb{R}_0^+]$ is a gl GCD domain and one can think of a host of other examples with non-discrete semi-groups.
$\endgroup$
