Given three coordinates (a,b,c), (d,e,f), and (l,m,n), what is the center of the circle in the 3D plane (h,k,i) that contains these three points. I have tried the following:
$$(a-h)^2+(b-k)^2+(c-i)^2=r^2$$
$$(d-h)^2+(e-k)^2+(f-i)^2=r^2$$
$$(l-h)^2+(m-k)^2+(n-i)^2=r^2$$
Subtracted equation 2 from 1, equation 3 from equation 2, and equation 3 from equation 1 to get:
$$2(a-d)h+2(b-e)k+2(c-f)i=a^2-d^2+b^2-e^2+c^2-f^2$$
$$2(d-l)h+2(e-m)k+2(f-n)i=d^2-l^2+e^2-m^2+f^2-n^2$$
$$2(a-l)h+2(b-m)k+2(c-n)i=a^2-l^2+b^2-m^2+c^2-n^2$$
Closer inspection revealed to me that the left side formed a 3 x 3 matrix and the left a 3 by 1 matrix.
Tried solving for $(h,k,i)$ by taking the inverse of the left side matrix and multiplying that by the right side 3 by 1 matrix, but it turns out the left side matrix's rows are dependent, so the determinant of the left side matrix is 0 (singular matrix); hence, no inverse.
Please guide me as to solve this problem.
Thank you in advance.
 A: The problem is that your equations don't capture the totality of conditions.
In fact, at present, you have three unknowns but, as you have remarked, only two independent equations.
This situation reflects that you have one degree of freedom. Indeed, what you have expressed is constraints on spheres, and there is an infinite number of spheres passing through 3 points in 3D space.
The supplementary condition you need to add is the coplanarity of the points and the center of the circle; this coplanarity is given by the nullity of the determinant:
$$\left|\begin{matrix}a&d&l&h\\b&e&m&k\\c&f&n&i\\1&1&1&1\end{matrix}\right|=0$$
(see http://mathworld.wolfram.com/Coplanar.html).
Solving this system of equations in 2016 does not rely on hand calculation but on the  use of  a computer algebra system (Maple, Mathematica...).
A: Maybe a faster way is to use Euler's theorem, $H=3G+2O$.
If we have to find the circumcenter of $ABC$, we may compute the projection of $C$ on the $AB$-line and the projection of $B$ on the $AC$-line through scalar products, or by finding the minimum of a squared distance (they are the same approach, indeed). So the orthocenter $H$ can be found by intersecting two lines. 
The centroid is very easy to compute by just considering the arithmetic mean of the coordinates of the vertices. So $O$ can be found through $O=\frac{H-3G}{2}$.

An alternative and faster approach is to exploit the fact that the barycentric coordinates of the circumcenter (i.e. the coefficients
  $(u,v,w)$ such that $O=uA+vB+wC$) just depends on the squared distances
  $AB^2,AC^2,BC^2$, that also in the 3d-case are straightforward to compute through the pythagorean theorem.

A: For points $(a,b,c)$ and $(d,e,f)$, their perpendicular bisector can be found by:
$$\begin{align*}
(x-a)^2 + (y-b)^2+(z-c)^2 &= (x-d)^2 + (y-e)^2 + (z - f)^2\\
a^2-2ax+b^2-2by+c^2 -2cz &= d^2 - 2dx + e^2 - 2ey + f^2 - 2fz\\
2(a-d)x +2(b-e)y + 2(c-f)z &= a^2+b^2+c^2-d^2-e^2-f^2
\end{align*}$$
Do the same and find the perpendicular bisector of $(d,e,f)$ and $(l,m,n)$.
For three points $(a,b,c)$, $(d,e,f)$ and $(l,m,n)$, their common plane $px+qy+r = z$ can be found by:
$$\begin{pmatrix}
a&b&1\\d&e&1\\l&m&1
\end{pmatrix}\begin{pmatrix}
p\\q\\r
\end{pmatrix} = \begin{pmatrix}
c\\f\\n
\end{pmatrix}$$
This might have some problem if their common plane is perpendicular with the $xy$-plane, but you can always change to use $x$ or $y$ to be the subject of the equation and try to solve again. You are using a computer anyway.
Now you have three equations in the form of $tx + uy + vz = w$, and you should be able to solve for their intersection $(x,y,z) = (h,k,i)$.
