[数值算法]线性方程组的求解---平方根法及改进平方根法

                                                                                                         By EmilMatthew

                                                                                                                       05/9/10

平方根法主要用于求解对称正定矩阵方程：

首先要提到的是有个关于正定矩阵的定理，说的是若A为n阶地称正定矩阵，则存在一个实的非奇异下三角矩阵L，使A=LL’（L’表示L的对称矩阵）

根据这个前提，在结合前面的LU分解算法，便有了这里的平方根算法：

 

平方根方法:

/*

squareRootMethod, coded by EmilMathew 05/9/10, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.

*/

void squareRootMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,int size)

{

       /*Maths Reason:

       L*U=A*x=b

       When you meet a duiCheng and zheng ding matrix:

       it could be pation like this:

       l11                           l11  l21  ...  ln1

       l21  l22               *          l22  ...  ln2

       ....                                           ...  ...

       ln1  ln2  ...  lnn                           lnn

      

       i=j:

       lij=sqrt(aij-sigma_k1...j-1(ljk^2))

      

       i>j:

       lij=(aij-sigma_k1...j-1(lik*ljk))/ljj

      

       the steps below is very easy :

       L*y=b;

       U*x=y;

      

       Enjoy!:)

                                                                      by EmilMatthew

                                                                             05/9/10.

       */

       Type** l_Matrix,* yAnsList;

       Type tmpData;

       int i,j;

      

       /*pointer data assertion*/

       assertF(inMatrixArr!=NULL,"in squareRootMethod,matrixArr is NULL/n");

       assertF(bList!=NULL,"in squareRootMethod,bList is NULL/n");

       assertF(xAnsList!=NULL,"in squareRootMethod,xAnsList is NULL/n");   

      

       /*correct pass in matrix assertion*/

       assertF(duiChengMatrixCheck(inMatrixArr,size),"in squareRootMethod,the pass in matrix is not dui cheng/n");

      

       /*Mem Apply*/            

       listArrMemApply(&yAnsList,size);

       twoDArrMemApply(&l_Matrix,size,size);

      

       assertF(l_Matrix!=NULL,"in squareRootMethod,l_Matrix is null/n");

       assertF(yAnsList!=NULL,"in squareRootMethod,yAnsList is null/n");

      

       /*Core Program*/

       for(i=0;i<size;i++)

              for(j=0;j<=i;j++)

              {

                     if(i==j)

                     {

                            tmpData=sumSomeRowPower(l_Matrix,j,0,j-1,2);

                     //     printf("tmpData:%f/n",tmpData);

                            l_Matrix[i][j]=(float)sqrt(inMatrixArr[i][i]-tmpData);

                     }

                     else

                     {

       l_Matrix[i][j]=(inMatrixArr[i][j]-sumTwoRowBy(l_Matrix,i,j,0,j-1))/l_Matrix[j][j];

                     }

              }

for(i=0;i<size;i++)

yAnsList[i]=(bList[i]-sumArr_JKByList_K(l_Matrix,yAnsList,i,0,i-1))/l_Matrix[i][i];

 

for(i=size-1;i>=0;i--)

xAnsList[i]=(yAnsList[i]-sumArr_KJByList_K(l_Matrix,xAnsList,i,i+1,size-1))/l_Matrix[i][i];

      

       twoDArrMemFree(&l_Matrix,size);

       free(yAnsList);

}

 

平方根算法的计算量约为普通三角分解算法的一半，但由于这里要用到开平方，效率不是很高，所以，便有了为效率而存在的改进版平方根算法:

改进的平方根算法:

/*

enhancedSquareRootMethod, coded by EmilMathew 05/9/10, you can modify and use these code as you wish , but there is no guarantee that it can fit all you need.

*/

 

如下:

void enhancedSquareRootMethod(Type** inMatrixArr,Type* bList,Type* xAnsList,int size)

{

       Type** l_Matrix,** t_Matrix;

       Type* yAnsList,* dList;

 

       int i,j;

      

       /*pointer data assertion*/

       assertF(inMatrixArr!=NULL,"in enhancedSquareRootMethod,matrixArr is NULL/n");

       assertF(bList!=NULL,"in enhancedSquareRootMethod,bList is NULL/n");

       assertF(xAnsList!=NULL,"in enhancedSquareRootMethod,xAnsList is NULL/n");   

      

       /*correct pass in matrix assertion*/

       assertF(duiChengMatrixCheck(inMatrixArr,size),"in enhancedSquareRootMethod,the pass in matrix is not dui cheng/n");

      

       /*Mem Apply*/            

       listArrMemApply(&yAnsList,size);

       listArrMemApply(&dList,size);

       twoDArrMemApply(&l_Matrix,size,size);

       twoDArrMemApply(&t_Matrix,size,size);

      

       assertF(t_Matrix!=NULL,"in enhancedSquareRootMethod,t_Matrix is null/n");

       assertF(l_Matrix!=NULL,"in enhancedSquareRootMethod,l_Matrix is null/n");

       assertF(yAnsList!=NULL,"in enhancedSquareRootMethod,yAnsList is null/n");

      

       for(i=0;i<size;i++)

              l_Matrix[i][i]=1;

      

       for(j=0;j<size;j++)

       {

              dList[j]=inMatrixArr[j][j]-sumArr1_IKByArr2_JK(t_Matrix,l_Matrix,j,j,0,j-1);

              for(i=j+1;i<size;i++)

              {

                     t_Matrix[i][j]=inMatrixArr[i][j]-sumArr1_IKByArr2_JK(inMatrixArr,l_Matrix,i,j,0,j-1);

                     l_Matrix[i][j]=t_Matrix[i][j]/dList[j];

              }

       }

      

      

       for(i=0;i<size;i++)

              yAnsList[i]=bList[i]-sumArr_JKByList_K(l_Matrix,yAnsList,i,0,i-1);

 

for(i=size-1;i>=0;i--)

xAnsList[i]=yAnsList[i]/dList[i]-sumArr_KJByList_K(l_Matrix,xAnsList,i,i+1,size-1);

      

       /*mem free*/

       twoDArrMemFree(&t_Matrix,size);

       twoDArrMemFree(&l_Matrix,size);

       free(yAnsList);

       free(dList);

}

 

