Determinant of arithmetic functions

Let $$f(a,b)=\sum_{d|a, \ d|b}d.$$ Let $M=(a_{i,j})_{i,j=1}^{n}$. Evaluate $\det M$.

Solution. The idea is utilization of Lower-upper triangular decomposition. Using a software (like Mathematica) one observes that for small $n$ the matrices $L$ and $U$ coulde be intepreted as - $L=(l_{i,j})_{i,j=1}^{n}$, where $$l_{i,j}=\begin{cases} 1, & \text{j divides i}\\
                           0, & \text{otherwise} \end{cases}$$ and $U=(u_{i,j})_{i,j=1}^{n}$ where $$u_{i,j}=\begin{cases} i, & \text{i divides j}\\                     
      0, & \text{otherwise} \end{cases}$$
It is directly verified that this holds in general (all $n$).
Obviously $\det L=1$, while $\det U=n!$ (the diagonal elements are $1,2,\ldots,n$). Thus $\det M=n!$.

The same idea could be employed for similar looking determinants, for example the one having $\gcd(i,j)$ at place $(i,j)$.