稀疏矩阵的快速转制c语言代码_有限元方法：从代码入手(4)

最新推荐文章于 2024-08-06 23:03:21 发布

原创最新推荐文章于 2024-08-06 23:03:21 发布 · 496 阅读

4 ·

CC 4.0 BY-SA版权

文章标签：

#稀疏矩阵的快速转制c语言代码

本文介绍了如何使用C语言封装稀疏矩阵，特别是行主元存储（CSR）格式。针对有限元方法中得到的稀疏矩阵，只存储非零元素以节省内存。文中详细解释了Value、ColI和RowI数组的作用和存储方式，并说明了矩阵加法、数乘和矩阵向量乘法等操作的实现。同时，提到了FEMatrix类作为连接有限元单元和稀疏矩阵的桥梁。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

稀疏矩阵的封装

由有限元方法获得的单元刚度矩阵通常与单个单元相关联。这样单元矩阵与有限个单元节点相互关联，由这些单元矩阵组装成的全局矩阵就含有大量的非零元，即稀疏矩阵。为了节省存储内存，通常仅仅存储非零元的信息。稀疏矩阵的具有多种存储格式，包括行主元存储，列主元存储，坐标存储格式。

在这里封装时，仅仅对CSR格式的稀疏阵进行封装。以上图所示的稀疏矩阵为例，在存储的时候，需要一个名为Value的数组存储非零元。对于CSR格式，按照行的优先顺序进行存储，相应的Value数组内存储的内容为：

此外还需另外一个数组ColI，用于存储各个非零元所对应的列索引。对于C++，数组的索引是从零开始的，相应列索引中的存储内容为：

随后，还需要一个RowI数组存储行索引。与列索引不同，这里行索引存储RowI[0]存储的是第零行的第一个元素在Value中的索引，索引RowI[0]中存储着0； RowI[1]存储着第1行的第一个元素在Value中的索引，这样RowI[1]中存储着3。对于其他元素，相应的存储内容依此类推。并且，RowI在最终还会与一个额外的存储内容，为所有非零元素的数目。这样，最终RowI的存储内容为。

在代码中，也封装了矩阵的加法，数乘，矩阵与向量的乘法运算。

在这个代码中，还具有另外一个FEMatrix，作为联系有限元单元和稀疏矩阵存储格式的桥梁，并初始化稀疏矩阵的ColI和RowI内容，这些内容将在下一小节详细论述。

//SparseMatrix.h

#pragma once
#include"mkl.h"
#include"VECTOR.h"
#include"FEMatrix.h"

class SparseMatrix
{
public:

	//稀疏矩阵所表示矩阵的行
	MKL_INT Row = 0;

	//稀疏矩阵所表示矩阵的列
	MKL_INT Col = 0;

	//稀疏矩阵的非零元素的数目
	MKL_INT Nonzeros = 0;

	//稀疏矩阵行索引数组
	MKL_INT* RowIndex = NULL;

	//稀疏矩阵的列索引数组
	MKL_INT* ColIndex = NULL;

	//稀疏矩阵的的压缩存储一维向量
	double* Value = NULL;


////////////////////////////////////////////////////////////////////

//初始化稀疏矩阵并分配内存
	void IniSP(MKL_INT row, MKL_INT col, MKL_INT nonzeros);

	//由FEMatrix初始化稀疏矩阵
	void IniSP(FEMatrix& fm);
/////////////////////////////////////////////////////////////////////
	//获得稀疏阵的信息
	//获得行数
	MKL_INT GetRow();

	//获得列数
	MKL_INT GetCol();

	//获得非零元数
	MKL_INT GetNonzeros();

	//获得行索引指针
	MKL_INT* GetPRow();

	//获得列索引指针
	MKL_INT* GetPCol();

	//获得非零元指针
	double* GetPNonzeros();
///////////////////////////////////////////////////////////////////////
	//设定sparse中的一个元素,该元素的行为row,列为column，更新为value
	void SetSingleValue(MKL_INT row, MKL_INT column, double value);

	//更新sparse中的一个元素,该元素的行为row,列为column，更新为value
	void AddSingleValue(MKL_INT row, MKL_INT column, double value);

	//获得稀疏矩阵存储区域第ii个索引的值
	double GetSingleValue(MKL_INT ii);

