I decided to post a video of what i get now:

You will see there that sometimes ray isn't computed properly (most of times when reaching corners of the screen)

And full code:

TRay RayFromScreen(int x, int y, int sw, int sh, float fov, float z_near, float z_far, float aspect)
{
mat4 mvm = CAM_MODEL * CAM_VIEW;
mvm.Transpose(); <- i have a code that gets up and right vector from oglmatrix modelview so i just transpose it here
vec3 dirX, dirY;
	dirX.x = mvm.m[0];
	dirX.y = mvm.m[4];
	dirX.z = mvm.m[8];
//
//
	dirY.x =	mvm.m[1];
	dirY.y =	mvm.m[5];
	dirY.z =	mvm.m[9];


	TRay res;
//	cot(x) = 1/tan(x)
//	x = z_near / ctg(a);

	float a = fov / 2.0;
float cotangent = 1.0 / tan( a * imopi );

float ax = z_near / cotangent;

float screen_w = 2.0*ax;

float yratio = 1.0 / aspect;

float screen_h = screen_w * yratio;


float scr_coord_x = float(x) / float(sw);
float scr_coord_y = float(sh - y) / float(sh);

vec3 dir = FPP_CAM->ReturnFrontVector();


//move to lower left corner
vec3 start_pos = (FPP_CAM->pos + dir * z_near) + (-dirX * (screen_w / 2.0)) + (-dirY * (screen_h/2.0));

res.start = start_pos + (dirX * (screen_w * scr_coord_x)) + (dirY * (screen_h * scr_coord_y));


//compute the world position on the other end since its a perspective projection
ax = z_far / cotangent;
screen_w = 2.0*ax;
screen_h = screen_w * yratio;

start_pos = (FPP_CAM->pos + dir * z_far) + (-dirX * (screen_w / 2.0)) + (-dirY * (screen_h/2.0));

res.end = start_pos + (dirX * (screen_w * scr_coord_x)) + (dirY * (screen_h * scr_coord_y));



return res;
}


//this function calls the frist one
t3dpoint<int> CellFromScreenCPU(int x, int y, int sw, int sh, float fov, float z_near, float z_far, float aspect)
		{
	TRay ray = RayFromScreen(x, y, sw, sh, fov, z_near, z_far, aspect); 
}




actual call:
t3dpoint<int> cell =	CellFromScreenCPU(x,y, SCREEN_WIDTH, SCREEN_HEIGHT, 90.0, 2.0, 28.0 * 1000.0, float(SCREEN_WIDTH) / float(SCREEN_HEIGHT));

lets say you have clicked on a point at (x,y) first you need to normalize it, which will be (x/width, y/height) next thing will be no move this to NDC coordinates, this point is considered on the near plane so what you do is:

(you need farZ, nearZ, camera position and projection and view matrices)

(I have used this method it works perfect for me)

vec2 NormalizedCoordinates=ScreenCoordinates/ScreenSize;
vec4 NDC=vec4( (NormalizedCoordinates-0.5)*2,-1,1); //z=-1 because its on near plane
 
TRay res;
 
res.start=vec3(inverse(ProjectionMatrix*ViewMatrix)*NDC);//this move our point from the NDC coordinates to world coordinates
 
vec3 RayDirection=vec3(res.start-CameraPosition);//ray direction is from camera to ray's starting point
RayDirection=RayDirection/RayDirection.z*(farZ-nearZ);//normalize by Z resize the z by farZ-nearZ so that it hits the far plane
 
res.end=res.start+RayDirection;

nah it doesnt work btw i think NormalizedCoordinates - 0.5 is wrong since to get ndc coord you need to do * 2.0 - 1.0

any ways function always shoots ray from the screen center to the 0,0,0 point...(or to be more precise from 0,0,0 point to direction(from 0,0,0 to cam_pos)*z_far... no matter what i set to normalized coords

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;

	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 res;
}

about the NDC cords "vec4( (NormalizedCoordinates-0.5)*2,-1,1);" is the same as vec4( NormalizedCoordinates*2-1.0,-1,1);.

About the function, it seems that your doing the exact thing I said, and if you agree the function mathematically should be working, one thing I suspect is that maybe your projection matrix has a different behavior, esp. on mapping z, i'd say take a look at the res.start and check if it's calculated correctly or if its w is 1, your matrix may result in other values for w so yes that's what I suspect.

