#include "denseBlas.h"

/* 创造一个矩阵，分配新的矩阵空间 */
demat* createDeMat(int m, int n) // m行，n列，nnz个非零元素
{
	int numEle;
	numEle = m * n;
	double* v = (double*)malloc(numEle * sizeof(double));		// 全部元素
	return assembleDeMat(m, n, v);
}

/* 利用已有空间组装一个矩阵 */
demat* assembleDeMat(int m, int n, double* v) // m行，n列，v为元素，
{
	demat* M = (demat*)malloc(sizeof(demat));
	M->m = m;
	M->n = n;
	M->v = v;
	return M;
}
/* 释放一个稠密矩阵空间 */
void freeDeMat(demat* M, int level)
{
	if (M->v) {
		free(M->v);
	}
	if (level == 2) {
		free(M);
	}
}

/* 2范数 */
double cmscp_norm2(double* v, int n) {
	int i;
	double normsquare = 0;
	for (i = 0; i < n; i++) {
		normsquare += v[i] * v[i];
	}
	return sqrt(normsquare);
}

/* 稠密矩阵与矩阵相乘 */
void denseMM(demat* matOut, demat* matLeftIn, demat* matRightIn) {
	/* 变量类型 */
	int mMatOut, nMatOut, mMatLeftIn, nMatLeftIn, mMatRightIn, nMatRightIn, i, j, k, p1, p2, p3;
	double* vMatOut, * vMatLeftIn, * vMatRightIn;

	/* 维数信息 */
	mMatOut = matLeftIn->m;
	nMatOut = matRightIn->n;
	mMatLeftIn = matLeftIn->m;
	nMatLeftIn = matLeftIn->n;
	mMatRightIn = matRightIn->m;
	nMatRightIn = matRightIn->n;

	/* 矩阵元素 */
	vMatOut = matOut->v;
	vMatLeftIn = matLeftIn->v;
	vMatRightIn = matRightIn->v;

	/* 运算 */
	for (i = 0; i < mMatOut; i++) {
		for (j = 0; j < nMatOut; j++) {
			p1 = i * mMatOut + j; // matOut中第（i,j）元素在向量存储格式中的位置
			vMatOut[p1] = 0;
			for (k = 0; k < nMatLeftIn; k++) {
				p2 = i * mMatLeftIn + k; // matLeftIn中第（i,k）元素在向量存储格式中的位置
				p3 = k * mMatRightIn + j; // matRightIn中第（k,j）元素在向量存储格式中的位置
				vMatOut[p1] += vMatLeftIn[p2] * vMatRightIn[p3];
			}
		}
	}
}


/* 稠密矩阵转置与矩阵相乘 */
void denseMtM(demat* matOut, demat* matLeftIn, demat* matRightIn) {
	/* 变量类型 */
	int mMatOut, nMatOut, mMatLeftIn, nMatLeftIn, mMatRightIn, nMatRightIn, i, j, k, p1, p2, p3;
	double* vMatOut, * vMatLeftIn, * vMatRightIn;

	/* 维数信息 */
	mMatOut = matLeftIn->m;
	nMatOut = matRightIn->n;
	mMatLeftIn = matLeftIn->m;
	nMatLeftIn = matLeftIn->n;
	mMatRightIn = matRightIn->m;
	nMatRightIn = matRightIn->n;

	/* 矩阵元素 */
	vMatOut = matOut->v;
	vMatLeftIn = matLeftIn->v;
	vMatRightIn = matRightIn->v;

	/* 运算 */
	for (i = 0; i < mMatOut; i++) {
		for (j = 0; j < nMatOut; j++) {
			p1 = i * mMatOut + j; // matOut中第（i,j）元素在向量存储格式中的位置
			vMatOut[p1] = 0;
			for (k = 0; k < nMatLeftIn; k++) {
				p2 = k * mMatLeftIn + i; // matLeftIn转置矩阵中第（i,k）元素在向量存储格式中的位置
				p3 = k * mMatRightIn + j; // matRightIn中第（k,j）元素在向量存储格式中的位置
				vMatOut[p1] += vMatLeftIn[p2] * vMatRightIn[p3];
			}
		}
	}
}


