docs: public static Vec solve(double eps, Matrix A, Vec x, Vec b)
      {
          if(!A.isSquare())
              throw new ArithmeticException("A must be a square (symmetric & positive definite) matrix");
          else if(A.rows() != b.length() || A.rows() != x.length())
              throw new ArithmeticException("Matrix A dimensions do not agree with x and b");
          int k = 0;
          Vec r_k = b.subtract(A.multiply(x));
          Vec p_k = r_k.clone();
          Vec Apk;
          
          double RdR = r_k.dot(r_k);
          do
          {
              Apk = A.multiply(p_k);
              double alpha_k = RdR /  p_k.dot(Apk) ;
              
              x.mutableAdd(alpha_k, p_k);
              
              r_k.mutableAdd(-alpha_k, Apk);
              
              double newRdR = r_k.dot(r_k);
              
              //Stop when we are close enough
              if(newRdR < eps*eps)
                  return x;
              
              double beta_k = newRdR/RdR;
              
              p_k.mutableMultiply(beta_k);
              p_k.mutableAdd(r_k);
              
              //Set up for next ste
              RdR = newRdR;
          }
          while(k++ < A.rows());
          
          return x;
      }