测试程序:

/*Square Method  Algorithm test program*/

#include "Global.h"

#include "Ulti.h"

#include "MyAssert.h"

#include "Matrix.h"

#include <time.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

 

 

char *inFileName="inputData.txt";

/*

       input data specification

       len,

       a00,a01,...,a0n-1,b0;

       .....

       an-10,an-11,...,an-1n-1,bn-1;

*/

 

 

char *outFileName="outputData.txt";

#define DEBUG 1

 

void main(int argc,char* argv[])

{

       FILE *inputFile;/*input file*/

       FILE *outputFile;/*output file*/

 

       double startTime,endTime,tweenTime;/*time callopsed info*/

      

       /*The read in data*/

       int len,methodIndex;

       Type** matrixArr;

       Type* bList,* xAnsList;

 

      

       int i,j;/*iterator index*/

      

       /*input file open*/

       if(argc>1)strcpy(inFileName,argv[1]);

       assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");

       printf("input file open success/n");

      

       /*outpout file open*/

       if(argc>2)strcpy(outFileName,argv[2]);

       assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");

       printf("output file open success/n");    

      

       fscanf(inputFile,"size=%d;/r/n",&len);

       fscanf(inputFile,"method=%d;/r/n",&methodIndex);

      

       /*Memory apply*/

       matrixArr=(Type**)malloc(sizeof(Type*)*len);

              for(i=0;i<len;i++)

                     matrixArr[i]=(Type*)malloc(sizeof(Type)*len);

      

       bList=(Type*)malloc(sizeof(Type)*len);

       xAnsList=(Type*)malloc(sizeof(Type)*len);

             

       /*Read info data*/

       for(i=0;i<len;i++)

       {

              for(j=0;j<len;j++)

                     fscanf(inputFile,"%f,",&matrixArr[i][j]);

              fscanf(inputFile,"%f;",&bList[i]);

       }

      

      

       /*Check the input data*/

       showArrListFloat(bList,0,len);

       show2DArrFloat(matrixArr,len,len);

      

#if  DEBUG

       printf("/n*******start of test program******/n");

       printf("now is runnig,please wait.../n");

       startTime=(double)clock()/(double)CLOCKS_PER_SEC;

       /******************Core program code*************/

              switch(methodIndex)

              {

              case 1:

                            enhancedSquareRootMethod(matrixArr,bList,xAnsList,len);

                            printf("after the enhancedSquareRootMethod:the ans x rows is:/n");

fprintf(outputFile,"after the enhancedSquareRootMethod:the ans x rows is:(from x0 to xn-1)/r/n");

                            break;

 

              case 2:

                    

                            squareRootMethod(matrixArr,bList,xAnsList,len);

                            printf("after the SquartRootPationMethod:the ans x rows is:/n");

fprintf(outputFile,"after the SquartRootPationMethod:the ans x rows is:(from x0 to xn-1)/r/n");

                            break;

                    

              default:

                            printf("input method index error/n");

                            break;

              }

              showArrListFloat(xAnsList,0,len);

              outputListArrFloat(xAnsList,0,len,outputFile);

       /******************End of Core program**********/

       endTime=(double)clock()/(double)CLOCKS_PER_SEC;

       tweenTime=endTime-startTime;/*Get the time collapsed*/

       /*Time collapsed output*/

       printf("the collapsed time in this algorithm implement is:%f/n",tweenTime);

       fprintf(outputFile,"the collapsed time in this algorithm implement is:%f/r/n",tweenTime); 

       printf("/n*******end of test program******/n");

#endif

       for(i=0;i<len;i++)

              free(matrixArr[i]);

       free(matrixArr);

      

       free(xAnsList);

       free(bList);

      

       printf("program end successfully,/n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer./n");

       getchar();/*Screen Delay Control*/

       return;

}

 

测试结果:

平方根法:

test1:

输入:

size=3;

5,-4,1,1;

-4,6,-4,2;

1,-4,6,-3;

输出:

after the SquartRootPationMethod:the ans x rows is:(from x0 to xn-1)

   1.00000    1.00000    0.00000

 

改进平方根法:

test1:

输入:

size=3;

method=1;

5,-4,1,1;

-4,6,-4,2;

1,-4,6,-3;

 

输出:

after the enhancedSquareRootMethod:the ans x rows is:(from x0 to xn-1)

   1.00000    1.00000    0.00000

 

Test2:

after the enhancedSquareRootMethod:the ans x rows is:(from x0 to xn-1)

   2.00000    1.00000   -1.00000

 

网上关于这个主题的相关参考:

http://jpkc.ecnu.edu.cn/gdds/xsxz/ZhangGengYun.htm

http://www.ascc.net/pd-man/linpack/node14.html