As long as DFT and the relevant exchange correlations are accepted, the converged results with respect to all the free parameters (Kmax, Ecut, relaxation, etc.) will be reasonable I guess. 

It is difficult to verify the correctness of concrete elastic moduli. Stiffness tensor however, must be positively-definite for stable configuration. So, calculate eigenvalues of obtained stiffness tensor

One thing you can do is to calculate the bulk modulus in two ways. One from the energy versus volume curve and the other from the appropriate expression using the elastic constants (for cubic B = 1/3(c11+2C12)). If they agree chances are good your

elastic constants are good to. This is no guarantee, of course.

Many thanks for the helpful suggestions .

You can use VASP to calculate the elastic moduli with IBRION=5 or 6， see the vasp manual.

There is no real way of knowing whether a calculation is accurate without using experimental data as a basis for checking. Of course this is not always possible, in which case comparing with other calculations is a possible method. If you are the first to do a calculation and there is no experimental data, then it would be worth doing experiments or alternatively see if there are theoretical models for the calculation you've done using ab initio methods. Calculating using another method, e.g. molecular dynamics, is possible. Really models are there to model reality, which is why experimental data is used as a basis. Note you are calculating at T = 0K, or a slightly higher temperature if you are using Fermi smearing and a set of assumptions are made in DFT. This means the results may not match 100% with experimental data (often conducted at room temperature along a fixed direction or plane).

In computational calculations it is always suggested to go for benchmark calculations.  First of all you have to know weather correct parameters are used in INCAR file to get proper elastic modulli. I suggest you to take a system which is already has reported using other codes (not using VASP) and otain the elastic constants. Then check your VASP values and get a consistency (by varying different parameters or optimization). If you already have done this, you might have followed the proper method to obtain the elastic constants for your unknown system. Besides, obtain these values with different exchange functional like GGA (PBE & PW91) and LDA. If your results are consistent in mechanical stability, Bulk modulus, etc  with using different functional then you can relay on your results. Note: There will be little different values from different functional but the trends should be same.

I would contend with the last answer. Just because a set of answers seem similar using different exchange-correlation functionals, doesn't mean the result is "correct". It just means the result is converged within the applied physical model. You may still be "off" relative to what is observed in nature. As Richard Feynmann put it, "It doesn't matter how beautiful your theory is, it doesn't matter how smart you are. If it doesn't agree with experiment, it's wrong.". It really is a tough tricky question to answer. Some compare results to all-electron calculations too, something I forgot to mention above. I don't think you can get more than a ball park figure using DFT, though the stability of the calculation is really a function of the material you are using. Good luck!

Agreement with the experiments is the ultimate evidence for proper theoretical calculations. Since, the experimental results are not available to compare it is always better to check the obtained values with other theoretical methods. Experimental agreement always does not mean to correctness of theory, rather a complement for both the experiments and theory to bring out the physics/chemistry of the materials. I feel one no need to wait for theory or experiments to be done for consistency, let your work know the scientific world then that may motivate to initiate theory/experiments.

Agreed! Though other xc functionals is perhaps a weak form of "other theoretical methods" (limited to the accuracy of the DFT method still). I suggested alternative theoretical methods too, perhaps an analytical method or modelling in a different way, e.g. molecular dynamics or all-electron calculations.

The standard calculation of the elastic constants is by changing the unit cell dimensions or angles and look at the energy change.

The other possibility is through the interatomic force constants (second derivatives of the energy wrt to atomic displacements, from which phonons are calculated).

There is a formula which gives the elastic constants as a function of the force constants:

C=sum_R phi(R ) R_i R_j where phi's are the force constants, R is a nearest neighbor vector and i,j=x,y,z

You'd have to lookup the exact formula. A good reference is "Thermodynamics of crystals" by Wallace or "Dynamical theory of crystal lattices" by Born and Huang.

 This igives you a consistency check. Although the total energy calculations are rather fast, but using the force constants gives access to all the elastic constants AND phonons in one shot.