/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "gctl/core.h"
#include "gctl/algorithm.h"
#include "lcg/lcg.h"
#include "ctime"

#define M 100
#define N 80

// 普通二维数组做核矩阵
gctl::matrix<double> kernel(M, N, 0.0);
// 稀疏矩阵为核矩阵
gctl::spmat<double> sp_kernel(M, N, 0.0);
// 中间结果数组
gctl::array<double> tmp_arr(M, 0.0);

// 计算核矩阵乘向量的乘积
void CalAx(void* instance, const lcg_float* x, lcg_float* prod_Ax, const int n_s)
{
	for (int i = 0; i < M; i++)
	{
		tmp_arr[i] = 0.0;
		for (int j = 0; j < n_s; j++)
		{
			tmp_arr[i] += kernel[i][j] * x[j];
		}
	}

	for (int j = 0; j < n_s; j++)
	{
		prod_Ax[j] = 0.0;
		for (int i = 0; i < M; i++)
		{
			prod_Ax[j] += kernel[i][j] * tmp_arr[i];
		}
	}
	return;
}

// 计算核矩阵乘向量的乘积
void CalAx_Spmat(void* instance, const lcg_float* x, lcg_float* prod_Ax, const int n_s)
{
	// 直接调用稀疏矩阵与向量的乘法
	// 注意第二次为向量乘矩阵 相当于矩阵的转置与向量相乘
	sp_kernel.multiply_vector(x, n_s, tmp_arr.get(), M);
	sp_kernel.multiply_vector(tmp_arr.get(), M, prod_Ax, n_s, gctl::Trans);
	return;
}

//定义共轭梯度监控函数
int Prog(void* instance, const lcg_float* m, const lcg_float converge, const lcg_para* param, const int n_s, const int k)
{
	std::clog << "Iteration-times: " << k << "\tconvergence: " << converge << std::endl;
	if (converge > param->epsilon) std::clog << "\033[1A\033[K";
	return 0;
}

int main(int argc, char const *argv[])
{
	srand(time(0));
	// 添加一些大数
	int tmp_id, tmp_size;
	double tmp_val;
	for (int i = 0; i < M; i++)
	{
		tmp_size = gctl::random(25, 35);
		for (int j = 0; j < tmp_size; j++)
		{
			tmp_id = gctl::random(0, N);
			tmp_val = gctl::random(-10.0, 10.0);

			kernel[i][tmp_id] = tmp_val;
			sp_kernel.insert(i, tmp_id, tmp_val);
		}
	}

	// 生成一组正演解
	gctl::array<double> fm(N);
	for (int i = 0; i < N; i++)
	{
		fm[i] = gctl::random(1.0, 2.0);
	}

	// 计算共轭梯度B项
	gctl::array<double> B(N);
	sp_kernel.multiply_vector(fm.get(), N, tmp_arr.get(), M);
	sp_kernel.multiply_vector(tmp_arr.get(), M, B.get(), N, gctl::Trans);
	/*
	for (int i = 0; i < M; i++)
	{
		tmp_arr[i] = 0.0;
		for (int j = 0; j < N; j++)
		{
			tmp_arr[i] += kernel[i][j]*fm[j];
		}
	}

	for (int j = 0; j < N; j++)
	{
		B[j] = 0.0;
		for (int i = 0; i < M; i++)
		{
			B[j] += kernel[i][j]*tmp_arr[i];
		}
	}
	*/

	/********************准备工作完成************************/
	lcg_para self_para = lcg_default_parameters();
	self_para.max_iterations = 1000;
	self_para.epsilon = 1e-10;

	// 声明两组解
	gctl::array<double> m(N, 0.0);
	gctl::array<double> m_sp(N, 0.0);

	clock_t start = clock();
	int ret = lcg_solver(CalAx, Prog, m.get(), B.get(), N, &self_para, NULL, LCG_CG);
	clock_t end = clock();
	if (ret < 0) lcg_error_str(ret);
	std::cout << "array2d's time: " << 1000.0*(end - start)/(double)CLOCKS_PER_SEC << " ms" << std::endl;

	start = clock();
	ret = lcg_solver(CalAx_Spmat, Prog, m_sp.get(), B.get(), N, &self_para, NULL, LCG_CG);
	if (ret < 0) lcg_error_str(ret);
	end = clock();
	std::cout << "spmat's time: " << 1000.0*(end - start)/(double)CLOCKS_PER_SEC << " ms" << std::endl;

	for (int i = 0; i < N; i++)
	{
		std::cout << fm[i] << " " << m[i] << " " << m_sp[i] << std::endl;
	}
	return 0;
}