#include<iostream>
using namespace std;

#include<wxmaths.h>
using namespace wxmaths;

int main() {

Matrix a;
cout << "Enter a Matrix...";
cin >> a;
cout << a;


Matrix b(transpose(a));
cout << "transpose:\n" << b;

double det;
if (a.square()) {
   det = determinant(a);
   cout << "determinant = " << det << '\n';
}

b = inverse(a);
if (b.rows == 0) cout << "singular!\n";
else {
   cout << "inverse of a =\n" << b;
   Matrix c = a * b;
   cout << a << "*\n" << b << "=\n" << c;
}

Vector y;
cout << "Enter a Vector...";
cin >> y;
cout << "inverse of\n" << a << "*\n" << y << "=\n";
Vector x = b * y;
cout << x;
Vector y1 = a * x;
cout << a << "*\n" << x << "=\n" << y1;

cout << "solution of\n" << a << "* x =\n" << y << " is x=\n";
double res;
x = solve(a,y,res,det);
if (x.n == 0) cout << "singular!\n";
else {
   cout << x << "\nresidual = " << res << " determinant = " << det << '\n';
   y1 = a * x;
   cout << a << "*\n" << x << "=\n" << y1;
}

{
cout << "Enter tridiagonal matrix and rhs vector (rank, Al, Ad, Ah, y)...\n";
int n; cin >> n;
Vector_fixed_dimension Al(n-1), Ad(n), Ah(n-1), y(n);
cin >> Al >> Ad >> Ah >> y;
Vector x = solve_tridiag_thomas(Al,Ad,Ah,y);
cout << "x = "; x.Write_elements(cout); cout << endl;
Vector_fixed_dimension y1(n);
y1[0] = x[0] * Ad[0] + x[1] * Ah[0];
for (int i=1;i<n-1;i++) y1[i] = x[i-1] * Al[i-1] + x[i] * Ad[i] + x[i+1] * Ah[i];
y1[n-1] = x[n-2] * Al[n-2] + x[n-1] * Ad[n-1];
cout << "ycalc = " << y1 << endl;
}

}