[数值算法]列主元三角分解法
By EmilMatthew 05/9/13
列主元三角分角法是对直接三角分解法的一种改进,主要目的和列主元高斯消元法一样,就是避免小数作为分母项.
有一些麻烦是的是要对原矩阵的增广矩阵部分进行LU分解,而且每次要对当前未处不景行作求列主行的下标,然后做变换.
其它的注意点和LU分解法其本一致,而最后求出U矩阵后,可以直接用回代法求出解x,而不必先求y再用Ux=y来求解向量x了.
其实,这里根本没有必要用到数学上推导时用的L和U两个三角矩阵,只要用一个增广矩阵就可以解决战斗了,不过我实现时还是用了老方法,以后有待改进.
/*
LUColMainMethod, coded by EmilMathew 05/9/13, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.
*/
void LUColMainMethod(Type** matrixArr,Type* bList,Type* xAnsList,int len)
{
/*Core think of this algorithm:
matrix_L*yAnsList=bList;
matrix_U*xAnsList=yAnsList;
*/
Type** matrix_L,** matrix_U;/*The core two matrix object of the LU method.*/
Type* yAnsList;
int i,j,k;/*iterator num*/
int i2,j2;
Type sum;/*the temp tween data*/
/*LU Col
Main Spec Data*/
int maxRowIndex;
Type** matrixA_B;
/*input assertion*/
assertF(matrixArr!=NULL,"in LUPationMethod,matrixArr is NULL\n");
assertF(bList!=NULL,"in LUPationMethod,bList is NULL\n");
assertF(xAnsList!=NULL,"in LUPationMethod,xAnsList is NULL\n");
/*memory apply*/
matrix_L=(Type**)malloc(sizeof(Type*)*len);
twoDArrMemApply(&matrix_U,len,len+1);
for(i=0;i<len;i++)
{
matrix_L[i]=(Type*)malloc(sizeof(Type)*len);
// matrix_U[i]=(Type*)malloc(sizeof(Type)*len);
}
yAnsList=(Type*)malloc(sizeof(Type)*len);
/*==end of memory apply==*/
/*---assertion after apply the mem---*/
assertF(matrix_L!=NULL,"in LUPationMethod,matrix_L is NULL\n");
assertF(matrix_U!=NULL,"in LUPationMethod,matrix_U is NULL\n");
assertF(yAnsList!=NULL,"in LUPationMethod,yAnsList is NULL\n");
/*---end of assertion---*/
// printf("arr mem applyed\n");
/*----data initialization----*/
for(i=0;i<len;i++)
{
for(j=0;j<len;j++)
{
matrix_L[i][j]=0;
matrix_U[i][j]=0;
}
matrix_U[i][j]=0;
}
for(i=0;i<len;i++)
matrix_L[i][i]=1;
/*matrix A_B prepare*/
twoDArrMemApply(&matrixA_B,len,len+1);
matrixCopy(matrixArr,matrixA_B,len,len);
for(i=0;i<len;i++)
matrixA_B[i][len]=bList[i];
printf("show copied matrix:\n");
show2DArrFloat(matrixA_B,len,len+1);
for(i=0;i<len;i++)
{
maxRowIndex=getMaxRowIndexLU(matrixA_B,i,matrix_L,matrix_U,len);
if(i!=maxRowIndex)
{
swapTwoRow(matrixA_B,i,maxRowIndex,len+1);
swapTwoRow(matrix_L,i,maxRowIndex,len);
swapTwoRow(matrix_U,i,maxRowIndex,len);
}
/*matrixU make*/
for(j=i;j<len;j++)
{
sum=0;
for(k=0;k<i;k++)
sum+=matrix_L[i][k]*matrix_U[k][j];
matrix_U[i][j]=matrixA_B[i][j]-sum;
}
matrix_U[i][j]=matrixA_B[i][j]-sumArr1_IKByArr2_KJ(matrix_L,matrix_U,i,j,0,i-1);
matrixA_B[i][j]=matrix_U[i][j];
/*matrixL make*/
if(i<len-1)
{
j2=i;
for(i2=j2+1;i2<len;i2++)
{
sum=0;
for(k=0;k<j2;k++)
sum+=matrix_L[i2][k]*matrix_U[k][j2];
matrix_L[i2][j2]=(matrixA_B[i2][j2]-sum)/matrix_U[j2][j2];
}
}
}
//copyBack
for(i=0;i<len;i++)
bList[i]=matrixA_B[i][len];
printf("show bList\n");
showArrListFloat(bList,0,len);
/*Adjusting of matrix_L*/
for(i=0;i<len;i++)
{
for(j=i;j<len;j++)
{
if(i==j)
matrix_L[i][i]=1;
else
matrix_L[i][j]=0;
}
}
printf("matrix L\n");
show2DArrFloat(matrix_L,len,len);
printf("matrix U\n");
show2DArrFloat(matrix_U,len,len);
for(i=len-1;i>=0;i--)
{
xAnsList[i]=(bList[i]-sumArr_JKByList_K(matrix_U,xAnsList,i,i+1,len-1))/matrix_U[i][i];
}
/*
downTriangleMethod(matrix_L,bList,yAnsList,len);
upTriangleMethod(matrix_U,yAnsList,xAnsList,len);
*/
/*memory free*/
for(i=0;i<len;i++)
{
free(matrix_L[i]);
free(matrix_U[i]);
}
free(matrix_L);
free(matrix_U);
free(yAnsList);
free(matrixA_B);
}
//辅助函数,求列主元的行位置.
int getMaxRowIndexLU(Type** inMatrixArr,int currentI,Type** matrix_L,Type** matrix_U,int size)
{
int i,maxRowIndex;
Type tmpSMax,tmpCompare;
assertF(inMatrixArr!=NULL,"in getMaxRowIndexLU,inMatrixArr is NULL\n");
assertF(matrix_L!=NULL,"in getMaxRowIndexLU,matrix_L is NULL\n");
assertF(matrix_U!=NULL,"in getMaxRowIndexLU,matrix_U is NULL\n");
tmpSMax=0;
// maxRowIndex=currentI;
for(i=currentI;i<size;i++)
{
tmpCompare=(float)fabs(inMatrixArr[i][currentI]-sumArr1_IKByArr2_KJ(matrix_L,matrix_U,i,currentI,0,currentI-1));
if(tmpCompare>tmpSMax)
{
tmpSMax=tmpCompare;
maxRowIndex=i;
}
}
return maxRowIndex;
}
我所写的程序的源码下载:
http://emilmatthew.51.net/downloads/luColMain.rar
//如果上面这个链接无法响应下载（有可能是被网站给屏蔽掉了），则可使用下载工具（如迅雷等）下载
测试结果:
test1:
input:
3;
4,-2,-4,10;
-2,17,10,3;
-4,10,9,-7;
output:
after the LUPationMethod:the ans x rows is:(from x0 to xn-1)
2.00000 1.00000 -1.00000
test2:
input:
3;
12,-3,3,15;
18,-3,1,15;
-1,2,1,6;
after the LUPationMethod:the ans x rows is:(from x0 to xn-1)
1.00000 2.00000 3.00000
