[数值算法]线性方程组的求解---平方根法及改进平方根法
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
