Validity of Ampère's law in terms of $H$

If I understand correctly the question, the answer lies in the 6.3.3 paragraph of Griffiths' book. Yes, it does inherit the discontinuity. In this case, you could still use the differential equation just above and just below the discontinuity, and taking an integration path that crosses the boundary will give you boundary conditions for the magnetic field and for the auxiliary field. Another reference is this:

https://unlcms.unl.edu/cas/physics/tsymbal/teaching/EM-913/section6-Magnetostatics.pdf

Post-EDIT:

Yes, you can still use the theorem, if you are cautious with the domains of integration. The loop integral will be defined even in case of discontinuity because of what I said in the comments.

The contour integral across the boundary should not be considered as the direct application of Stoke's theorem in this area, but in a spirit similar to the principle of analytical continuation.

The Stokes' theorem

$$ oint_C mathbf{H}times dbf{l} = iint_A boldsymbolnabla times mathbf{H}cdot dmathbf{a} = mathbf{J}_{text{f,enclosed}}. $$

This relation is applicable when $boldsymbolnabla times bf{H}$ is existed.It is, of course, not hold when $boldsymbolnabla times bf{H}$ diverges.

In region across the boundary, even though $boldsymbolnabla times bf{H}$ diverges, the area integral failed, but the contour integral is still working. We then adopt the contour integral as the definition, extended to such regions.

Therefore, the relation

$$ oint_C bf{H} times dbf{l} = bf{J}_{text{f,enclosed}}. $$

is an analytical continuation into the across-boundary region, even the $nabla times bf{H}$ diverges therein.

Similarly, $boldsymbol{nabla}cdotbf{E}$ diverges at $r=0$ for $1/r^2$ field, but the surface integral still working there

$$ iint mathbf{E}cdot dbf{A}=frac{1}{epsilon_0} $$

We thus adopt the result from surface integral to define the strength of divergence.

$$ boldsymbol{nabla} cdot frac{hat r}{r^2} = -4pi delta(r). $$

Just use Ampère's law for $mathbf B$ to get the total current. Then subtract the bound current from the result to get the free current, if needed.

Is this integral equation valid if applied across a boundary of a magnetized material,(that is one end of the integral is inside the magnetised body and the other is out of it)? Or does it inherit any problems from its differential form?

There is an assumption regarding this boundary condition that is often overlooked for various reasons, but can sometimes become problematic. The limitations of this assumption are discussed in detail in Sections I.5 and I.6 of Classical Electrodynamics, Third Edition by John D. Jackson (i.e., blue cover version). So we start with: $$ oint_{C} mathbf{H} cdot dmathbf{l} = int_{S'} da left[ mathbf{j} + partial_{t} mathbf{D} right] cdot mathbf{n}' tag{0} $$ where $mathbf{H}$ is the magnetic field (technically, $mathbf{B}$ is the magnetic induction), $mathbf{D}$ is the electric displacement, $mathbf{j}$ is the current density (specifically the macroscopic average current density, see pages 248--258 in Jackson's book for definition and derivation), $S'$ is a closed surface with an outward unit normal $mathbf{n}'$, and $partial_{t} = tfrac{ partial }{ partial t }$.

Typically one sweeps the following under the rug, as it were. There can be a surface current density, $mathbf{K}$, that exists in a thin layer no thicker than one electron skin depth on the surface of the conducting material either caused by time-varying fields or merely present due to some source. In such scenarios, the right-hand side of Equation 0 changes to: $$ int_{S'} da left[ mathbf{j} + partial_{t} mathbf{D} right] cdot mathbf{t} = mathbf{K} cdot mathbf{t} Delta l tag{1} $$ where $mathbf{t}$ is the unit vector transverse to the surface $S'$ and $Delta l$ is the scale length of the pill box transverse to the surface $S'$.

There is another issue that arises at such a boundary involving what Jackson refers to as truly microscopic surface charge densities, i.e., $rholeft( x right) = sigma deltaleft( x right)$ where $sigma$ is the average surface charge density and $deltaleft( x right)$ is the Dirac delta function. The idealized scenario we are taught in class is that $sigma$ exists on the surface and that it has zero thickness, i.e., no charge density inside conductors, only on the surface. However, the truth is that $rholeft( x right)$ is confined to within $pm 2$ angstroms of the ''surface'' of the ionic distribution. For nearly all purposes, there is a discontinuity in the electric field at this boundary but in reality it likely varies over a finite length comparable to a few atomic widths or so.

However the curl is not defined on the boundary so is the loop integral well defined? Another way of saying is that: can I use stokes theorem when the curl of function isn't defined on some region?

It depends on what you are doing. Are you testing what Jackson would call the microscopic or the macroscopic? If the latter, things are much easier and the right-hand side of Equation 0 is can vanish if there is no significant surface current density, $mathbf{K}$, at the boundary. Even if present, one can still work with Equation 0 or 1, depending on scenario, and get meaningful results for macroscopic approximations.

Side Note
I did not rederive Jackson's definition of microscopic vs macroscopic because it's 10 pages in his book and not really necessary here. It basically involves noting the difference between spatial and temporal ansemble averages and why spatial are the correct choice then lots of details as to why XYZ is okay under WUV limits. The distinction in regards to this question is whether the OP wants to properly model $mathbf{H}$ across the boundary on scales down to the atomic or if they are okay with the typical, larger scale approximations of micrometers and up (give or take).