rlsmatrix.h

很好的矩阵链的计算代码

void AddSMatrix(RLSMatrix a,RLSMatrix b,RLSMatrix &c)
{
	//求两个用三元组顺序表表示的二维矩阵a与b的和，并用新的矩阵c返回

	int i,j,s;
	int t,tp,arow,ccol;
	int ctemp[21];

	if((a.mu!=b.mu) || (b.nu!=a.nu))                  //判断是否可以进行加法运算
	{
		printf("add error\n");
		exit(0);
	}
	c.mu=a.mu;c.nu=a.nu;c.tu=0;
	if(a.tu+b.tu==0)
		return;
	
	for(arow=1;arow<=a.mu;arow++)
	{
		for(i=1;i<=a.nu;i++)
			ctemp[i]=0;

		c.rpos[arow]=c.tu+1;
		if(arow<a.mu)
		{
			tp=a.rpos[arow+1];
			t =b.rpos[arow+1];
		}
		else
		{
			tp=a.tu+1;
			t =b.tu+1;
		}
		j=a.rpos[arow];s=b.rpos[arow];
		//在一行中进行同时遍历，当a中的一个元素的列数更小时赋给ctemp数组，更大则b的赋给数组，否则运算后赋给数组
		while(j<tp && s<t)
		{
			if(a.data[j].col<b.data[s].col)
			{
				ccol       =a.data[j].col;
				ctemp[ccol]=a.data[j].number;
				j++;
			}
			else if(a.data[j].col>b.data[s].col)
			{
				ccol       =b.data[s].col;
				ctemp[ccol]=b.data[s].number;
				s++;
			}
			else
			{
				ccol       =a.data[j].col;
				ctemp[ccol]=a.data[j].number+b.data[s].number;
				++j;s++;
			}
			
		}
		while(j<tp)
		{
			ccol       =a.data[j].col;
			ctemp[ccol]=a.data[j].number;
			j++;
		}
		while(s<t)
		{
			ccol       =b.data[s].col;
			ctemp[ccol]=b.data[s].number;
			s++;
		}

		for(ccol=1;ccol<=c.nu;ccol++)
			if(ctemp[ccol])
			{
				if(++c.tu>MAXSIZE)
				{
					printf("add bounds error\n");
					exit(0);
				}
				c.data[c.tu].row   =arow;
				c.data[c.tu].col   =ccol;
				c.data[c.tu].number=ctemp[ccol];
			}
	}
}

void SubSMatrix(RLSMatrix a,RLSMatrix b,RLSMatrix &c)
{
	//求两个用三元组顺序表表示的二维矩阵a与b的差，并用新的矩阵c返回

	int i,j,s;
	int t,tp,arow,ccol;
	int ctemp[21];

	if((a.mu!=b.mu) || (b.nu!=a.nu))
	{
		printf("sub error\n");
		exit(0);
	}
	c.mu=a.mu;c.nu=a.nu;c.tu=0;
	if(a.tu+b.tu==0)
		return;
	
	for(arow=1;arow<=a.mu;arow++)
	{
		for(i=1;i<=a.nu;i++)
			ctemp[i]=0;

		c.rpos[arow]=c.tu+1;
		if(arow<a.mu)
		{
			tp=a.rpos[arow+1];
			t =b.rpos[arow+1];
		}
		else
		{
			tp=a.tu+1;
			t =b.tu+1;
		}
		j=a.rpos[arow],s=b.rpos[arow];
		while(j<tp && s<t)
		{
			if(a.data[j].col<b.data[s].col)
			{
				ccol       =a.data[j].col;
				ctemp[ccol]=a.data[j].number;
				j++;
			}
			else if(a.data[j].col>b.data[s].col)
			{
				ccol       =b.data[s].col;
				ctemp[ccol]=0-b.data[s].number;
				s++;
			}
			else
			{
				ccol       =a.data[j].col;
				ctemp[ccol]=a.data[j].number-b.data[s].number;
				j++;s++;
			}
			
		}
		while(j<tp)
		{
			ccol       =a.data[j].col;
			ctemp[ccol]=a.data[j].number;
			j++;
		}
		while(s<t)
		{
			ccol       =b.data[s].col;
			ctemp[ccol]=0-b.data[s].number;
			s++;
		}

		for(ccol=1;ccol<=c.nu;ccol++)
			if(ctemp[ccol])
			{
				if(++c.tu>MAXSIZE)
				{
					printf("sub bounds error\n");
					exit(0);
				}
				c.data[c.tu].row   =arow;
				c.data[c.tu].col   =ccol;
				c.data[c.tu].number=ctemp[ccol];
			}
	}
}

