docs: public static Vec solve(double eps, Matrix A, Vec x, Vec b, Matrix Minv)
      {
          if(!A.isSquare() || !Minv.isSquare())
              throw new ArithmeticException("A and Minv must be 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");
          else if(A.rows() != Minv.rows() || A.cols() != Minv.cols())
              throw new ArithmeticException("Matrix A and Minv do not have the same dimmentions");
          
          int k = 0;
          Vec r_k = b.subtract(A.multiply(x));
          Vec z_k = Minv.multiply(r_k);
          Vec p_k = z_k.clone();
          Vec Apk;
          double rkzk = r_k.dot(z_k);
          
          do
          {
              Apk = A.multiply(p_k);
              
              double alpha = rkzk/p_k.dot(Apk);
              x.mutableAdd(alpha, p_k);
              r_k.mutableSubtract(alpha, Apk);
              
              if(r_k.dot(r_k) < eps*eps)
                  return x;
              
              z_k = Minv.multiply(r_k);//TODO implement method so we can reuse z_k space, instead of creating a new vector
              
              double newRkZk = r_k.dot(z_k);
              double beta = newRkZk/rkzk;
              rkzk = newRkZk;
              
              p_k.mutableMultiply(beta);
              p_k.mutableAdd(z_k);
              
              
          }
          while(k++ < A.rows());
          
          return x;
      }