Cholesky法（汇总）

爱被打了一巴掌 2022-12-28 00:48 117阅读 0赞

//封装函数一
    double **two_array_malloc(int m, int n) 
    //二维数组分配 
    {  
    	double **a; 
    	int i, j; 
    	a = (double **)malloc(m * sizeof(double *)); 
    	for (i = 0; i < m; i++) 
    	{  
    		a[i] = (double *)malloc(n * sizeof(double)); 
    		for (j = 0; j < n; j++) 
    		{  
    			a[i][j] = 0;
    		} 
    	}
    	return a; 
    }
    void two_array_clear(double **a, int m) 
    //二维数组释放 
    {  
    	int i; 
    	for (i = 0; i < m; i++) 
    	{  
    		delete[] a[i]; 
    	}
    	delete[] a; 
    }
    void SolverEqCholesky(double **A, double *b, int n, double *x, double **L) 
    {  
    	int i, j, k; 
    	double ** a; 
    	double temp; 
    	a = two_array_malloc(n, n + 1); 
    	for (i = 0; i < n; i++) 
    	{  
    		for (j = 0; j < n; j++) 
    		{ 	 
    			a[i][j] = A[i][j]; 
    		}
    		a[i][n] = b[i]; 
    	}
    	//分解过程 
    	for (k = 0; k < n; k++) 
    	{  
    		temp = 0; 
    		for (i = 0; i < k; i++) 
    		{  
    			temp = temp + a[i][k] * a[i][k]; 
    		}
    		a[k][k] = sqrt(a[k][k] - temp); 
    		for (i = k + 1; i < n + 1; i++) 
    		{  
    			temp = 0; 
    			for (j = 0; j < k; j++) 
    			{  
    				temp = temp + a[j][k] * a[j][i];
    			}
    			a[k][i] = (a[k][i] - temp) / a[k][k]; 
    		} 
    	}
    	//回代过程 
    	a[n - 1][n] = a[n - 1][n] / a[n - 1][n - 1]; 
    	for (i = 0; i < n; i++) 
    	{  
    		for (j = 0; j < n; j++) 
    		{  
    			L[j][i]= a[i][j]; 
    		} 
    	}
    	for (k = n - 2; k >= 0; k--) 
    	{  
    		temp = 0; 
    		for (i = k + 1; i < n; i++) 
    		{  
    			temp = temp + a[k][i] * a[i][n]; 
    		}
    		a[k][n] = (a[k][n] - temp) / a[k][k]; 
    	}
    	//输出解 
    	for (i = 0; i < n; i++) 
    	{  
    		x[i] = a[i][n]; 
    	}
    	two_array_clear(a, n); 
    }

//封装函数二
    //较大阶数时可能报错
    #include<bits/stdc++.h>
    #include "windows.h" 
    using namespace std;
    int SolverEqCholesky(double **A,double *b,int n,double eps,double *x,double **L)
    { 
    	int m,i,j,k;
    	double **a =new double *[n];
    	double *btemp= new double[n];
    	double s;
    	for(i=0;i<n;i++)
    	{ 
    		a[i]=new double[n+1];
    	}
    	m=n;
    	//赋值
    	for(i=0;i<n;i++)
    	{ 
    		for(j=0;j<n;j++)
    		{ 
    			a[i][j]=A[i][j];
    		}
    		a[i][n]=b[i];
    	}  
    	if(m<=0||n<=0)
    	{ 
    		cout<<"Input Mistake"<<endl;
    	}
    	for(k=0;k<n;k++)
    	{ 
    		for(i=k;i<n;i++)
    		{ 
    			for(j=0;j<k;j++)
    			{ 
    				a[i][k]=a[i][k]-a[i][j]*a[k][j];
    			}
    		}
    		if(fabs(a[k][k])<eps)
    		{ 
    			printf("%f",a[k][k]);
    			return 1;
    		}
    		a[k][k]=sqrt(a[k][k]);
    		for(i=k+1;i<n;i++)
    		{ 
    			a[i][k]=a[i][k]/a[k][k];
    		}
    	} 
    	for(i=0;i<n;i++)
    	{ 
    		for(j=0;j<i;j++)
    		{ 
    			a[i][n]=a[i][n]-a[i][j]*a[j][n];
    		}
    		if(fabs(a[i][i])<eps)
    		{ 
    			return 1;
    		}
    		else
    		{ 
    			a[i][n]=a[i][n]/a[i][i];
    		}
    	}
    	for(i=n-1;i>=0;i--)
    	{ 
    		for(j=i+1;j<n;j++)
    		{ 
    			a[i][n]=a[i][n]-a[j][i]*a[j][n];
    		}
    		a[i][n]=a[i][n]/a[i][i];
    	}
    	for(i=0;i<m;i++)
    	{ 
    		x[i]=a[i][n];
    	}
    	for(i=0;i<n;i++)
    	{ 
    		delete[]a[i];
    	}
    	delete[]a;
    	delete[]btemp;
    }
    int main()
    { 
    	LARGE_INTEGER lFrequency, lEndCount, lBeginCount;//LARGE_INTEGER为数据类型
    	QueryPerformanceFrequency(&lFrequency);//获取CPU的时钟频率
    	QueryPerformanceCounter(&lBeginCount);//获取CPU计数器数字，放在计时开头
    	double CompuTime;
    	double **A,*b,*x,**L;
    	int kk=2;
    	A=new double *[kk];
    	b=new double[kk];
    	x=new double[kk];
    	L=new double *[kk];
    	for(int i=0;i<kk;i++)
    	{ 
    		A[i]=new double[kk];
    		L[i]=new double[kk];
    	}
    	int Q_num=1;
    	//test
    	//自己决定赋值或输入
    	int n=2;
    	b[0]=2;b[1]=2;
    	A[0][0]=1.001;A[0][1]=0.999;
    	A[1][0]=0.999;A[1][1]=1.001;
    	SolverEqCholesky(A,b,n,1.0e-10,x,L);
    	for(int i=0;i<n;i++)
    	{ 
    		cout<<"x"<<i+1<<"=";
    		printf("%f\n",x[i]);
    	}
    	QueryPerformanceCounter(&lEndCount);//获得CPU计数器数字，放在计时结尾
    	CompuTime = (double)(lEndCount.QuadPart - lBeginCount.QuadPart) / (double)lFrequency.QuadPart;//得到运行时间，单位微秒
    	cout<<"Cholesky法计算时间为"<<CompuTime; 
    }

//封装函数三
    void SolverEqCholesky(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" 
    using namespace std;
    const int N = 1005;
    typedef double Type;
    Type A[N][N], L[N][N];
    /** 分解A得到A = L * L^T */
    void Cholesky(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;
    }
     
    void Print(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;
    }
    int main()
    { 
        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; 
        return 0;
    }

//改进版本
    #include <iostream>
    #include <string.h>
    #include <stdio.h>
    #include <vector>
    #include <math.h>
    #include "windows.h"
    using namespace std;
    const int N = 1005;
    typedef double Type;
    Type A[N][N], L[N][N], D[N][N];
    /** 分解A得到A = LDL^T */
    void Cholesky(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];
        }
    }
    void Transposition(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]);
        }
    }
    void Multi(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;
    }
    void Print(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;
    }
    int main()
    { 
    	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; 
        return 0;
    }