the code i wrote is working good but i had to disable the apsect ration for y screen, looks like fov is for screen HEIGHT not for WIDTH, i am confused

and aspect is width / height

	float a = fov / 2.0;
float cotangent = 1.0 / tan( a * imopi );

float ax = z_near / cotangent;

float screen_w = 2.0*ax;

float screen_h = screen_w;// * yratio; <--- now i dont multiple by 1.0/aspect

screen_w = screen_w * aspect;  <--- here

float scr_coord_x = float(x) / float(sw);
float scr_coord_y = float(sh - y) / float(sh);

well it depend what the fov is for, it mostly for Y, depend on your projection matrix function

i did something liek this:

		float sx = x;
		float sy = y;
		float width = sw;
		float height= sh;

	    float d = 1.0f/tan(fov*pi/360.0f);
	  //  float aspect = width/height;
	    vec3 viewPoint = vec3( 2.0f*aspect*sx/width - aspect, -2.0f*sy/height + 1.0f, -d);//// : d );
mat4 invV = CAM_VIEW;
invV.Inverse();
	   vec3 viewPointI = invV * viewPoint; //i assume thats viewToWorldMatrix

	   res.start 	= FPP_CAM->pos;
	   res.end 		= res.start + Normalize(viewPointI-FPP_CAM->pos)*4000.0;
return res;

that really doesnt work :/

HELP!

nah it doesnt work btw i think NormalizedCoordinates - 0.5 is wrong since to get ndc coord you need to do * 2.0 - 1.0

any ways function always shoots ray from the screen center to the 0,0,0 point...(or to be more precise from 0,0,0 point to direction(from 0,0,0 to cam_pos)*z_far... no matter what i set to normalized coords

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;

	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 res;
}

the function is correct here, it should be working, though your new function is not, maybe send the code that makes the projection and view matrices, that may help, also I really am not getting what the problem is, maybe some videos or pics help

ray isnt hitting desired part of screen

void Translate(T x, T y, T z)
{
	Matrix44<T> tran( 1, 0, 0, x,
					  0, 1, 0, y,
					  0, 0, 1, z,
					  0, 0, 0, 1);
	(*this) = (*this) * tran;
}

template <class T> void
glLookAt(Matrix44<T> &matrix, t3dpoint<T> eyePosition3D, t3dpoint<T> center3D, t3dpoint<T> upVector3D )
{
   t3dpoint<T>  forward, side, up;
   forward = Normalize( vectorAB(eyePosition3D, center3D) );
   side = Normalize( forward * upVector3D );
   up = side * forward;
  matrix.LoadIdentity();

	matrix.m[0] = side.x;
	matrix.m[1] = side.y;
	matrix.m[2] = side.z;

	matrix.m[4] = up.x;
	matrix.m[5] = up.y;
	matrix.m[6] = up.z;

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



Matrix44<T> transgender;
transgender.Translate(-eyePosition3D.x, -eyePosition3D.y, -eyePosition3D.z);


matrix = transgender * matrix;
}



template <class T>
void gluPerspectiveA(Matrix44<T> & matrix, T fovy, T aspect, T zmin, T zmax)
{
 T xmin, xmax, ymin, ymax;
 ymax = zmin * tan(fovy * M_PI / 360.0);
 ymin = -ymax;
 xmin = ymin * aspect;
 xmax = ymax * aspect;
 //mglFrustum(m, xmin, xmax, ymin, ymax, zmin, zmax);
 //mglFrustum(Matrix44<T> & matrix, T l, T r, T b, T t, T n, T f)


matrix.m[0] = 	(2.0*zmin)/(xmax-xmin);
matrix.m[1] = 	0.0;
matrix.m[2] =  	(xmax + xmin) / (xmax - xmin);
matrix.m[3] =	0.0;

matrix.m[4] = 	0.0;
matrix.m[5] = 	(2.0*zmin) / (ymax - ymin);
matrix.m[6] =   (ymax + ymin) / (ymax - ymin);
matrix.m[7] =   0.0;

matrix.m[8] = 	0.0;
matrix.m[9] =  	0.0;
matrix.m[10] =	-(zmax + zmin) / (zmax-zmin);
matrix.m[11] =  (-2.0*zmax*zmin) / (zmax-zmin);

matrix.m[12] =  0.0;
matrix.m[13] =  0.0;
matrix.m[14] =  -1.0;
matrix.m[15] =  0.0;


}