//封装函数三voidSolverEqCholesky(double** A,double* b,int n,double* x,double** L){double omiga ,lamda;for(int i =1;i <= n;i++){
lamda =0;
omiga =0;for(int k =1;k < i;k++)
lamda += L[i][k]* L[i][k];
L[i][i]=sqrt(A[i][i]- lamda);for(int j = i +1;j <= n;j++){for(int k =1;k < i;k++)
omiga += L[k][i]* L[k][j];
L[j][i]=L[i][j]=(A[i][j]- omiga)/ L[i][i];
omiga =0;}}double y[100];for(int u =1;u <= n;u++){for(int v =1;v < u;v++){
b[u]-= L[u][v]* y[v];}
y[u]= b[u]/ L[u][u];}for(int u = n;u >0;u--){for(int v = n;v > u;v--){
y[u]-= L[u][v]* x[v];}
x[u]= y[u]/ L[u][u];}}
//开方可能影响精度并增加运算量#include<iostream>#include<string.h>#include<stdio.h>#include<vector>#include<math.h>#include"windows.h"usingnamespace std;constint N =1005;typedefdouble Type;
Type A[N][N], L[N][N];/** 分解A得到A = L * L^T */voidCholesky(Type A[][N], Type L[][N],int n){for(int k =0; k < n; k++){
Type sum =0;for(int i =0; i < k; i++)
sum += L[k][i]* L[k][i];
sum = A[k][k]- sum;
L[k][k]=sqrt(sum >0? sum :0);for(int i = k +1; i < n; i++){
sum =0;for(int j =0; j < k; j++)
sum += L[i][j]* L[k][j];
L[i][k]=(A[i][k]- sum)/ L[k][k];}for(int j =0; j < k; j++)
L[j][k]=0;}}/** 回带过程 */
vector<Type>Solve(Type L[][N], vector<Type> X,int n){/** LY = B => Y */for(int k =0; k < n; k++){for(int i =0; i < k; i++)
X[k]-= X[i]* L[k][i];
X[k]/= L[k][k];}/** L^TX = Y => X */for(int k = n -1; k >=0; k--){for(int i = k +1; i < n; i++)
X[k]-= X[i]* L[i][k];
X[k]/= L[k][k];}return X;}voidPrint(Type L[][N],const vector<Type> B,int n){
vector<Type> X =Solve(L, B, n);
vector<Type>::iterator it;int op;for(it = X.begin(),op=1; it != X.end(); op++,it++)printf("x%d=%.6f\n",op,*it);
cout<<endl;}intmain(){
LARGE_INTEGER lFrequency, lEndCount, lBeginCount;//LARGE_INTEGER为数据类型QueryPerformanceFrequency(&lFrequency);//获取CPU的时钟频率QueryPerformanceCounter(&lBeginCount);//获取CPU计数器数字,放在计时开头double CompuTime;memset(L,0,sizeof(L));int n=20;for(int i=0;i<n;i++){for(int j=0;j<n;j++){
A[i][j]=1.0/(double)(i+j+1);}}
vector<Type> B;for(int i =1; i <=n; i++){
Type y;
y=(double)i;
B.push_back(y);}Cholesky(A, L, n);Print(L, B, n);QueryPerformanceCounter(&lEndCount);//获得CPU计数器数字,放在计时结尾
CompuTime =(double)(lEndCount.QuadPart - lBeginCount.QuadPart)/(double)lFrequency.QuadPart;//得到运行时间,单位微秒
cout<<"Cholesky法计算时间为"<<CompuTime;return0;}
//改进版本#include<iostream>#include<string.h>#include<stdio.h>#include<vector>#include<math.h>#include"windows.h"usingnamespace std;constint N =1005;typedefdouble Type;
Type A[N][N], L[N][N], D[N][N];/** 分解A得到A = LDL^T */voidCholesky(Type A[][N], Type L[][N], Type D[][N],int n){for(int k =0; k < n; k++){for(int i =0; i < k; i++)
A[k][k]-= A[i][i]* A[k][i]* A[k][i];for(int j = k +1; j < n; j++){for(int i =0; i < k; i++)
A[j][k]-= A[j][i]* A[i][i]* A[k][i];
A[j][k]/= A[k][k];}}memset(L,0,sizeof(L));memset(D,0,sizeof(D));for(int i =0; i < n; i++){
D[i][i]= A[i][i];
L[i][i]=1;}for(int i =0; i < n; i++){for(int j =0; j < i; j++)
L[i][j]= A[i][j];}}voidTransposition(Type L[][N],int n){for(int i =0; i < n; i++){for(int j =0; j < i; j++)swap(L[i][j], L[j][i]);}}voidMulti(Type A[][N], Type B[][N],int n){
Type **C =new Type*[n];for(int i =0; i < n; i++)
C[i]=new Type[n];for(int i =0; i < n; i++){for(int j =0; j < n; j++){
C[i][j]=0;for(int k =0; k < n; k++)
C[i][j]+= A[i][k]* B[k][j];}}for(int i =0; i < n; i++){for(int j =0; j < n; j++)
B[i][j]= C[i][j];}for(int i =0; i < n; i++){delete[] C[i];
C[i]=NULL;}delete C;
C =NULL;}/** 回带过程 */
vector<Type>Solve(Type L[][N], Type D[][N], vector<Type> X,int n){/** LY = B => Y */for(int k =0; k < n; k++){for(int i =0; i < k; i++)
X[k]-= X[i]* L[k][i];
X[k]/= L[k][k];}/** DL^TX = Y => X */Transposition(L, n);Multi(D, L, n);for(int k = n -1; k >=0; k--){for(int i = k +1; i < n; i++)
X[k]-= X[i]* L[k][i];
X[k]/= L[k][k];}return X;}voidPrint(Type L[][N], Type D[][N],const vector<Type> B,int n){int op;
vector<Type> X =Solve(L, D, B, n);
vector<Type>::iterator it;for(it = X.begin(),op=1; it != X.end(); op++,it++)printf("x%d=%.6f\n",op,*it);
cout<<endl;}intmain(){
LARGE_INTEGER lFrequency, lEndCount, lBeginCount;//LARGE_INTEGER为数据类型QueryPerformanceFrequency(&lFrequency);//获取CPU的时钟频率QueryPerformanceCounter(&lBeginCount);//获取CPU计数器数字,放在计时开头double CompuTime;int n=10;memset(L,0,sizeof(L));for(int i=0;i<n;i++){for(int j=0;j<n;j++){
A[i][j]=1.0/(double)(i+j+1);}}
vector<Type> B;for(int i =0; i < n; i++){
Type y;
y=(double)i;
B.push_back(y);}Cholesky(A, L, D, n);Print(L, D, B, n);QueryPerformanceCounter(&lEndCount);//获得CPU计数器数字,放在计时结尾
CompuTime =(double)(lEndCount.QuadPart - lBeginCount.QuadPart)/(double)lFrequency.QuadPart;//得到运行时间,单位微秒
cout<<"Cholesky法计算时间为"<<CompuTime;return0;}