	//获得稀疏矩阵第ii行的起始索引值
	MKL_INT GetRowIndex(MKL_INT ii);

	//获得稀疏阵第ii个存储位置的列索引
	MKL_INT GetColIndex(MKL_INT ii);

	//设定稀疏矩阵存储区域第ii个索引的值为Value
	void SetSingleValue(MKL_INT ii, double value);
///////////////////////////////////////////////////////////////////////
	//清零稀疏阵中的所有元素
	void Zeros();

	//将稀疏矩阵的格式复制给另外一个矩阵
	void CopyFormat(char c, SparseMatrix & sm);

	//将稀疏矩阵复制给另外一个矩阵
	void CopyMatrix(char c, SparseMatrix& sm);

	//将稀疏矩阵加上一个稀疏阵
	void AddMatrix(SparseMatrix& sm);

	//将稀疏矩阵加上一个稀疏阵
	void MinusMatrix(SparseMatrix& sm);

	//将当前稀疏矩阵乘上一个常数
	void MutiConst(double b);

	//稀疏矩阵乘以矩阵并返回一个矩阵
	void MutiVec(VECTOR& v, VECTOR& res);


	//求解稀疏矩阵
	void SolverPARDISO(VECTOR& v, VECTOR& res);
	void SolverDSS(VECTOR& v, VECTOR& res);

	//收回分配的内存
	void Clear();
	SparseMatrix();
	SparseMatrix(MKL_INT row, MKL_INT col, MKL_INT nonzeros);
	~SparseMatrix();
};

//SparseMatrix.cpp

#include "SparseMatrix.h"
#include <iostream>
#include"mkl.h"
#include"mkl_dss.h"
#include <fstream>
#include <cmath>
/////////////////////////////////////////////////////////////////////
//初始化稀疏矩阵并分配内存
void SparseMatrix::IniSP(MKL_INT row, MKL_INT col, MKL_INT nonzeros)
{
	Clear();
	Row = row;
	Col = col;
	Nonzeros = nonzeros;
	if (RowIndex == NULL)
	{
		RowIndex = new MKL_INT[Row + 1];
		memset(RowIndex, 0, (Row + 1) * sizeof(MKL_INT));
	}
	if (ColIndex == NULL)
	{
		ColIndex = new MKL_INT[Nonzeros];
		memset(ColIndex, 0, (Nonzeros) * sizeof(MKL_INT));
	}
	if (Value == NULL)
	{
		Value = new double[Nonzeros];
		memset(Value, 0, Nonzeros * sizeof(double));
	}
}

//由FEMatrix初始化稀疏矩阵
void SparseMatrix::IniSP(FEMatrix& fm)
{
	IniSP(fm.GetRow(), fm.GetCol(), fm.GetNonzeros());

	RowIndex[0] = 0;
	for (MKL_INT rowii = 0; rowii < Row; rowii++)
	{
		MKL_INT col_start = RowIndex[rowii];
		MKL_INT row_length = fm.GetRowLength(rowii);
		for (MKL_INT coljj = 0; coljj < row_length; coljj++)
		{
			ColIndex[col_start + coljj] = fm.GetMatrixElement(rowii, coljj);
		}
		RowIndex[rowii + 1] = RowIndex[rowii] + row_length;
	}
}
/////////////////////////////////////////////////////////////////////
//获得稀疏阵的信息
	//获得行数
MKL_INT SparseMatrix::GetRow()
{
	return Row;
}

//获得列数
MKL_INT SparseMatrix::GetCol()
{
	return Col;
}

//获得非零元数
MKL_INT SparseMatrix::GetNonzeros()
{
	return Nonzeros;
}

//获得行索引指针
MKL_INT* SparseMatrix::GetPRow()
{
	return RowIndex;
}

//获得列索引指针
MKL_INT* SparseMatrix::GetPCol()
{
	return ColIndex;
}

//获得非零元指针
double* SparseMatrix::GetPNonzeros()
{
	return Value;
}

///////////////////////////////////////////////////////////////////////
//设定sparse中的一个元素,该元素的行为row,列为column，更新为value
void SparseMatrix::SetSingleValue(MKL_INT row, MKL_INT column, double value)
{
	MKL_INT s = 0;
	MKL_INT row_loc = RowIndex[row];
	MKL_INT row_length = RowIndex[row + 1] - row_loc;
	for (MKL_INT ii = 0; ii < row_length; ii++)
	{
		s = row_loc + ii;
		if (ColIndex[s] == column)
		{
			Value[s] = value;
			return;
		}
		else continue;
	}
	std::cout.clear();
	std::cout << "SparseMatrix::SetSingleValue ERROR" << std::endl;
	getchar();
	exit(0);
}

