ok I guess I found you're problem:

matrix.m[8] = -forward.x;
matrix.m[9] = -forward.y;
matrix.m[10] = -forward.z;

you have the forward in the opposite direction (this is the right way I do it like this too but references and actual related math don't do like this) so you need to negate the z value of the starting point after calculating it

code will be like:

TRay RayFromScreen(int x, int y, int sw, int sh, float fov, float z_near, float z_far, float aspect)
{

	vec2 NormalizedCoordinates;
	NormalizedCoordinates.x = (float(x) / float(sw) ) * 2.0 - 1.0;
	NormalizedCoordinates.y = (float(y) / float(sh) ) * 2.0 - 1.0;

	vec4 NDC=vec4(NormalizedCoordinates.x, NormalizedCoordinates.y,-1.0, 1.0); //z=-1 because its on near plane

	TRay res;
	mat4 ProjectionMatrix = CAM_PROJECTION;
	ProjectionMatrix.Transpose();
	mat4 ViewMatrix = CAM_VIEW;
	ViewMatrix.Transpose();
	mat4 pvm = ProjectionMatrix*ViewMatrix;
	pvm.Inverse();
	res.start=pvm*NDC;
	res.start.z=res.start.z*-1.0f;

	vec3 RayDirection=res.start-FPP_CAM->pos;
	RayDirection=(RayDirection/RayDirection.z)*(z_far-z_near);//normalize by Z resize the z by farZ-nearZ so that it hits the far plane

	res.end=res.start+RayDirection;

return 

gluperspective function accepts fov as fovy not fov for xz plane

I'm not getting your point, maybe chat'd be better

for gluPerspecitve you define a fov angle which is between top and bottom clipplanes,

yes I understand that, but where are you getting with that? it's supposed to be like that in other functions too, you have fovy, D3D does the same too

proper results....