Mariusz:
#include <iostream>
#include <complex>
#include <conio.h>
using namespace std;
complex<double> det(complex<double>**,int,double);
void swapcols(int,int,int,complex<double>**,complex<double>**);
void cramer(complex<double>**,int,double,complex<double>*);
void gauss(complex<double>**,int,double,complex<double>*);
void inverse(complex<double>**,int,double,int&);
void multiply(complex<double>**,complex<double>**,complex<double>**,int,int,int);
void ludcmp(complex<double>**,int,int*,complex<double>&);
void lubksb(complex<double>**,int,int*,complex<double>*);
int main(){
char esc;
int i,j,n;
int err;
int c;
int* indx;
complex<double> **A,**B,**C,*X;
double eps=1e−12;
complex<double> d;
//clrscr();
do{
cout<<"Rozwiazywanie ukladow rownan liniowych "<<endl;
cout<<"1 Metoda Cramera "<<endl;
cout<<"2 Metoda Gaussa "<<endl;
cout<<"3 Metoda X=A
(−1)*B "<<endl;
cout<<"4 Metoda rozkladu LU "<<endl;
cout<<endl;
cin>>c;
cout<<endl;
cout.precision(12);
switch(c){
case 1:{
cout<<"Podaj n=";
cin>>n;
A=new complex<double>*[n+2];
for(i=0;i<=n+1;i++)
A[i]=new complex<double>[n+2];
X=new complex<double>[n+1];
for(i=1;i<=n;i++)
for(j=1;j<=n+1;j++){
cout<<"a["<<i<<"]["<<j<<"]=";
cin>>A[i][j];
}
cramer(A,n,eps,X);
if(X[0]!=0.0)
for(i=1;i<=n;i++)
cout<<"X["<<i<<"]="<<X[i]<<endl;
cout<<endl;
for(i=0;i<=n+1;i++)
delete[] A[i];
delete[] A;
delete[] X;
break;
}
case 2:{
cout<<"Podaj n=";
cin>>n;
A=new complex<double>*[n+2];
for(i=0;i<=n+1;i++)
A[i]=new complex<double>[n+2];
X=new complex<double>[n+1];
for(i=1;i<=n;i++)
for(j=1;j<=n+1;j++){
cout<<"A["<<i<<"]["<<j<<"]=";
cin>>A[i][j];
}
gauss(A,n,eps,X);
for(i=1;i<=n;i++)
cout<<"X["<<i<<"]="<<X[i]<<endl;
cout<<endl;
for(i=0;i<=n+1;i++)
delete[] A[i];
delete[] A;
delete[] X;
break;
}
case 3:{
cout<<"Podaj n=";
cin>>n;
A=new complex<double>*[n+1];
for(i=0;i<=n;i++)
A[i]=new complex<double>[n+1];
B=new complex<double>*[2];
for(i=0;i<=n;i++)
B[i]=new complex<double>[2];
C=new complex<double>*[2];
for(i=0;i<=n;i++)
C[i]=new complex<double>[2];
for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
cout<<"A["<<i<<"]["<<j<<"]=";
cin>>A[i][j];
}
cout<<"B["<<i<<"]=";
cin>>B[i][1];
}
inverse(A,n,eps,err);
multiply(C,A,B,n,1,n);
if(err==0)
for(i=1;i<=n;i++)
cout<<"X["<<i<<"]="<<C[i][1]<<endl;
for(i=0;i<=n;i++)
delete[] A[i];
delete[] A;
for(i=0;i<=1;i++)
delete[] B[i];
delete[] B;
for(i=0;i<=1;i++)
delete[] C[i];
delete[] C;
break;
}
case 4:{
cout<<"Podaj n=";
cin>>n;
A=new complex<double>*[n+2];
for(i=0;i<=n+1;i++)
A[i]=new complex<double>[n+2];
X=new complex<double>[n+1];
indx=new int[n+1];
for(i=1;i<=n;i++)
for(j=1;j<=n+1;j++){
cout<<"A["<<i<<"]["<<j<<"]=";
cin>>A[i][j];
}
ludcmp(A,n,indx,d);
for(i=1;i<=n;i++)
X[i]=A[i][n+1];
if(d!=0.0){
lubksb(A,n,indx,X);
for(i=1;i<=n;i++)
cout<<"X["<<i<<"]="<<X[i]<<endl;
cout<<endl;
}
for(i=0;i<=n+1;i++)
delete[] A[i];
delete[] A;
delete[] X;
delete[] indx;
break;
}
default: cout<<"Wprowadz liczbe od 1 do 4 "<<endl;
}
esc=getch();
}
while (esc!=27);
return 0;
}
void inverse(complex<double>** A,int n,double eps,int& blad){
int i,j,k,l;
complex<double> maxA,d,e;
int* M;
int b;
M=new int [n+1];
b=0;
for(i=1;i<=n;i++){
maxA=0.0;
for(j=i;j<=n;j++){
d=A[j][i];
if(abs(maxA)<abs(d)){
maxA=d;
k=j;
}
}
if(abs(maxA)<eps) b=1;
if(!b){
M[i]=k;
A[k][i]=1.0;
for(j=1;j<=n;j++){
d=complex<double>(A[k][j]/maxA);
A[k][j]=A[i][j];
A[i][j]=d;
}
for(j=1;j<=n;j++)
if(j!=i){
d=A[j][i];
A[j][i]=0.0;
for(l=1;l<=n;l++){
e=d*A[i][l];
A[j][l]−=e;
}
}
}
}
if(!b){
for(i=n;i>=1;i−−){
k=M[i];
if(k!=i)
for(j=1;j<=n;j++){
d=A[j][i];
A[j][i]=A[j][k];
A[j][k]=d;
}
}
}
blad=b;
delete[] M;
}
void multiply(complex<double>** D,complex<double>** A,complex<double>** B,int m,int p,int n){
int i,j,k;
complex<double> s;
for(i=1;i<=m;i++){
for(k=1;k<=p;k++) {
s=0.0;
for(j=1;j<=n;j++)
s+=A[i][j]*B[j][k];
D[i][k]=s;
}
}
}
void gauss(complex<double>** a,int n,double eps,complex<double>* b){
int i,j,k,wm;
complex<double> p,w;
double max;
w=1.0;
for (i=1;i<n;i++){
max=abs(a[i][i]);
wm=i;
for(j=i+1;j<=n;j++)
if(abs(a[j][i])>max){
max=abs(a[j][i]);
wm=j;
}
if(wm>i) {
for(k=i;k<=n+1;k++){
p=a[i][k];
a[i][k]=a[wm][k];
a[wm][k]=p;
}
w*=−1.0;
}
if (abs(a[i][i])<eps) w=0;
if (w!=0.0)
for(j=i+1;j<=n;j++){
p=complex<double>(a[j][i]/a[i][i]);
for(k=n+1;k>i;k−−)
a[j][k]−=p*a[i][k];
}
w*=a[i][i];
}
w*=a[n][n];
if(abs(w)<eps) w=0.0;
b[0]=w;
if(b[0]!=0.0)
for(i=n;i>=1;i−−){
p=a[i][n+1];
for(j=i+1;j<=n;j++)
p−=a[i][j]*b[j];
b[i]=complex<double>(p/a[i][i]);
}
}
complex<double> det(complex<double>** a,int n,double eps){
int i,j,k,wm;
complex<double> p,w;
double max;
w=1.0;
for (i=1;i<n;i++){
max=abs(a[i][i]);
wm=i;
for(j=i+1;j<=n;j++)
if(abs(a[j][i])>max){
max=abs(a[j][i]);
wm=j;
}
if(wm>i) {
for(k=i;k<=n+1;k++){
p=a[i][k];
a[i][k]=a[wm][k];
a[wm][k]=p;
}
w*=−1.0;
}
if (abs(a[i][i])<eps) w=0;
if (w!=0.0)
for(j=i+1;j<=n;j++){
p=complex<double>(a[j][i]/a[i][i]);
for(k=n+1;k>i;k−−)
a[j][k]−=p*a[i][k];
}
w*=a[i][i];
}
w*=a[n][n];
return w;
}
void swapcols(int n,int p,int q,complex<double>** a,complex<double>** b){
int i,j;
complex<double> t;
for(i=1;i<=n;i++)
for(j=1;j<=n+1;j++)
b[i][j]=a[i][j];
for(i=1;i<=n;i++){
t=b[i][p];
b[i][p]=b[i][q];
b[i][q]=t;
}
}
void cramer(complex<double>** A,int n,double eps,complex<double>* X){
int i;
complex<double>** B;
B=new complex<double>*[n+2];
for(i=0;i<=n+1;i++)
B[i]=new complex<double>[n+2];
for(i=1;i<=n;i++){
swapcols(n,i,n+1,A,B);
X[i]=det(B,n,eps);
}
X[0]=det(A,n,eps);
if(X[0]!=0.0)
for(i=1;i<=n;i++)
X[i]=complex<double>(X[i]/X[0]);
for(i=0;i<=n+1;i++)
delete[] B[i];
delete[] B;
}
void ludcmp(complex<double>**a,int n,int* indx,complex<double>&d)
{
int i,imax,j,k;
complex<double> big,dum,sum,temp;
complex<double>* vv;
double tiny=1e−20;
vv=new complex<double>[n+1];
d=1.0;
for(i=1;i<=n;i++){
big=0.0;
for(j=1;j<=n;j++)
if(abs((temp=abs(a[i][j])))>abs(big)) big=temp;
if(big==0.0) {
d=0.0;
return;
}
vv[i]=complex<double>(1.0/big);
}
for(j=1;j<=n;j++){
for(i=1;i<j;i++){
sum=a[i][j];
for(k=1;k<i;k++) sum−=a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for(i=j;i<=n;i++){
sum=a[i][j];
for(k=1;k<j;k++)
sum−=a[i][k]*a[k][j];
a[i][j]=sum;
if(abs((dum=vv[i]*abs(sum)))>=abs(big)){
big=dum;
imax=i;
}
}
if(j!=imax){
for(k=1;k<=n;k++){
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
d=−d;
vv[imax]=vv[j];
}
indx[j]=imax;
if(a[j][j]==0.0) a[j][j]=tiny;
if(j!=n){
dum=1.0/(a[j][j]);
for(i=j+1;i<=n;i++) a[i][j]*=dum;
}
}
delete[] vv;
}
void lubksb(complex<double>** a,int n,int*indx,complex<double> *b)
{
int i,ip,j,ii=0;
complex<double> sum;
for(i=1;i<=n;i++) {
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
if(ii)
for(j=ii;j<=i−1;j++) sum−=a[i][j]*b[j];
else if(abs(sum)) ii=i;
b[i]=sum;
}
for(i=n;i>=1;i−−){
sum=b[i];
for(j=i+1;j<=n;j++) sum−= a[i][j]*b[j];
b[i]=complex<double>(sum/a[i][i]);
}
}
Spróbuj coś z tego wyciągnąć