//更新sparse中的一个元素,该元素的行为row,列为column，更新为value
void SparseMatrix::AddSingleValue(MKL_INT row, MKL_INT column, double value)
{
	MKL_INT s = 0;
	MKL_INT row_loc = RowIndex[row];
	MKL_INT row_length = RowIndex[row + 1] - row_loc;
	for (MKL_INT ii = 0; ii < row_length; ii++)
	{
		s = row_loc + ii;
		if (ColIndex[s] == column)
		{
			Value[s] += value;
			return;
		}
		else continue;
	}
	std::cout.clear();
	std::cout << "SparseMatrix::SetSingleValue ERROR" << std::endl;
	getchar();
	exit(0);
}


//获得稀疏矩阵存储区域第ii个索引的值
double SparseMatrix::GetSingleValue(MKL_INT ii)
{
	if (ii < Nonzeros)
		return Value[ii];
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::GetSingleVale ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//赋予稀疏矩阵存储区域第ii个索引的值为value
void SparseMatrix::SetSingleValue(MKL_INT ii, double value)
{
	if (ii < Nonzeros)
		Value[ii] = value;
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::GetSingleVale ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//获得稀疏矩阵第ii行的起始索引值
MKL_INT SparseMatrix::GetRowIndex(MKL_INT ii)
{
	if (ii <= Row)
		return RowIndex[ii];
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::GetRowIndex ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//获得稀疏阵第ii个存储位置的列索引
MKL_INT SparseMatrix::GetColIndex(MKL_INT ii)
{
	return ColIndex[ii];
}



//清零稀疏阵中的所有元素
void SparseMatrix::Zeros()
{
	memset(Value, 0, Nonzeros * sizeof(double));
}

//将稀疏矩阵的格式复制给另外一个矩阵
void SparseMatrix::CopyFormat(char c, SparseMatrix& sm)
{
	if (c == '<')
	{
		IniSP(sm.GetRow(), sm.GetCol(), sm.GetNonzeros());
		memcpy(ColIndex, sm.GetPCol(), Nonzeros * sizeof(MKL_INT));
		memcpy(RowIndex, sm.GetPRow(), (Row + 1) * sizeof(MKL_INT));
		Zeros();
	}
	else if (c == '>')
	{
		sm.IniSP(Row, Col, Nonzeros);
		memcpy(sm.GetPCol(), ColIndex, Nonzeros * sizeof(MKL_INT));
		memcpy(sm.GetPRow(), RowIndex, (Row + 1) * sizeof(MKL_INT));
		sm.Zeros();
	}
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::CopyFormat ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//将稀疏矩阵复制给另外一个矩阵
void SparseMatrix::CopyMatrix(char c, SparseMatrix& sm)
{
	if (c == '<')
	{
		CopyFormat('<', sm);
		memcpy(Value, sm.GetPNonzeros(), Nonzeros * sizeof(double));
	}
	else if (c == '>')
	{
		CopyFormat('>', sm);
		memcpy(sm.GetPNonzeros(), Value, Nonzeros * sizeof(double));
	}
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::CopyMatrix ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//将稀疏矩阵加上一个稀疏阵
void SparseMatrix::AddMatrix(SparseMatrix& sm)
{
	if ((sm.GetRow() == Row) && (sm.GetCol() == Col) && sm.GetNonzeros() == Nonzeros)
	{
		cblas_daxpy(Nonzeros, 1.0, sm.GetPNonzeros(), 1, Value, 1);
	}
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::AddMatrix ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//将稀疏矩阵减去一个稀疏阵
void SparseMatrix::MinusMatrix(SparseMatrix& sm)
{
	if ((sm.GetRow() == Row) && (sm.GetCol() == Col) && sm.GetNonzeros() == Nonzeros)
	{
		cblas_daxpy(Nonzeros, -1.0, sm.GetPNonzeros(), 1, Value, 1);
	}
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::AddMatrix ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//将当前稀疏矩阵乘上一个常数
void SparseMatrix::MutiConst(double b)
{
	cblas_dscal(Nonzeros, b, Value, 1);
}

//稀疏矩阵乘以矩阵并返回一个矩阵
void SparseMatrix::MutiVec(VECTOR& v, VECTOR& res)
{
	if (Col == v.GetVecDim())
	{
		VECTOR nv;
		nv.InitVec(Row);
		// Descriptor of main sparse matrix properties
		struct matrix_descr descrA;
		// // Structure with sparse matrix stored in CSR format
		sparse_matrix_t       csrA;

		// Create handle with matrix stored in CSR format
		mkl_sparse_d_create_csr(&csrA, SPARSE_INDEX_BASE_ZERO,
			Row,  // number of rows
			Col,  // number of cols
			RowIndex,
			RowIndex + 1,
			ColIndex,
			Value);

		// Create matrix descriptor
		descrA.type = SPARSE_MATRIX_TYPE_GENERAL;
		// Analyze sparse matrix; choose proper kernels and workload balancing strategy
		mkl_sparse_optimize(csrA);
		// Compute y = alpha * A * x + beta * y
		mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE,
			1.0,
			csrA,
			descrA,
			v.GetVecP(),
			0,
			nv.GetVecP());
		// Release matrix handle and deallocate matrix
		mkl_sparse_destroy(csrA);
		res.Clear();
		res.InitVec(nv.GetVecDim());
		res.CopyVec('<', nv);
		nv.Clear();
	}
	else
	{
		std::cout.clear();
		std::cout << "SparseMatrix::MutiVec ERROR" << std::endl;
		getchar();
		exit(0);
	}
}

//求解稀疏矩阵
void SparseMatrix::SolverPARDISO(VECTOR& v, VECTOR& res)
{

	MKL_INT mtype = 11;   /* Real unsymmetric matrix */
	/* RHS and solution vectors. */
	MKL_INT nrhs = 1;     /* Number of right hand sides. */
	/* Internal solver memory pointer pt, */
	/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
	/* or void *pt[64] should be OK on both architectures */
	void* pt[64];
	/* Pardiso control parameters. */
	MKL_INT iparm[64];
	MKL_INT maxfct, mnum, phase, error, msglvl;
	/* Auxiliary variables. */
	MKL_INT idum;         /* Integer dummy. */
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
	memset(iparm, 0, 64 * sizeof(MKL_INT));
	iparm[0] = 1;         /* No solver default */
	iparm[1] = 2;         /* Fill-in reordering from METIS */
	iparm[3] = 0;         /* No iterative-direct algorithm */
	iparm[4] = 0;         /* No user fill-in reducing permutation */
	iparm[5] = 0;         /* Write solution into x */
	iparm[6] = 0;         /* Not in use */
	iparm[7] = 5;         /* Max numbers of iterative refinement steps */
	iparm[8] = 0;         /* Not in use */
	iparm[9] = 13;        /* Perturb the pivot elements with 1E-13 */
	iparm[10] = 1;        /* Use nonsymmetric permutation and scaling MPS */
	iparm[11] = 0;        /* Conjugate transposed/transpose solve */
	iparm[12] = 1;        /* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
	iparm[13] = 0;        /* Output: Number of perturbed pivots */
	iparm[14] = 0;        /* Not in use */
	iparm[15] = 0;        /* Not in use */
	iparm[16] = 0;        /* Not in use */
	iparm[17] = -1;       /* Output: Number of nonzeros in the factor LU */
	iparm[18] = -1;       /* Output: Mflops for LU factorization */
	iparm[19] = 0;        /* Output: Numbers of CG Iterations */
	iparm[34] = 1;		  //Zero based
	maxfct = 1;           /* Maximum number of numerical factorizations. */
	mnum = 1;         /* Which factorization to use. */
	msglvl = 0;           /* Print statistical information  */
	error = 0;            /* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
	memset(pt, 0, 64 * sizeof(void*));
	/* -------------------------------------------------------------------- */
	/* .. Reordering and Symbolic Factorization. This step also allocates */
	/* all memory that is necessary for the factorization. */
	/* -------------------------------------------------------------------- */
	phase = 11;
	double ddum;
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&(this->Row), this->Value, this->RowIndex, this->ColIndex, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
	//if (error != 0)
	//{
	//	printf("nERROR during symbolic factorization: %d", error);
	//	exit(1);
	//}
	//printf("nReordering completed ... ");
	//printf("nNumber of nonzeros in factors = %d", iparm[17]);
	//printf("nNumber of factorization MFLOPS = %d", iparm[18]);
	/* -------------------------------------------------------------------- */
	/* .. Numerical factorization. */
	/* -------------------------------------------------------------------- */
	phase = 22;
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&(this->Row), this->Value, this->RowIndex, this->ColIndex, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0)
	{
		printf("nERROR during numerical factorization: %d", error);
		exit(2);
	}
	/* -------------------------------------------------------------------- */
	/* .. Back substitution and iterative refinement. */
	/* -------------------------------------------------------------------- */
	phase = 33;
	// solving steps: Ax=b
	iparm[11] = 0;        /* Conjugate transposed/transpose solve */
	//printf("nnSolving system with iparm[11] = %d ...n", (int)iparm[11]);
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&(this->Row), this->Value, this->RowIndex, this->ColIndex, &idum, &nrhs, iparm, &msglvl, v.GetVecP(), res.GetVecP(), &error);

	/* -------------------------------------------------------------------- */
	/* .. Termination and release of memory. */
	/* -------------------------------------------------------------------- */
	phase = -1;           /* Release internal memory. */
	PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
		&(this->Row), &ddum, this->RowIndex, this->ColIndex, &idum, &nrhs,
		iparm, &msglvl, &ddum, &ddum, &error);
	return;
}

void SparseMatrix::SolverDSS(VECTOR& v, VECTOR& res)
{
	/* Allocate storage for the solver handle and the right-hand side. */
	MKL_INT nRhs = 1;
	_MKL_DSS_HANDLE_t handle;
	_INTEGER_t error;
	_DOUBLE_PRECISION_t statOut[5], eps = 1e-6;
	MKL_INT opt = MKL_DSS_DEFAULTS;
	opt = opt + MKL_DSS_ZERO_BASED_INDEXING;
	MKL_INT  opt1;
	MKL_INT sym = MKL_DSS_NON_SYMMETRIC;
	MKL_INT type = MKL_DSS_INDEFINITE;
	/* --------------------- */
	/* Initialize the solver */
	/* --------------------- */
	error = dss_create(handle, opt);

	/* ------------------------------------------- */
	/* Define the non-zero structure of the matrix */
	/* ------------------------------------------- */
	error = dss_define_structure(handle, sym, RowIndex, Row, Col, ColIndex, Nonzeros);

	/* ------------------ */
	/* Reorder the matrix */
	/* ------------------ */
	//error = dss_reorder(handle, opt, 0);
	opt1 = MKL_DSS_AUTO_ORDER;
	error = dss_reorder(handle, opt1, 0);

	/* ------------------ */
	/* Factor the matrix  */
	/* ------------------ */
	error = dss_factor_real(handle, type, Value);

	/* ------------------------ */
	/* Get the solution vector for Ax=b and ATx=b and check correctness */
	/* ------------------------ */
	// Apply trans or non-trans option, solve system
	opt = MKL_DSS_DEFAULTS;
	
	error = dss_solve_real(handle, opt, v.GetVecP(), nRhs, res.GetVecP());


	/* -------------------------- */
	/* Deallocate solver storage  */
	/* -------------------------- */
	opt = MKL_DSS_MSG_LVL_WARNING + MKL_DSS_TERM_LVL_ERROR;
	error = dss_delete(handle, opt);


}


//////////////////////////////////////////////////////////////////////
//收回分配的内存
//构造与析构函数
void SparseMatrix::Clear()
{
	Row = 0;
	Col = 0;
	Nonzeros = 0;
	if (ColIndex != NULL)
	{
		delete[] ColIndex;
		ColIndex = NULL;
	}
	if (RowIndex != NULL)
	{
		delete[] RowIndex;
		RowIndex = NULL;
	}
	if (Value != NULL)
	{
		delete[] Value;
		Value = NULL;
	}
}
SparseMatrix::SparseMatrix(MKL_INT row, MKL_INT col, MKL_INT nonzeros)
{
	Clear();
	IniSP(row, col, nonzeros);
}
SparseMatrix::SparseMatrix()
{

}
SparseMatrix::~SparseMatrix()
{
	Clear();
}