/* 稠密矩阵与矩阵转置相乘 */
void denseMMt(demat* matOut, demat* matLeftIn, demat* matRightIn) {
	/* 变量类型 */
	int mMatOut, nMatOut, mMatLeftIn, nMatLeftIn, mMatRightIn, nMatRightIn, i, j, k, p1, p2, p3;
	double* vMatOut, * vMatLeftIn, * vMatRightIn;

	/* 维数信息 */
	mMatOut = matLeftIn->m;
	nMatOut = matRightIn->n;
	mMatLeftIn = matLeftIn->m;
	nMatLeftIn = matLeftIn->n;
	mMatRightIn = matRightIn->m;
	nMatRightIn = matRightIn->n;

	/* 矩阵元素 */
	vMatOut = matOut->v;
	vMatLeftIn = matLeftIn->v;
	vMatRightIn = matRightIn->v;

	/* 运算 */
	for (i = 0; i < mMatOut; i++) {
		for (j = 0; j < nMatOut; j++) {
			p1 = i * mMatOut + j; // matOut中第（i,j）元素在向量存储格式中的位置
			vMatOut[p1] = 0;
			for (k = 0; k < nMatLeftIn; k++) {
				p2 = i * mMatLeftIn + k; // matLeftIn矩阵中第（i,k）元素在向量存储格式中的位置
				p3 = j * mMatRightIn + k; // matRightIn转置中第（k,j）元素在向量存储格式中的位置
				vMatOut[p1] += vMatLeftIn[p2] * vMatRightIn[p3];
			}
		}
	}
}

/* 稠密矩阵与向量相乘 */
void denseMV(double* vecOut, demat* matLeftIn, double* vecRightIn, int sign, int newVecOut) {
	/* 变量类型 */
	int mMatLeftIn, nMatLeftIn, i, j, p;
	double* vMatLeftIn;

	/* 维数信息 */
	mMatLeftIn = matLeftIn->m;
	nMatLeftIn = matLeftIn->n;

	/* 输入矩阵元素 */
	vMatLeftIn = matLeftIn->v;

	if (newVecOut == 1) // 
	{
		for (i = 0; i < mMatLeftIn; i++)
		{
			vecOut[i] = 0;
		}
	}
	if (sign == 1) // 加法
	{
		for (i = 0; i < mMatLeftIn; i++)
		{
			for (j = 0; j < nMatLeftIn; j++)
			{
				p = i * mMatLeftIn + j;
				vecOut[i] += vMatLeftIn[p] * vecRightIn[j];
			}
		}
	}
	else // 减法
	{
		for (i = 0; i < mMatLeftIn; i++)
		{
			for (j = 0; j < nMatLeftIn; j++)
			{
				p = i * mMatLeftIn + j;
				vecOut[i] -= vMatLeftIn[p] * vecRightIn[j];
			}
		}
	}
}


/* 稠密矩阵与向量相乘 */
void denseMtV(double* vecOut, demat* matLeftIn, double* vecRightIn, int sign, int newVecOut) {
	/* 变量类型 */
	int mMatLeftIn, nMatLeftIn, i, j, p;
	double* vMatLeftIn;

	/* 维数信息 */
	mMatLeftIn = matLeftIn->m;
	nMatLeftIn = matLeftIn->n;

	/* 输入矩阵元素 */
	vMatLeftIn = matLeftIn->v;

	if (newVecOut == 1) // 
	{
		for (j = 0; j < mMatLeftIn; j++)
		{
			vecOut[j] = 0;
		}
	}
	if (sign == 1) // 加法
	{
		for (i = 0; i < mMatLeftIn; i++)
		{
			for (j = 0; j < nMatLeftIn; j++)
			{
				p = i * mMatLeftIn + j; // matLeftIn的转置矩阵的(j,i)元素
				vecOut[j] += vMatLeftIn[p] * vecRightIn[i];
			}
		}
	}
	else // 减法
	{
		for (i = 0; i < mMatLeftIn; i++)
		{
			for (j = 0; j < nMatLeftIn; j++)
			{
				p = i * mMatLeftIn + j; // matLeftIn的转置矩阵的(j,i)元素
				vecOut[j] -= vMatLeftIn[p] * vecRightIn[i];
			}
		}
	}
}


/* 常值与向量相乘 */
void denseScaV(double* vecOut, double scaIn, double* vecIn, int length) {
	/* 变量类型 */
	int i;
	for (i = 0; i < length; i++)
	{
		vecOut[i] = scaIn * vecIn[i];
	}
}

/* 常值与向量相加 */
void denseScaPlusVec(double* vecOut, double scaIn, double* vecIn, int length) {
	/* 变量类型 */
	int i;
	for (i = 0; i < length; i++)
	{
		vecOut[i] = scaIn + vecIn[i];
	}
}