Why is $ \nabla \cdot \vec{E} \neq 4 \pi\ k$ times the charge density? For a line charge with density $\lambda$, we wish to compute the divergence at a point on line charge. We start from definition of divergence:
$ \nabla \cdot \vec{E}=\dfrac{d \iint \vec{E} \cdot d\vec{S}}{dV}
=4 \pi\ k\ \dfrac{q_{\text{ enclosed}}}{dV}=4 \pi\ k\ \dfrac{dq}{dV}$
But this is not equal to $4 \pi\ k\ \lambda = 4 \pi\ k\ \dfrac{dq}{dl}$, i.e. this expression does not contain line charge density.
But my book says $\nabla \cdot \vec{E}$ at source points is always $4 \pi\ k$ times charge density.
 A: Like Phillip Wood mentioned in his answer, a line charge would have infinite charge density along the line and zero outside it. However there is a way around this problem. What we can do is to introduce delta function (Dirac delta to be precise) $\delta(x)$ to solve this.
Basically a delta function,according to Physicists, is defined as 
$$ \delta(x) = \begin{cases}
    \infty,& \text{if } x = 0\\
    0,              & \text{otherwise}
\end{cases}
$$
Also ,
$$
\int_{-\epsilon}^{+\epsilon}{\delta(x)\mathop{}\!\mathrm{d}x} = 1, \qquad \epsilon > 0
$$
Using this we can define total charge density for line charge along $z$-axis as follows 
$$\rho(x,y,z) = \lambda \delta(x)\delta(y)$$
This is consistent with the observation that charge density is infinite along the line (i.e., $x=y=0$) and zero everywhere else.
To get line charge density along $z$-axis we integrate $\rho$ through all space at each $z$
$$
\text{Line charge density} = \int_{x}\int_{y}\rho(x,y,z)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}y=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\lambda \delta(x)\delta(y) \mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}y
$$
Use the definition of delta function
$$
\text{Line charge density} = \lambda
$$
This exactly what we expected. 
Now you can try to integrate $\rho$ to get total charge (Ans: $\int\rho \mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}y\mathop{}\!\mathrm{d}z = l \lambda$ where, $l$ is length of line element).
Thus the divergence can written as 
$$
 \nabla\cdot \mathbf{E} = 4\pi k \rho = 4\pi k \lambda \delta(x) \delta(y)
$$
This what you wanted in the question.
To get the integral form 
$$
\int \nabla\cdot\mathbf{E}\mathop{}\!\mathrm{d}V = \int 4\pi k \rho(x,y,z) \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}y \mathop{}\!\mathrm{d}z
\\
\implies \oint \mathbf{E}\cdot\mathop{}\!\mathrm{d}\mathbf{S} = 4 \pi k \lambda l
$$
Note:
The delta function is extremely useful. Imagine you have a point particle with charge $q$ at origin.
It will have charge density
$$
\rho(x,y,z) = q \delta(x)\delta(y)\delta(z)
$$
Thus, Gauss' law will be
$$
 \nabla\cdot\mathbf{E} = 4\pi k \rho = 4\pi k q \delta(x) \delta(y) \delta(y)
$$
Integrating this (using property of delta function) gives 
$$
\oint \mathbf{E}\cdot\mathop{}\!\mathrm{d}\mathbf{S} = 4 \pi k q
$$
Now let the surface be surface of sphere with radius $R$. Now, by spherical symmetry (the problem has no preferred direction: if you rotate a point charge it will still look the same) electric field will be function of distance from charge, $r$ only. Also because electric force is a central force, electric field will be along the radial direction (i.e., $\mathbf{E} = |E(r)| \hat r $).Thus
$$
\oint \mathbf{E}\cdot\mathop{}\!\mathrm{d}\mathbf{S} = \oint |E(r)|\mathop{}\!\mathrm{d}A = |E| \oint_{sphere} \mathop{}\!\mathrm{d}A = 4\pi R^{2}|E| \\
\implies 4\pi R^{2}|E| = 4 \pi k q \\
\implies \mathbf{E} = \frac{k q}{R^{2}} \hat{r}
$$
Just what you expected from a point charge!! 
Even though there are simpler ways to see the obvious result above, it becomes extremely useful in complicated problems and is applied almost everywhere in Physics. Try exploring more about this!
A: A line charge is a convenient fiction: the line (which has no thickness) is supposed to be charged uniformly all along its length. Because the line has no thickness, the charge density, $\rho,$ that is the $\frac{dq}{dV}$ that appears in the version of Gauss's law that you have quoted, is infinite all along the line. To side-step the difficulty we define $\lambda,$ the linear charge density, the charge per unit length of line. This will be finite.
This doesn't help you to apply Gauss's law to a point actually on the line. $\rho$ is infinite, so the divergence is infinite. But you can apply Gauss's law in the form you quoted to all points outside the line itself, because $\rho$ is zero at all these points. 
More usefully, you can apply the integral from of Gauss's law to an imaginary cylindrical surface of radius r with the line of charges along its axis. This is because you know that for a cylinder of axial length $l,$ the charge contained is $\lambda l.$ Thus (with your notation for the constant)
 $$\iint \textbf{E}\cdot d \textbf {S} = 4\pi k \lambda l.$$ 
Because of the symmetry, $\textbf{E}$ is normal to the cylindrical surface and of the same magnitude all along and all round the surface, so the surface integral evaluates to $4\pi r^2 lE$.  
