ldlt.cpp

结构矩阵分析

//采用一维数组时p=i*(2*n-i+1)/2+j-i;
#include<iostream.h>
#include"error.h"
int LDLT(double *A,int n,double *B,int l);
//A传入左端矩阵，n为是方程维数，B是右端矩阵，l工况数
error LDLTsolve(double *A,int n,double *B,int l)
{
		int i,j;
		if(LDLT(A,n,B,l)){
			for(i=0;i<n;i++){
				cout<<"x("<<i+1<<")="<<B[i];
				for(j=1;j<l;j++){
					cout<<","<<B[i+j*n];
				}
				cout<<endl;
			}
			return yes;
		}
		else return zero;
}
int LDLT(double *A,int n,double *B,int l){
	double *a;int i,j,k,p,q,t,m;
	a=new double[n];
	for(k=0;k<n;k++){
		p=k*(2*n-k+1)/2;
		for(m=0;m<k;m++){
			q=m*(2*n-m+1)/2;
			a[m]=A[q+k-m]*A[q];
			A[p]=A[p]-a[m]*A[q+k-m];
			if(A[p]==0)return 0;
		}
		for(j=k+1;j<n;j++){
			t=p+j-k;
			for(m=0;m<k;m++){
				q=m*(2*n-m+1)/2;
				A[t]=A[t]-A[q+j-m]*a[m];
			}
			A[t]=A[t]/A[p];
		}
	}
	for(k=0;k<l;k++){
		p=k*n;
		for(i=1;i<n;i++){
			for(j=0;j<i;j++)
				B[i+p]=B[i+p]-A[j*(2*n-j+1)/2+i-j]*B[j+p];
		}
		for(i=0;i<n;i++)
			B[i+p]=B[i+p]/A[i*(2*n-i+1)/2];
		for(i=n-2;i>-1;i--){
			for(j=n-1;j>i;j--)
				B[i+p]=B[i+p]-A[i*(2*n-i+1)/2+j-i]*B[j+p];
		}
	}
	delete[] a;
	return 1;
}