matrix_utils.h

Matrix operations solution of AX=B Jordan and newton Methods

//////////////////////////////////////////////////////
// programmer : tufan biyiktas
// email      : [e-mail]
// date       : March, 7 2005
//////////////////////////////////////////////////////

//////////////////////////////////////////////////////
// HEADER FILES
//////////////////////////////////////////////////////
#include<stdio.h>
#include<stdlib.h>

#include "matrix_op.h"

//////////////////////////////////////////////////////
// FUNCTIONS
//////////////////////////////////////////////////////

/* 1 */ matrix_t least_square(matrix_t matrix_x, matrix_t matrix_y, int degree);
/* 2 */ matrix_t matrix_eq_solve(matrix_t matrix_a,matrix_t matrix_b);
/* 3 */ matrix_t matrix_eq_solve2(matrix_t matrix_a,matrix_t matrix_b);
/* 4 */ int matrix_select_piv(matrix_t matrix,int prow);
/* 5 */ matrix_t matrix_jordan_iterate(matrix_t *matrix, int prow,int pcol);
/* 6 */ matrix_t matrix_jordan_iterate2(matrix_t *matrix, int prow,int pcol);
/* 7 */ matrix_t matrix_jordan(matrix_t matrix_a);
/* 8 */ matrix_t matrix_jordan2(matrix_t matrix_a);
/* 9 */ matrix_t matrix_jordan_solve(matrix_t matrix_a, matrix_t matrix_b);
/* 10*/ matrix_t matrix_jordan_solve2(matrix_t matrix_a, matrix_t matrix_b);
//////////////////////////////////////////////////////

//////////////////////////////////////////////////////
//  L E A S T S  S Q U A R E S
//////////////////////////////////////////////////////
////////////////////////////////////////////////////////
//  Definition   : least square method
////////////////////////////////////////////////////////
matrix_t least_square(matrix_t matrix_x, matrix_t matrix_y, int degree){
    int row,col,row_size,col_size,count;
    float data,mul;
    matrix_t matrix_a,matrix_c;
    
    row_size=matrix_x.row_size;
    col_size=degree+1;
    
    matrix_a=matrix_build_zero(row_size,degree+1);
    matrix_c=matrix_build_zero(row_size,1);
    
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
            data=matrix_get_data(matrix_x,row,0);
            mul=1;
            for(count=0;count<col;count++){
                mul*=data;
            }    
                             
            matrix_set_data(&matrix_a,row,col,mul);
        }
    }            
    
    return matrix_mul(matrix_mul(matrix_inverse2(matrix_mul(matrix_transpose(matrix_a),matrix_a)),matrix_transpose(matrix_a)),matrix_y);
}
//////////////////////////////////////////////////////
//  L E A S T S  S Q U A R E S   ---END---
//////////////////////////////////////////////////////

//////////////////////////////////////////////////////
//  L I N E A R  E Q U A T I O N  S O L V E 
//////////////////////////////////////////////////////     
//////////////////////////////////////////////////////
// Definition   : AX=B, X=inv2(A)B
//////////////////////////////////////////////////////
matrix_t matrix_eq_solve(matrix_t matrix_a,matrix_t matrix_b){
    
    matrix_t inv=matrix_inverse2(matrix_a);
    
    return matrix_mul(inv,matrix_b);
}    
    
////////////////////////////////////////////////////////
//  Definition   : AX=B, PA=PB, LUX=PB
////////////////////////////////////////////////////////     
matrix_t matrix_eq_solve2(matrix_t matrix_a,matrix_t matrix_b){
    int row,xcol,row_size,col_size;
    float data,sum;
    matrix_t lower,upper,perm,matrix_x,matrix_y,matrix_pb;
    
    row_size=matrix_a.row_size;
    col_size=matrix_a.col_size;
    
    matrix_build_LU2(matrix_a,&lower,&upper);
    
    perm=matrix_build_P(matrix_a);
    matrix_x=matrix_build_zero(row_size,1);
    matrix_y=matrix_build_zero(row_size,1);
    
    matrix_pb=matrix_mul(perm,matrix_b);
    
    for(row=0;row<row_size;row++){
        sum=0;
        for(xcol=0;xcol<row;xcol++){
            sum+=matrix_get_data(lower,row,xcol)*matrix_get_data(matrix_y,xcol,0);
        }
        data=(matrix_get_data(matrix_pb,row,0)-sum);
        matrix_set_data(&matrix_y,row,0,data);
    }
    
    for(row=row_size-1;row>=0;row--){
        sum=0;
        for(xcol=row;xcol<row_size;xcol++){
            sum+=matrix_get_data(upper,row,xcol)*matrix_get_data(matrix_x,xcol,0);
        }
        data=(matrix_get_data(matrix_y,row,0)-sum)/matrix_get_data(upper,row,row);
        matrix_set_data(&matrix_x,row,0,data);
    }
    
    return matrix_x;
}
////////////////////////////////////////////////////////
//  L I N E A R  S O L V E   ---END--- 
//////////////////////////////////////////////////////// 

