#include "matwekop.h"
#include "memfunk.h"
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>


void gcw(long double **A, int m, int n){

   long double l;
   int wmax;

	for(int i = 0; i < m; i++)
   	A[i][m + n] = i;

   for(int i = 0; i < m; i++){
		wmax = absmax(A, m, m, i, i);
      swapmat(A, m, m + n + 1, i, wmax);
      for(int j = i + 1; j < m; j++){
         if(A[i][i])
		      l = A[j][i]/A[i][i];
	      else{
   	   	cout << "Blad w module matwkop.h w funkcji gcw,\n macierz jest osobliwa, dzielenie przez zero\n";
      	   getchar();
         	exit(-1);
      	}
         odejmij(A, m, m + n, j, i, l);
         A[j][i] = l;
      }
   }
}
//------------------------------------------------------------------------------
void odwroc(long double **A, long double **A_1, int m){

   long double **L;
   long double **U;
   long double **LU;
	long double *x;
   int nr_kolumny;

   L = alokmat(m, m);
   U = alokmat(m, m);
   LU = alokmat(m, m + m + 1);
   x = alokwek(m);

   zerujmat(L, m, m);
   zerujmat(U, m, m);
   zerujmat(LU, m, m + m + 1);
   zerujwek(x, m);

	gcw(A, m, 0);

   for(int i = 0; i < m; i++)
   	for(int j = 0; j < m; j++)
      	if(i == j){
         	L[i][j] = 1;
            U[i][j] = A[i][j];
         }
         else if (i > j){
         	L[i][j] = A[i][j];
            U[i][j] = 0;
         }
         else{
         	U[i][j] = A[i][j];
         	L[i][j] = 0;
         }

   mnozmat(L, m, m, U, m, m, LU);

   for(int i = 0; i < m; i++)
   	for(int j = m; j < m + m; j++)
      	if( (i + m) == j)
         	LU[i][j] = 1;

   gcw(LU, m, m);

   for(int i = 1; i <= m; i++){
	   rozwiaz(LU, x, m, i);
      nr_kolumny = A[i - 1][m];
      for(int j = 0; j < m; j++)
      	A_1[j][nr_kolumny] = x[j];
   }

	dealokmat(L, m);
   dealokmat(U, m);
   dealokmat(LU, m);
   dealokwek(x);

}
//------------------------------------------------------------------------------
void rozwiaz(long double **A, long double *x, int m, int n){

	long double suma;

	for(int i = (m - 1); i >= 0; i--){
   	suma = 0.0;
		for(int k = i + 1; k < m; k++)
      	suma += A[i][k] * x[k];
      if(A[i][i])
	      x[i] = (A[i][m - 1 + n] - suma) / A[i][i];
      else{
      	cout << "Blad w module matwkop.h w funkcji rozwiaz,\n macierz jest osobliwa, dzielenie przez zero\n";
         getchar();
         exit(-1);
      }
   }
}
//------------------------------------------------------------------------------
void mnozmat(long double **A,int am, int an, long double **B, int bm, int bn,
				 long double **C){

   long double suma;

	if(bm == an){
   	for(int i = 0; i < am; i++)
      	for(int j = 0; j < bn; j++){
         	suma = 0.0;
         	for(int k = 0; k < an; k++)
            	suma += A[i][k] * B[k][j];
         	C[i][j] = suma;
         }
   }
   else{
   	cout << "Blad w module matwekop.h, w funkcji mnozmat,\nnie zgadzaja sie wiersze i kolumny macierzy\n";
      getchar();
      exit(-1);
   }
}
//------------------------------------------------------------------------------
void matwek(long double **A, int am, int an, long double *W, int wm,
				long double *W2){

   long double suma;

	if(an == wm){
   	for(int i = 0; i < am; i++){
        	suma = 0.0;
        	for(int k = 0; k < an; k++)
           	suma += A[i][k] * W[k];
        	W2[i] = suma;
      }
   }
   else{
   	cout << "Blad w module matwekop.h, w funkcji matwek,\nnie zgadzaja sie wiersze i kolumny wektora\n";
      getchar();
      exit(-1);
   }
}
//------------------------------------------------------------------------------
void transp(long double **A, long double **T, int am, int an){

	for(int i = 0; i < an; i++)
   	for(int j = 0; j < am; j++)
      	*(*(T + i) +j) = *(*(A + j) +i);
}
//------------------------------------------------------------------------------
void swapmat(long double **A, int am, int an, int wk, int wj){

	long double temp;
   if(wk < am && wj < am){
	   for(int i = 0; i < an; i++){
   	   temp = A[wk][i];
      	A[wk][i] = A[wj][i];
      	A[wj][i] = temp;
   	}
   }
   else{
   	cout << "Blad w module matwekop.h w funkcji swapmat,\n chcesz przestawic wiersze ktore nie istnieja\n";
      getchar();
      exit(-1);
   }
}
//------------------------------------------------------------------------------
void swapwek(long double *W1, int am ,int wk, int wj){

	long double temp;
   if(wk < am && wj < am){
   	temp = W1[wk];
      W1[wk] = W1[wj];
      W1[wj] = temp;
   }
   else{
   	cout << "Blad w module matwekop.h w funkcji swapwek,\n chcesz przestawic wiersze ktore nie istnieja\n";
      getchar();
      exit(-1);
   }
}
//------------------------------------------------------------------------------
int absmax(long double **A, int m, int n, int wm, int kn){

   long double max = 0.0;
   int wm_max;

   if(wm < m && kn < n){
   	max = A[wm][kn];
        wm_max = wm;
		for(int i = wm; i < m; i++)
   		if(A[i][kn] >= max){
      		max = A[i][kn];
         	wm_max = i;
      	}
      	return wm_max;
   }
   else{
   	cout << "Blad w module matwekop.h w funkcji absmax,\n chesz szukac poza zakresem\nkolumn lub wierszy\n";
      getchar();
      exit(-1);
   }
}
//------------------------------------------------------------------------------
long double ilskal(long double *W1, int w1m, long double *W2){

	int k;
   int m;
   long double suma = 0.0;

   k = w1m / 4;
   m = w1m % 4;

   while(k--){
   	suma += *W1++ * *W2++;
   	suma += *W1++ * *W2++;
   	suma += *W1++ * *W2++;
   	suma += *W1++ * *W2++;
   }
   while(m--)
   	suma += *W1++ * *W2++;

   return suma;
}
//------------------------------------------------------------------------------
void odejmij(long double **A, int am, int an, int wj, int wk, long double lb){

   if(wj < am){
		for(int i = wk; i < an; i++)
      	A[wj][i] -= A[wk][i] * lb;
   }
   else{
   	cout << "Blad w module matwekop.h w funkcji wmatlb,\n chcesz przemozyc wiersz ktory nie istnieja\n";
      getchar();
      exit(-1);
   }
}
//------------------------------------------------------------------------------
void zerujmat(long double **A, int am, int an){

	for(int i = 0; i < am; i++)
   	for(int j = 0; j < an; j++)
      	A[i][j] = 0.00;
}
//------------------------------------------------------------------------------
void zerujwek(long double *W, int wm){

   	for(int j = 0; j < wm; j++)
      	W[j] = 0.000;
}
//------------------------------------------------------------------------------
void drukmat(long double **A, int am, int an){

	for(int i = 0; i < am; i++){
   	for(int j = 0; j < an; j++)
      	cout << *(*(A + i) +j) << " ";
      cout << "\n";
   }
}
//------------------------------------------------------------------------------
void drukwek(long double *W, int wm){

	for(int i = 0; i < wm; i++)
     	cout << *(W + i) << " ";
}
//------------------------------------------------------------------------------

