It seems like the sum $S(n)$ should be possible to upperbound by an expression of the form ${\mathcal O}(n^a\cdot \log^b(n))$ as $n \rightarrow \infty$:
$$
S(n)\stackrel{\triangle}{=}\sum\_{1\leq x \neq y\leq n} \frac{gcd(x,y)^2}{x y}.
$$
Any pointers, ideas, appreciated. Since $gcd(x,y) \leq \min(x,y)$ an upper bound is given by
$$
S(n) \leq \sum\_{1\leq x \neq y\leq n} \frac{\min(x,y)^2}{x y} = 2 \sum\_{1\leq x < y\leq n} \frac{x}{y}
$$
which can be rewritten as
$$
S(n) \leq 2 \sum\_{1\leq x \leq n-1} x \sum\_{x+1 \leq y \leq n} \frac{1}{y} = 
2 \sum\_{1\leq x \leq n-1} x (H(n)-H(x)),
$$
where $H(n)$ denotes the harmonic sum. This then gives
$$
S(n) \leq 2 H(n) \left[\sum\_{1\leq x \leq n-1} x\right] -2 \left[\sum\_{1\leq x \leq n-1} x H(x)\right]
$$
or
$$
S(n) \leq n(n-1) H(n) - 2 (1\times H(1)+2 \times H(2)+\cdots+(n-1) \times H(n-1))
$$
or to a first order approximation
$$
S(n) \leq n(n-1) H(n) - 2 \sum\_{1\leq x \leq n-1} x \log x \approx n^2 \log n - 2 (1+\sum\_{2\leq x \leq n-1} x \log x) .
$$
An application of the arithmetic geometric mean inequality an Stirling seems to give an estimate of $(n-1)^2$ for the subtracted sum while an integral approximation says that the subtracted sum is ${\cal O}(n^2 \log n)$ which would cancel the first term! Using the estimate $\log n + \gamma$ for the Harmonic sum then leaves us with $S(n) \leq \gamma n(n-1).$


Is there a more accurate estimate?



An upper bound to sum of ratios of gcd's and products
Suppose that $x <y$. Write $x=ga$ and $y=gb$ where $g$ is the gcd of $x$ and $y$; thus $a$ and $b$ are coprime with $a<b$, and $g\le n/b$. So 
$$ 
S(n) = 2 \sum\_{\substack{{a<b \le n} \\ {(a,b)=1}}} \frac{1}{ab} \sum\_{g\le n/b} 1 \le 2\sum\_{\substack{{a<b\le n} \\ {(a,b)=1}}} \frac{1}{ab} \frac{n}{b}.
$$
Now ignore the condition that $(a,b)=1$, getting (using $\sum\_{b>a} 1/b^2 \le 1/a$)
$$ 
S(n) \le 2n \sum\_{a\le n} \frac{1}{a} \sum\_{a<b\le n} \frac{1}{b^2} \le 2n\sum\_{a\le n} \frac{1}{a^2} \le \frac{\pi^2}{3} n. 
$$
With more effort you can prove that $S(n) \sim Cn$ for some constant $C$. 


Edit: The OP Kodlu noted that numerically the constant $C$ seemed to be $2$, which was surprising to me. So here's a sketch as to why the constant comes out to be $2$. The argument above is pretty accurate, and is essentially an asymptotic. Thus from above one finds that 
$$ 
C = 2\sum\_{\substack{{1\le a<b } \\ {(a,b)=1} }} \frac{1}{ab^2}. 
$$ 
Now use Mobius inversion $\sum\_{d|(a,b)} \mu(d) =1$ if $(a,b)=1$ and zero otherwise. Thus, with $a=dr$ and $b=ds$,  

$$ 
C= 2\sum\_{1\le a<b} \frac{1}{ab^2} \sum\_{d|(a,b)} \mu(d) =2 \sum\_{d=1}^{\infty} \frac{\mu(d)}{d^3} \sum\_{1\le r<s} \frac{1}{rs^2}. 
$$ 
The sum over $d$ is $1/\zeta(3)$, and the sum over $r$ and $s$ is the multizeta value $\zeta(2,1)$ which Euler showed equals $\zeta(3)$. Thus $C=2$ and $S(n) \sim 2n$ (and in fact the argument above gives this as an upper bound always, as seen in the numerics). 