void MulSMatrix(RLSMatrix a,RLSMatrix b,RLSMatrix &c)
{
	//求两个用三元组顺序表表示的二维矩阵a与b的乘积，并用新的矩阵c返回

	int i,j,t,tp;
	int ccol,arow,brow;
	int ctemp[21];

	if(a.nu!=b.mu)
	{
		printf("mul error\n");
		exit(0);
	}
	c.mu=a.mu;c.nu=b.nu;c.tu=0;
	if(a.tu*b.tu==0)
		return;
	for(arow=1;arow<=a.mu;arow++)
	{
		for(i=1;i<=b.nu;i++)
			ctemp[i]=0;
		c.rpos[arow]=c.tu+1;
		if(arow<a.mu)
			tp=a.rpos[arow+1];
		else 
			tp=a.tu+1;
		for(i=a.rpos[arow];i<tp;i++)
		{
			brow=a.data[i].col;
			if(brow<b.mu)
				t=b.rpos[brow+1];
			else
				t=b.tu+1;
			for(j=b.rpos[brow];j<t;j++)
			{
				ccol=b.data[j].col;
				ctemp[ccol]+=a.data[i].number*b.data[j].number;

			}

		}
		for(ccol=1;ccol<=c.nu;ccol++)
			if(ctemp[ccol])
			{
				if(++c.tu>MAXSIZE)
				{
					printf("bounds error\n");
					exit(0);
				}
				c.data[c.tu].row   =arow;
				c.data[c.tu].col   =ccol;
				c.data[c.tu].number=ctemp[ccol];
			}
	}
}

void TransposeSMatrix(RLSMatrix a,RLSMatrix &t)
{
	//求矩阵a的转置矩阵，并用新的矩阵t返回

	int i,j,acol;
	int num[21],cpot[21];

	t.mu=a.nu; t.nu=a.mu; t.tu=a.tu;
	if((a.mu==0) || (a.nu==0))
		return;
	if(t.tu==0)
		return;
	for(acol=1;acol<=a.nu;acol++)
		num[acol]=0;
	for(i=1;i<=a.tu;i++)
		++num[a.data[i].col];
	cpot[1]=1;
	for(acol=2;acol<=a.nu;acol++)
		cpot[acol]=cpot[acol-1]+num[acol-1];
	for(i=1;i<=a.tu;i++)
	{
		acol=a.data[i].col;
		j   =cpot[acol];

		t.data[j].row   =a.data[i].col;
		t.data[j].col   =a.data[i].row;
		t.data[j].number=a.data[i].number;
		cpot[acol]++;
	}

}