//////////////////////////////////////////////////////
// J O R D A N   E L E M I N A T I O N
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
// Name         : 
// Definition   : pivotun selmesi
//////////////////////////////////////////////////////
int matrix_select_piv(matrix_t matrix,int prow){
    float piv;
    int col,sel;
            
    for(col=prow;col<matrix.col_size;col++){
        piv=matrix_get_data(matrix,prow,col);
    
        if(!piv){
            continue;
        }
        else{
            sel=col;
            return sel;
        }
    }        
    return -1;
}

//////////////////////////////////////////////////////
// Name         : 
// Definition   : alisilmis jordan elemesi icin
//////////////////////////////////////////////////////
matrix_t matrix_jordan_iterate(matrix_t *matrix, int prow,int pcol){
    int row,col,row_size,col_size;
    float piv,data;
        
    row_size=matrix->row_size;
    col_size=matrix->col_size;
        
    piv=matrix_get_data(*matrix,prow,pcol);
                                       
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
                if(row==prow || col==pcol){
                    continue;
                }
                else{
                    data=matrix_get_data(*matrix,row,col)*matrix_get_data(*matrix,prow,pcol);
                    data-=matrix_get_data(*matrix,row,pcol)*matrix_get_data(*matrix,prow,col);
                    matrix_set_data(matrix,row,col,data);
                }    
            }//col
        }//row
                       
        for(col=0;col<col_size;col++){
            if(col==pcol){
                continue;
            }    
            data = matrix_get_data(*matrix,prow,col)*(-1);
            matrix_set_data(matrix,prow,col,data);
        } //col
        
        matrix_set_data(matrix,prow,pcol,1);
        *matrix=matrix_smul(*matrix,1/piv);
}
                               
//////////////////////////////////////////////////////
// Name         : 
// Definition   : alisilmis jordan elemesi
//////////////////////////////////////////////////////
matrix_t matrix_jordan(matrix_t matrix_a){
    int row,col,prow,pcol,row_size,col_size;
    matrix_t matrix, perm, col_trace;
        
    row_size=matrix_a.row_size;
    col_size=matrix_a.col_size;
    
    matrix_copy(matrix_a,&matrix);
    col_trace=matrix_build_zero(row_size,1);
    
    for(row=0;row<row_size;row++){
        matrix_set_data(&col_trace,row,0,row);
    }    
        
    for(row=0;row<row_size;row++){
                
        pcol=matrix_select_piv(matrix,row);
               
        if(pcol==-1){
            error_num=3;
            matrix_error_print();
            exit(1);
        }//if
        
        matrix_set_data(&col_trace,row,0,pcol);
        matrix_jordan_iterate(&matrix,row,pcol);
        
    }
                              	
	for(row=0;row<row_size;row++){
	    prow=matrix_get_data(col_trace,row,0);
	    
	    if(prow==-1){
	       continue;
        }
        
        matrix_set_data(&col_trace,prow,0,-1);
	        
        if(prow!=row){
	       
           matrix_exchange_rows(&matrix,row,prow);
        }
    }
           
    for(col=0;col<col_size;col++){
	    pcol=matrix_get_data(col_trace,col,0);
	    	    
        if(pcol==-1){
	       continue;
        }
        
        matrix_set_data(&col_trace,pcol,0,-1);
        
        if(pcol!=col){
	        
            matrix_exchange_cols(&matrix,col,pcol);
        }
    }
   
   return matrix;
}

//////////////////////////////////////////////////////
// Name         : 
// Definition   : degistirilmis jordan elemesi icin
//////////////////////////////////////////////////////
matrix_t matrix_jordan_iterate2(matrix_t *matrix, int prow,int pcol){
    int row,col,row_size,col_size;
    float piv,data;
        
    row_size=matrix->row_size;
    col_size=matrix->col_size;
        
    piv=matrix_get_data(*matrix,prow,pcol);
                                       
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
                if(row==prow || col==pcol){
                    continue;
                }
                else{
                    data=matrix_get_data(*matrix,row,col)*matrix_get_data(*matrix,prow,pcol);
                    data-=matrix_get_data(*matrix,row,pcol)*matrix_get_data(*matrix,prow,col);
                    matrix_set_data(matrix,row,col,data);
                }    
            }//col
        }//row
                       
        for(row=0;row<row_size;row++){
            if(row==prow){
                continue;
            }    
            data = matrix_get_data(*matrix,row,pcol)*(-1);
            matrix_set_data(matrix,row,pcol,data);
        } //col
        
        matrix_set_data(matrix,prow,pcol,1);
        *matrix=matrix_smul(*matrix,1/piv);
}

