Code:

int  ConjGradientFiltered(float3N &X, const float3Nx3N &A, const float3N &B,const float3Nx3N &S){       // Solves for unknown X in equation AX=B       conjgrad_loopcount=0;       int n=B.count;       float3N q(n),d(n),tmp(n),r(n);       r = B - Mul(tmp,A,X);    // just use B if X known to be zero       filter(r,S);       d = r;       float s = dot(r,r);       float starget = s * sqr(conjgrad_epsilon);       while( s>starget && conjgrad_loopcount++ < conjgrad_looplimit)       {               Mul(q,A,d); // q = A*d;               filter(q,S);               float a = s/dot(d,q);               X = X + d*a;               if(conjgrad_loopcount%50==0)               {                       r = B - Mul(tmp,A,X);                       filter(r,S);               }               else               {                       r = r - q*a;               }               float s_prev = s;               s = dot(r,r);               d = r+d*(s/s_prev);               filter(d,S);       }       conjgrad_lasterror = s;       return conjgrad_loopcount<conjgrad_looplimit;  // true means we reached desired accuracy in given time - ie stable}