hardila

/**gaussElim Exercise***/
#include<stdio.h>
#include<math.h>
/*******
 Function that performs Gauss-Elimination and returns the Upper triangular matrix and solution of equations:
There are two options to do this in C.
1. Pass the augmented matrix (a) as the parameter, and calculate and store the upperTriangular(Gauss-Eliminated Matrix) in it.
2. Use malloc and make the function of pointer type and return the pointer.
This program uses the first option.
********/
void gaussEliminationLS(int m, int n, double a[m][n], double x[n-1][1]){
    int i,j,k;
    for(i=0;i<m-1;i++){
        //Partial Pivoting
        for(k=i+1;k<m;k++){
            //If diagonal element(absolute vallue) is smaller than any of the terms below it
            if(fabs(a[i][i])<fabs(a[k][i])){
                //Swap the rows
                for(j=0;j<n;j++){                
                    double temp;
                    temp=a[i][j];
                    a[i][j]=a[k][j];
                    a[k][j]=temp;
                }
            }
        }
        //Begin Gauss Elimination
        for(k=i+1;k<m;k++){
            double  term=a[k][i]/ a[i][i];
            for(j=0;j<n;j++){
                a[k][j]=a[k][j]-term*a[i][j];
            }
        }
         
    }
    //Begin Back-substitution
    for(i=m-1;i>=0;i--){
        x[i][0]=a[i][n-1];
        for(j=i+1;j<n-1;j++){
            x[i][0]=x[i][0]-a[i][j]*x[j][0];
        }
        x[i][0]=x[i][0]/a[i][i];
    }
             
}
/*******
Function that generates the Hilbert matrix
Parameters: order (n) ,matrix[n][n] 
*******/
void genMatrix(int n, double matrix[n][n]){
    int i,j;
    //Initialize Coefficients
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            matrix[i][j]=(double)1.0/((i+1)+(j+1)-1);
        }
    }
}
/*******
Function that generates the Augmented Hilbert matrix
Parameters: order (n) ,matrix[n][n+1] 
*******/
void genAugMatrix(int n, double matrix[n][n+1]){
    int i,j;
    //Initialize Coefficients
    for(i=0;i<n;i++){
        for(j=0;j<n;j++){
            matrix[i][j]=(double)1.0/((i+1)+(j+1)-1);
        }
    }
    //Initialize RHS part 
    for(i=0;i<n;i++){
        matrix[i][n]=1;
    }
}
/*******
Function that prints the elements of a matrix row-wise
Parameters: rows(m),columns(n),matrix[m][n] 
*******/
void printMatrix(int m, int n, double matrix[m][n]){
    int i,j;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++){
            printf("%lf\t",matrix[i][j]);
        }
        printf("\n");
    } 
}
/*******
Function that copies the elements of a matrix to another matrix
Parameters: rows(m),columns(n),matrix1[m][n] , matrix2[m][n]
*******/
void copyMatrix(int m, int n, double matrix1[m][n], double matrix2[m][n]){
    int i,j;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++){
            matrix2[i][j]=matrix1[i][j];
        }
    } 
}
/*******
 Function that calculates the product of two matrices:
There are two options to do this in C.
1. Pass a matrix (prod) as the parameter, and calculate and store the product in it.
2. Use malloc and make the function of pointer type and return the pointer.
This program uses the first option.
********/
void matProduct(int m, int n, int n1,double a[m][n1], double b[n1][n], double prod[m][n]){
    int i,j,k;
    for(i=0;i<m;i++){
        for(j=0;j<n;j++){
            prod[i][j]=0;
            for(k=0;k<n1;k++){
                prod[i][j]=prod[i][j]+a[i][k]*b[k][j];  
            }
        }
    }   
}
 
int main(){
    int n,i,j;
    printf("Enter the order:\n(n)\n");
    scanf("%d",&n);
    //Declare a matrix to store the augmented Hilbert matrix for the problem
    double a[n][n+1];
    //Declare another matrix to store the resultant matrix obtained after Gauss Elimination
    double U[n][n+1];
    //Declare an array to store the solution of equations
    double x[n][1];
    genAugMatrix(n,a);
    printf("The auto-generated augmented Hilbert matrix for the problem is:\n\n");
    printMatrix(n,n+1,a);
    copyMatrix(n,n+1,a,U);
    //Perform Gauss Elimination 
    gaussEliminationLS(n,n+1,U,x);
    printf("\nThe Upper Triangular matrix after Gauss Eliminiation is:\n\n");
    printMatrix(n,n+1,U);
    printf("\nThe solution of linear equations is:\n\n");
    for(i=0;i<n;i++){
        printf("x[%d]=\t%lf\n",i+1,x[i][0]);
    }
    //Now we will verify if the answer is correct by multiplying X (solution) with the Hilbert matrix and see if we get 1.
    double B[n][1]; //matrix to store product: A.X=B 
    //Declare a matrix to store the un-augmented(square) Hilbert matrix for the problem
    double a1[n][n];
    genMatrix(n,a1);
    matProduct(n,1,n,a1,x,B);
    //Print the product to verify
    printf("\nThe product of matrix A.X=B:\n\n");
    printMatrix(n,1,B);
}