But if those c's are independent of your unknown c's you may want to know if there exists a unique and exact solution for your problem.

In this case define a modified (augmented) matrix

mm2 = Transpose[Join[Transpose[mat], {vec}]];
MatrixForm[mm2]

and perform a Gauss-Elimination (sort of)

mnn = mm2;
Do[
 x = mnn[[js, js]];
 If[ x == 0,
  jz = js;
  While[
   And[jz <= 21, x == 0],
   jz = jz + 1;
   x = mnn[[jz, js]]];
  If[x =!= 0,
   r = mnn[[js]];
   mnn[[js]] = mnn[[jz]];
   mnn[[jz]] = r]];
 Do[
  If[jz != js,
   x = mnn[[js, js]];
   If[x =!= 0,
    mnn[[jz]] = mnn[[jz]] - mnn[[jz, js]]/x mnn[[js]]]],
  {jz, 1, 21}],
 {js, 1, 9}]
MatrixForm[mnn];

crit = Drop[Transpose[mnn][[-1]], 9] // FullSimplify

Then your problem has a solution if and only if all components of crit are zero.

Your system has more equations than unknowns. Solve only gives generic solutions. You may try Reduce, but it is probably too complicated. Here I replace the coefficients with random numbers and check if there are solutions:

eqs = {uj == aj*c1 + bj*c2 + cj*c3 + ui, 
   vj == aj*c4 + bj*c5 + cj*c6 + vi, tj == aj*c7 + bj*c8 + cj*c9 + ti,
    uk == ak*c1 + bk*c2 + ck*c3 + ui, 
   vk == ak*c4 + bk*c5 + ck*c6 + vi, tk == ak*c7 + bk*c8 + ck*c9 + ti,
    ul == al*c1 + bl*c2 + cl*c3 + ui, 
   vl == al*c4 + bl*c5 + cl*c6 + vi, tl == al*c7 + bl*c8 + cl*c9 + ti,
    um == am*c1 + bm*c2 + cm*c3 + ui, 
   vm == am*c4 + bm*c5 + cm*c6 + vi, tm == am*c7 + bm*c8 + cm*c9 + ti,
    un == an*c1 + bn*c2 + cn*c3 + ui, 
   vn == an*c4 + bn*c5 + cn*c6 + vi, tn == an*c7 + bn*c8 + cn*c9 + ti,
    uo == ao*c1 + bo*c2 + co*c3 + ui, 
   vo == ao*c4 + bo*c5 + co*c6 + vi, to == ao*c7 + bo*c8 + co*c9 + ti,
    up == ap*c1 + bp*c2 + cp*c3 + ui, 
   vp == ap*c4 + bp*c5 + cp*c6 + vi, tp == ap*c7 + bp*c8 + cp*c9 + ti};
vars = {c1, c2, c3, c4, c5, c6, c7, c8, c9};
coeffs = Complement[Cases[eqs, _Symbol, All], vars];
values = RandomInteger[{1, 2}, Length[coeffs]]
eqs /. Thread[coeffs -> values]
Reduce[%]