matrixinv.cpp

电力系统发输电系统基于非序贯的蒙特卡罗抽样的风险评估程序

//矩阵求逆函数
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
void matrixinv(double *pI,double *B,int N_num)
{   
   double *pxsh;
   pxsh=new double[N_num*N_num];
   int i,k,j;
   for(i=0;i<N_num;i++)
   {
	   for(k=0;k<N_num;k++)
		   *(pxsh+i*N_num+k)=*(B+i*N_num+k);
   }

	for(i=0;i<N_num;i++)
	{
		for(k=0;k<N_num;k++)
			if(i==k)
				*(pI+i*N_num+k)=1.0;
			else
				*(pI+i*N_num+k)=0.0;
	}
	for(i=0;i<N_num-1;i++)
	{
		int line=i;//标志绝对值最大的元素所在地行
		double max=fabs(*(pxsh+i*N_num+i));
		for(j=i;j<N_num;j++)//求取每一列的绝对值最大的
		{
			if(fabs(*(pxsh+j*N_num+i))>max)
			{
				max=fabs(*(pxsh+j*N_num+i));
				line=j;
			}
			if(max==0)
			{
				printf("第%d列元素为0,方程组无解!\n",j+1);
				exit(-1);
			}
		}
		while(1)
		{
			if(line==i)
				break;
			double temp;
			for(int j=0;j<N_num;j++)//交换对角线行和绝对值最大的行
			{
				temp=*(pxsh+i*N_num+j);
				*(pxsh+i*N_num+j)=*(pxsh+line*N_num+j);
				*(pxsh+line*N_num+j)=temp;

				temp=*(pI+i*N_num+j);
				*(pI+i*N_num+j)=*(pI+line*N_num+j);
				*(pI+line*N_num+j)=temp;
			}
			break;
		}
		//将非对角线元素置为0
		for(j=0;j<N_num;j++)
		{
			if(i==j)
				continue;
			if(*(pxsh+j*N_num+i)!=0)
			{
				double value=*(pxsh+j*N_num+i)/(*(pxsh+i*N_num+i));
				for(int m=0;m<N_num;m++)
				{
					*(pxsh+j*N_num+m)-=*(pxsh+i*N_num+m)*value;
					*(pI+j*N_num+m)-=*(pI+i*N_num+m)*value;
					//*(pxsh+j*N_num+i)=0.0;
				}
			}
		}
	}
	if(*(pxsh+i*N_num+i)==0)
	{
		printf("最后一个元素为0!该方程组无解!\n");
		exit(-1);
	}
	//对最后一列进行操作
	for(int ii=0;ii<N_num-1;ii++)
	{
		if(*(pxsh+ii*N_num+N_num-1)!=0)
		{
            double value=*(pxsh+ii*N_num+N_num-1)/(*(pxsh+i*N_num+i));
			*(pxsh+ii*N_num+N_num-1)=0.0;
			for(int m=0;m<N_num;m++)
				*(pI+ii*N_num+m)-=*(pI+i*N_num+m)*value;			
		}
	}
	//以下将对角线元素变为1，从而得到方程的解
	for(i=0;i<N_num;i++)
	{
		double value=*(pxsh+i*N_num+i);
		for(int j=0;j<N_num;j++)
		{
			
			*(pxsh+i*N_num+j)/=value;
			*(pI+i*N_num+j)/=value;
		}
	}
	delete []pxsh;
	pxsh=NULL;
 }