#include <stdio.h>
#include <iostream>
#include <math.h>
void create_mat(float*** _mat, int _row, int _col);
void init_mat(float** _mat, int _row, int _col, float _val);
void print_mat(float** _mat, int _row, int _col);
void get_mat_val(float** _mat, int _row, int _col);
void pivoting(float** _matA, float** _matb, float **pivot, int _row, int _col);
void gausselimination(float ** _matA, float ** _matb, float **_mult, int _row, int _col);
void backward_substitution(float **_matA, float **_matb, float **matx, int _row, int _col);//A와 b를 넣으면 행렬 x를 구함.
/*
void LUdecomp(float** X)
void matinv(float** X)
*/
void main() //main 시작---------------
{
//1. declare variables; get matrix size from user;
int M;
int _row, _col;
std::cout << "Ax=b; A(M*M); x(M*1); b(M*1);\n Enter Matrix size M :";
std::cin >> M;
float **matA;
float **matb;
float **matx;
float **mult;
float **pivot;
_row = M;
_col = M;
create_mat(&matA, _row, _col); // function create_mat 실행
create_mat(&matb, _row, 1);
create_mat(&matx, _row, 1);
create_mat(&mult, _row, _col);
create_mat(&pivot, 1, _col+1);
init_mat(matA, _row, _col, 0); // function init_mat 실행
init_mat(matb, _row, 1, 0);
init_mat(matx, _row, 1, 0);
init_mat(mult, _row, _col, 0);
init_mat(pivot, 1, _col+1, 0);
//2. enter matrix values; get matrix values from user;
std::cout << "\nInput matrix A\n";
get_mat_val(matA, _row, _col);
std::cout << "Printing matrix A : \n";
print_mat(matA, _row, _col);
std::cout << "\nInput matrix b\n";
get_mat_val(matb, _row, 1);
std::cout << "Printing matrix b : \n";
print_mat(matb, _row, 1);
print_mat(matb, _row, 1);
std::cout << "\n";
//3. pivoting
pivoting(matA, matb, pivot, _row, _col);
std::cout << "\nPivoted A,b\n";
print_mat(matA, _row, _col);
print_mat(matb, _row, 1);
print_mat(pivot, 1, _col + 1);
//GaussElimination
gausselimination(matA, matb, mult, _row, _col);
//std::cout << "\nGauss ellimination A\n";
//print_mat(matA, _row, _col);
//std::cout << "\nmult array\n";
//print_mat(mult, _row, _col);
//std::cout << "permutaion b\n";
//print_mat(matb, _row, 1);
//backwardsubstitution
backward_substitution(matA, matb, matx, _row, _col);
//std::cout << "x\n";
//print_mat(matx, _row, 1);
//LU decomposition 쉽다 나중에 만들자
//void LUdecomp(float** X)
//Inverse matrix
//void matinv(float** X)
system("pause");
free(matA); // Free allocated memory when program ends
free(matb);
free(matx);
free(mult);
free(pivot);
} //main 끝-----------------------------
// 사용자에게 받은 row와 col값 이용, 삼중포인터 mat에 배열을 만들어서 주소를 저장시킴
void create_mat(float*** _mat, int _row, int _col)
{
int i;
/* 1. allocate row first */
*_mat = (float**)malloc(sizeof(float) * (_row));
/* 2. allocate column */
for (i = 0; i < _row; i++)
(*_mat)[i] = (float*)malloc(sizeof(float) * (_col));
}// matrix 생성 함수; m*m 크기로 생성
// 만들어진 배열에 저장되어있는 쓰레기값들을 버리고 지정한 어떤 수 _val로 초기화
void init_mat(float** _mat, int _row, int _col, float _val)
{
int i, j;
for (i = 0; i < _row; i++)
{
for (j = 0; j < _col; j++)
(_mat)[i][j] = _val;
}
}
// 배열에 저장되어있는 값들을 출력; 행렬 출력
void print_mat(float** _mat, int _row, int _col)
{
int i, j;
for (i = 0; i < _row; i++)
{
for (j = 0; j < _col; j++)
printf("%f, ", _mat[i][j]);
printf("\n");
}
printf("==== ==== ==== ==== ==== ==== ==== ====\n");
}// matrix 프린트 함수;
// 사용자에게서 입력을 받아서 만들어진 배열에 값을 저장
void get_mat_val(float** _mat, int _row, int _col)
{
printf("Enter the matrix values\n");
int i,j;
for (i=0; i<_row; i++)
{
printf("row[%d]\n",i);
for(j = 0; j <_col; j++)
{
printf("[%d][%d]=",i,j);
std::cin >> _mat[i][j];
}
}
printf("\n");
}
void pivoting(float **_matA, float **_matb, float **pivot, int _row, int _col)
{
int i, j, K,R;
for (K = 0; K < _row; K++)
{
if (_matA[K][K] == 0)//<pivoting
{
for (j = 0; j < _col; j++)
{
pivot[0][j] = _matA[K][j];
}
pivot[0][_col] = _matb[K][0];
for (i = 1; i < _col - K && _matA[K + i][K] != 0; i++)
{
R = i + K;
printf("%d", R);
}
for (j = 0; j < _col; j++)
{
_matA[K][j] = _matA[R][j];
_matb[K][0] = _matb[R][0];
_matA[R][j] = pivot[0][j];
_matb[R][0] = pivot[0][_col];
}
}//pivoting>
}
}
// Ax=b에서 'A와 'b'에 해당하는 행렬의 주소값을 받아서, A와 b에가우스 소거법을 적용
void gausselimination(float **_matA, float **_matb, float **_mult, int _row, int _col)
{
int K;
int i, j;
for (K = 0; K < _row; K++)
{
for (i = K + 1; i < _row; i++)
{
_mult[i][K] = (_matA[i][K] / _matA[K][K]); // mult 값은 열과 상관없으므로 연산을 줄이기위해 앞으로 뺌
_matb[i][0] = _matb[i][0] - _mult[i][K] * _matb[K][0]; //A행렬의 연산을 b행렬에도 적용
for (j = K + 1; j < _col; j++)
{
_matA[i][j] = _matA[i][j] - _mult[i][K] * _matA[K][j];//row[i]의 2번째 3...번째 값들
}
_matA[i][K] = 0; //A[K][K]와 같은 열에 있는 값들 중 더 낮은 row의 값들은 연산 후 0이 될 것임.
}
}
}
// Ax=b에서 행렬 A가 upper triangular matrix일 때, backward substitution을 사용하여 x를 계산
void backward_substitution(float **_matA, float** _matb, float** _matx, int _row, int _col)
{
int i, j;
for (i = 0; i < _row; i++)
_matx[i][0] = _matb[i][0];// 현재 행렬 b의 데이터를 건드리지 않기 위해 행렬 x로 옮김
for (i = _row-1; i >=0; i--) //backward substituion 진행
{
_matx[i][0] = _matx[i][0] / _matA[i][i];
for (j = 1;j< _row-i; j++)
{
_matx[i][0] = _matx[i][0] - (_matA[i][i + j] * _matx[i + j][0]) / _matA[i][i];
}
}
}
/*Ax=b A=LU; should be acting after GaussElimination L=mult+E(_row) U=Gausseliminated 'A'
void LUdecomp(float** X)
{
mult+E -> L
Gausseliminated A -> U
}
*/
/*
void matinv(float** X)
{
}
*/
글쓰기