//////////////////////////////////////////////////////
// Name         : 
// Definition   : degistirilmis jordan elemesi 
//////////////////////////////////////////////////////
matrix_t matrix_jordan2(matrix_t matrix_a){
    int row,col,prow,pcol,row_size,col_size;
    matrix_t matrix, perm, col_trace;
        
    row_size=matrix_a.row_size;
    col_size=matrix_a.col_size;
    
    matrix_copy(matrix_a,&matrix);
    col_trace=matrix_build_zero(row_size,1);
    
    for(row=0;row<row_size;row++){
        matrix_set_data(&col_trace,row,0,row);
    }    
        
    for(row=0;row<row_size;row++){
                
        pcol=matrix_select_piv(matrix,row);
               
        if(pcol==-1){
            error_num=3;
            matrix_error_print();
            exit(1);
        }//if
        
        matrix_set_data(&col_trace,row,0,pcol);
        matrix_jordan_iterate2(&matrix,row,pcol);
        
    }
                              	
	for(row=0;row<row_size;row++){
	    prow=matrix_get_data(col_trace,row,0);
	    
	    if(prow==-1){
	       continue;
        }
        
        matrix_set_data(&col_trace,prow,0,-1);
	        
        if(prow!=row){
	       
           matrix_exchange_rows(&matrix,row,prow);
        }
    }
           
    for(col=0;col<col_size;col++){
	    pcol=matrix_get_data(col_trace,col,0);
	    	    
        if(pcol==-1){
	       continue;
        }
        
        matrix_set_data(&col_trace,pcol,0,-1);
        
        if(pcol!=col){
	        
            matrix_exchange_cols(&matrix,col,pcol);
        }
    }
   
   return matrix;
}
//////////////////////////////////////////////////////
// J O R D A N   E L E M I N A T I O N   ---END---
////////////////////////////////////////////////////// 

//////////////////////////////////////////////////////
// L I N E A R  S O L V E  W I T H  J O R D A N     
////////////////////////////////////////////////////// 
////////////////////////////////////////////////////////
//  Definition   : AX=B, alisilmis jordan elemesi ile cozum
////////////////////////////////////////////////////////   
matrix_t matrix_jordan_solve(matrix_t matrix_a, matrix_t matrix_b){
    int row,col,row_size,col_size;
    float data;
    matrix_t matrix_x;
    
    row_size = matrix_a.row_size;
    col_size = matrix_a.col_size;
    
    matrix_x = matrix_build_zero(row_size,1);
    
    matrix_a = matrix_jordan(matrix_a);
    matrix_b = matrix_transpose(matrix_b);
        
    for(row=0;row<row_size;row++){
        data=0;
        for(col=0;col<col_size;col++){
            data+=matrix_get_data(matrix_a,row,col)*matrix_get_data(matrix_b,0,col);
        }
        matrix_set_data(&matrix_x,row,0,data);
    }
    
    return matrix_x;
}

////////////////////////////////////////////////////////
//  Definition   : AX=B, degistirilmis jordan elemesi ile cozum
////////////////////////////////////////////////////////   
matrix_t matrix_jordan_solve2(matrix_t matrix_a, matrix_t matrix_b){
    int row,col,row_size,col_size;
    float data;
    matrix_t matrix_x;
    
    row_size = matrix_a.row_size;
    col_size = matrix_a.col_size;
    
    matrix_x = matrix_build_zero(row_size,1);
    
    matrix_a = matrix_jordan2(matrix_a);
    matrix_b = matrix_transpose(matrix_b);
        
    for(row=0;row<row_size;row++){
        data=0;
        for(col=0;col<col_size;col++){
            data+=matrix_get_data(matrix_a,row,col)*matrix_get_data(matrix_b,0,col);
        }
        matrix_set_data(&matrix_x,row,0,data);
    }
    
    return matrix_x;
}
//////////////////////////////////////////////////////
// L I N E A R  S O L V E  W I T H  J O R D A N  ---END---     
//////////////////////////////////////////////////////