计算方法工程数学第一节课Doolittle

Doolitte实验代码

#include <map>
#include <queue>
#include <string>
#include<iostream>
#include<stdio.h>
#include<string.h>
#include <algorithm>
#include <math.h>
typedef long long ll;
typedef unsigned long long ull;
using namespace std;
typedef pair<ll,ll> pii;
#define mem(a,x) memset(a,x,sizeof(a))
#define debug(x) cout << #x << ": " << x << endl;
#define rep(i,n) for(int i=0;i<(n);++i)
#define repi(i,a,b) for(int i=int(a);i<=(b);++i)
#define repr(i,b,a) for(int i=int(b);i>=(a);--i)
const int maxn=2e5+1010;
#define inf 0x3f3f3f3f
#define sf scanf
#define pf printf
const int mod=998244353;
const int MOD=10007;

inline int read() {
	int x=0;
	bool t=false;
	char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}

double l[10][10],u[10][10];
double b[maxn],x[maxn],y[maxn];
double a[10][10];
ll n;

double sum(int k,int j){
	double s=0;
	for(int i=1;i<=k-1;i++){
		s+=l[k][i]*u[i][j];
	}
	return s;
}
double sum2(int i,int k){
	double ss=0;
	for(int s=1;s<k;s++){
		ss+=l[i][s]*u[s][k];
	}
	return ss;
	
}
double sum3(int k){
	double r=0;
	for(int i=1;i<k;i++){
		r+=l[k][i]*y[i];
	}
	return r;	
}
double sum4(int k){
	double r=0;
	for(int i=k+1;i<=n;i++){
		r+=u[k][i]*x[i];
	}
	return r;	
}
int main() {
	cin>>n;
	for(int i=1; i<=n; i++) {
		for(int j=1; j<=n; j++) {
			cin>>a[i][j];
		}
	}
	for(int i=1;i<=n;i++) l[i][i]=1;
	for(int i=1;i<=n;i++) u[1][i]=a[1][i];
	
	for(int i=2;i<=n;i++) l[i][1]=a[i][1]/u[1][1];
	
	for(int k=2; k<=n; k++) {
		for(int i=k;i<=n;i++){
			u[k][i]=a[k][i]-sum(k,i);
		}
		for(int i=k+1;i<=n;i++){
			l[i][k]=(a[i][k]-sum2(i,k))/u[k][k];
		}
	}
	/*
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			cout<<u[i][j]<<"　"; 
		} 
		cout<<endl; 
		
	}
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			cout<<l[i][j]<<"　";
		} 
		cout<<endl; 		
	}
	*/
	for(int i=1;i<=n;i++) cin>>b[i];
	y[1]=b[1];
	for(int i=2;i<=n;i++){
		y[i]=b[i]-sum3(i);
	}
	///for(int i=1;i<=n;i++) cout<<y[i]<<" ";
	x[n]=y[n]/u[n][n];
	for(int i=n-1;i>=1;i--){
		x[i]=(y[i]-sum4(i))/u[i][i];
	}
	for(int i=1;i<=n;i++) cout<<x[i]<<" ";
	return 0;
}