RLSMatrix RemainderSMatrix(RLSMatrix a,int delrow,int delcol)
{
	//求除去第delrow行和第delcol列的所有元素后的剩余矩阵，用矩阵b返回

	int i,tp,arow;
	RLSMatrix b;

	if(a.mu!=a.nu)
	{
		printf("输入矩阵应为n阶方阵\n");
		exit(0);
	}
	if((delrow>a.mu) || (delcol>a.nu))          //当要删除的行数或列数不合法时的处理
	{
		printf("remainder error\n");
		exit(0);
	}
	if(a.mu==0)
	{
		printf("矩阵为空\n");
		exit(0);
	}
	b.mu=a.mu-1;b.nu=a.nu-1;b.tu=0;
	//先处理一行，根据该行中元素的列数和要删除列数比较来进行不同的处理
	for(arow=1;arow<=a.mu;arow++)
	{
		if(arow<a.mu)                   
			tp=a.rpos[arow+1];
		else
			tp=a.tu+1;
		if(arow==delrow)continue;       //列数与目标列数相等则不做处理
		if(arow<delrow)
		{
			b.rpos[arow]=b.tu+1;          //列数比目标列数更小，则直接复制，否则行数或列数需减一
		}
		else
		{
			b.rpos[arow-1]=b.tu+1;
		}
		for(i=a.rpos[arow];i<tp;i++)
		{
			if(a.data[i].col==delcol)continue;
			if(arow<delrow)
			{
				b.data[b.tu+1].row   =arow;
			}
			else
			{
				b.data[b.tu+1].row   =arow-1;
			}
			if(a.data[i].col<delcol)
			{
					
				b.data[b.tu+1].col=a.data[i].col;
					
			}
			else
			{
				b.data[b.tu+1].col=a.data[i].col-1;

			}
			b.data[b.tu+1].number=a.data[i].number;
			b.tu++;
		}
	}
	return b;
}

int ValueSMatrix(RLSMatrix a)
{
	//求矩阵a的行列式的值,并用value返回
	int value,i,arow,tp,m;
	
	value=0;
	arow=1;
	if(a.mu!=a.nu)
	{
		printf("请输入方阵\n");
		exit(0);
	}
	if(a.tu==0)
		return 0;
	if(a.mu==1)
		return a.data[1].number;
	if(arow<a.mu)
		tp=a.rpos[arow+1];
	else 
		tp=a.tu+1;
	//直接按第一行展开来算行列式
	for(i=a.rpos[arow];i<tp;i++)
	{
		
		m=(a.data[i].row+a.data[i].col)%2;

		if(m==1)         //该元素的行数和列数之和为奇数
		{
			value+=a.data[i].number*(-1)*ValueSMatrix(RemainderSMatrix(a,arow,a.data[i].col));

		}
		else             //该元素的行数和列数之和为偶数
			value+=a.data[i].number*ValueSMatrix(RemainderSMatrix(a,arow,a.data[i].col));
	}
	return value;
}

void GoWithSMatrix(RLSMatrix a,RLSMatrix &b)
{
	//求矩阵a的伴随矩阵，并用矩阵b返回

	int ctemp[21][21];
	int arow,acol,m;

	for(arow=1;arow<=a.mu;arow++)
	{
		for(acol=1;acol<=a.nu;acol++)
		{
			m=(arow+acol)%2;
			if(m==0)           //该元素的行数和列数之和为偶数
			{
				ctemp[acol][arow]=ValueSMatrix(RemainderSMatrix(a,arow,acol));
			}
			else                 //该元素的行数和列数之和为奇数
				ctemp[acol][arow]=(-1)*ValueSMatrix(RemainderSMatrix(a,arow,acol));
		}
	}
	b.mu=a.mu;b.nu=a.nu;b.tu=0;
	//将ctemp中保存的伴随矩阵中的元素压入稀疏矩阵进行存储
	for(arow=1;arow<=a.mu;arow++)
	{
		b.rpos[arow]=b.tu+1;
		for(acol=1;acol<=a.nu;acol++)
			if(ctemp[arow][acol])
			{
				if(++b.tu>MAXSIZE)
				{
					printf("create bounds error\n");
					exit(0);
				}
				b.data[b.tu].row   =arow;
				b.data[b.tu].col   =acol;
				b.data[b.tu].number=ctemp[arow][acol];
			}
	}
}

void AthwartSMatrix(RLSMatrix a,RLSMatrix &b)
{
	//求矩阵a的逆矩阵，并用新的矩阵b返回

	int value,i;

	value=ValueSMatrix(a);
	if(value==0)
	{
		printf("该矩阵无逆矩阵\n");
		exit(0);
	}
	GoWithSMatrix(a,b);         //用b返回a的伴随矩阵         
	for(i=1;i<=b.tu;i++)
	{
		b.data[i].number=b.data[i].number/value;          //将伴随矩阵的每个元素再除以行列式的值即可
	}
}