[数值算法]列主元三角分解法

[数值算法]列主元三角分解法

                                                                                                  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

最后

以上就是欢呼冰淇淋为你收集整理的[数值算法]列主元三角分解法的全部内容，希望文章能够帮你解决[数值算法]列主元三角分解法所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。