#pragma once
#include "getmatrix.h"
#include <vector>

using namespace std;

/**     (A) Inner product of two vectors (from skalar_stl)
	@param[in] x	vector
	@param[in] y	vector
	@return 	    resulting Euclidian inner product <x,y>
*/
double scalar(vector<double> const &x, vector<double> const &y);

/**     (A) 5.(b) Inner product of two vectors using the Kahan scalar product
	@param[in] x	vector
	@param[in] y	vector
	@return 	    resulting Euclidian inner product <x,y>
*/
double Kahan_skalar(vector<double> const &x,  vector<double> const &y);

/** 	(A) 6. cBLAS scalar product of two vectors
	@param[in] x	vector
	@param[in] y	vector
	@return 	    resulting Euclidian inner product <x,y>
*/
double scalar_cBLAS(vector<double> const &x, vector<double> const &y);


/** 	(B) Matrix-vector product (from intro_vector_densematrix)
 * 	@param[in] A	dense matrix (1D access)
 *  @param[in] u	vector
 *
 *	@return    resulting vector
*/
vector<double> MatVec(vector<double> const &A, vector<double> const &x);

/** 	(B) 6. cBLAS Matrix-vector product
 * 	@param[in] A	dense matrix (1D access)
 *  @param[in] u	vector
 *
 *	@return    resulting vector
*/
vector<double> MatVec_cBLAS(vector<double> const &A, vector<double> const &x);


/** 	(C) Matrix-matrix product
 * 	@param[in] A			MxL dense matrix (1D access)
 *  @param[in] B			LxN dense matrix (1D access)
 *  @param[in] shared_dim 	shared dimension L
 *
 *	@return    resulting MxN matrix
*/
vector<double> MatMat(vector<double> const &A, vector<double> const &B, size_t const &shared_dim);

/** 	(C) 6. cBLAS Matrix-matrix product
 * 	@param[in] A			MxL dense matrix (1D access)
 *  @param[in] B			LxN dense matrix (1D access)
 *  @param[in] shared_dim 	shared dimension L
 *
 *	@return    resulting MxN matrix
*/
vector<double> MatMat_cBLAS(vector<double> const &A, vector<double> const &B, size_t const &shared_dim);


/** 	(D) Evaluation of a polynomial function using Horner's scheme
 * 	@param[in] a	coefficient vector
 *  @param[in] x	vector with input values
 *
 *	@return    vector with output values
*/
vector<double> poly(vector<double> const &a, vector<double> const &x);


/**     (E) Solves linear system of equations  K @p u = @p f  via the Jacobi iteration (from jaboci_oo_stl)
 * We use a distributed symmetric  CSR matrix @p SK and initial guess of the
 * solution is set to 0.
 * @param[in] SK	CSR matrix
 * @param[in] f		distributed local vector storing the right hand side
 * @param[out] u	accumulated local vector storing the solution.
*/
void JacobiSolve(CRS_Matrix const &SK, vector<double> const &f, vector<double> &u);