matrix_op.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_def.h"

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

//////////////////////////////////////////////////////
/* 1 */ matrix_t matrix_build_zero(int row_size,int col_size);
/* 2 */ matrix_t matrix_smul(matrix_t matrix_a , float coef);
/* 3 */ matrix_t matrix_transpose(matrix_t matrix);
/* 4 */ matrix_t matrix_adj(matrix_t matrix);
/* 5 */ matrix_t matrix_inverse(matrix_t matrix);
/* 6 */ matrix_t matrix_minor(matrix_t matrix_a,int xrow,int xcol);
/* 7 */ matrix_t matrix_add(matrix_t matrix_a , matrix_t matrix_b);
/* 8 */ matrix_t matrix_sub(matrix_t matrix_a , matrix_t matrix_b);
/* 9 */ matrix_t matrix_mul(matrix_t matrix_a , matrix_t matrix_b);
/* 10*/ float flabs(float data);
/* 11*/ int matrix_cofactor(int row, int col);
/* 12*/ float matrix_det(matrix_t matrix_a);
/* 13*/ float matrix_det2(matrix_t matrix_a);
/* 14*/ void matrix_exchange_rows(matrix_t *matrix, int rowx, int rowy);
/* 15*/ void matrix_exchange_cols(matrix_t *matrix, int colx, int coly);
/* 16*/ matrix_t matrix_P_decomposition(matrix_t matrix_a);
/* 17*/ matrix_t matrix_build_Pv(matrix_t matrix_a);
/* 18*/ matrix_t matrix_build_P(matrix_t matrix_a);
/* 19*/ void matrix_build_LU(matrix_t matrix,matrix_t *lower_tm,matrix_t *upper_tm);
/* 20*/ void matrix_build_LU2(matrix_t matrix_a,matrix_t *lower_tm,matrix_t *upper_tm);
/* 21*/ matrix_t matrix_inverse_L(matrix_t matrix_lm);
/* 22*/ matrix_t matrix_inverse_U(matrix_t matrix_um);
/* 23*/ matrix_t matrix_inverse2(matrix_t matrix);
/* 24*/ matrix_t matrix_adj2(matrix_t matrix);

//////////////////////////////////////////////////////

//////////////////////////////////////////////////////
//  M A T R I X   O P E R A T I O N S  
//////////////////////////////////////////////////////
//////////////////////////////////////////////////////
// Definition   : build a zero matrix at the given dims  
//////////////////////////////////////////////////////
matrix_t matrix_build_zero(int row_size,int col_size){
    int row,col;
    matrix_t matrix_zero;
    
    if(!matrix_allocate(&matrix_zero,row_size,col_size)){
        printf("\nin function matrix_build_zero ...\n");
        exit(1);
    }    
    
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
                matrix_set_data(&matrix_zero,row,col,0);
        }
    }
    
    return matrix_zero;
}

//////////////////////////////////////////////////////
// Definition   : scalar mul 
//////////////////////////////////////////////////////
matrix_t matrix_smul(matrix_t matrix_a , float coef){
    int row,col;
    matrix_t matrix_b;
    
    matrix_b = matrix_build_zero(matrix_a.row_size,matrix_a.col_size);
        
    for(row=0;row<matrix_a.row_size;row++){
        for(col=0;col<matrix_a.col_size;col++){
                      
           matrix_set_data(&matrix_b,row,col,coef*matrix_get_data(matrix_a,row,col));
        }
    }
    
    return matrix_b;
}

//////////////////////////////////////////////////////
// Definition   : transpose of given matrix 
//////////////////////////////////////////////////////
matrix_t matrix_transpose(matrix_t matrix){
    int row,col;
    matrix_t matrix_tr;
    
    matrix_tr=matrix_build_zero(matrix.col_size,matrix.row_size);
    
    for(row=0;row<matrix.row_size;row++){
        for(col=0;col<matrix.col_size;col++){
           
           matrix_set_data(&matrix_tr,col,row,matrix_get_data(matrix,row,col));
        }
    }
    
    return matrix_tr;
}

//////////////////////////////////////////////////////
// Definition   : adjoint matrix
//////////////////////////////////////////////////////
matrix_t matrix_adj(matrix_t matrix){
    int row,col,row_size,col_size;
    float data;
    matrix_t adj;
    
    row_size=matrix.row_size;
    col_size=matrix.col_size;
            
    adj=matrix_build_zero(row_size,col_size);
    matrix = matrix_transpose(matrix);
          
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
            data=matrix_cofactor(row,col)*matrix_det2(matrix_minor(matrix,row,col));
            matrix_set_data(&adj,row,col,data);
        }
    }
      
    return adj;
}

//////////////////////////////////////////////////////
// Definition   : matrix inverse 
//////////////////////////////////////////////////////
matrix_t matrix_inverse(matrix_t matrix){
    matrix_t adj,inv;
    int row_size,col_size;
    float det;
    
    row_size=matrix.row_size;
    col_size=matrix.col_size;    
        
    inv = matrix_build_zero(row_size,col_size);
    det = matrix_det2(matrix);
    adj = matrix_adj(matrix);
    
    inv = matrix_smul(adj,1/det);
    return inv;
}

//////////////////////////////////////////////////////
// Definition   : build matrix minor 
//////////////////////////////////////////////////////
matrix_t matrix_minor(matrix_t matrix_a,int xrow,int xcol){
    
    float data;
    int row,col,mrow,mcol;
    matrix_t minor;
    
    mrow=matrix_a.row_size-1;
    mcol=matrix_a.col_size-1;
    
    minor=matrix_build_zero(mrow,mcol);   
    mrow=0;
    mcol=0;
     
    for(row=0;row<matrix_a.row_size;row++){
        for(col=0;col<matrix_a.col_size;col++){
            
            if(row==xrow || col==xcol){
               continue;
            }
            
            data=matrix_get_data(matrix_a,row,col);
            
            if(mcol==minor.col_size){
               mcol=0;
               mrow++;
            }
            matrix_set_data(&minor,mrow,mcol,data);
            mcol++;
                 
        }       
    }
    
    return minor;
}

//////////////////////////////////////////////////////
// Definition   : matrix addition 
//////////////////////////////////////////////////////
matrix_t matrix_add(matrix_t matrix_a , matrix_t matrix_b){
    int row,col;
    matrix_t matrix_c;
        
    matrix_c=matrix_build_zero(matrix_a.row_size,matrix_a.col_size);           
    
    for(row=0;row<matrix_c.row_size;row++){
        for(col=0;col<matrix_c.col_size;col++){
        matrix_set_data(&matrix_c,row,col,
                        matrix_get_data(matrix_a,row,col)+matrix_get_data(matrix_b,row,col)
                        );
        }
    }
    
    return matrix_c;
}

//////////////////////////////////////////////////////
// Definition   : matrix subs 
//////////////////////////////////////////////////////
matrix_t matrix_sub(matrix_t matrix_a , matrix_t matrix_b){
    int row,col;
    matrix_t matrix_c;
    
    matrix_c=matrix_build_zero(matrix_a.row_size,matrix_a.col_size); 
            
    for(row=0;row<matrix_c.row_size;row++){
        for(col=0;col<matrix_c.col_size;col++){
        matrix_set_data(&matrix_c,row,col,
                        matrix_get_data(matrix_a,row,col)-matrix_get_data(matrix_b,row,col)
                        );
        }
    }
    
    return matrix_c;
}

//////////////////////////////////////////////////////
// Definition   : matrix mul 
//////////////////////////////////////////////////////
matrix_t matrix_mul(matrix_t matrix_a , matrix_t matrix_b){
    int row,col,rcol;
    float sum;
    matrix_t matrix_c;
    
    matrix_c=matrix_build_zero(matrix_a.row_size,matrix_b.col_size);
     
    for(row=0;row<matrix_a.row_size;row++){
        for(col=0;col<matrix_b.col_size;col++){
           sum=0;
           for(rcol=0;rcol<matrix_a.col_size;rcol++){
               sum+=matrix_get_data(matrix_a,row,rcol)*matrix_get_data(matrix_b,rcol,col);
           }
           matrix_set_data(&matrix_c,row,col,sum);
        }
    }
    
    return matrix_c;
}

//////////////////////////////////////////////////////
// Definition   : simple abs function 
//////////////////////////////////////////////////////
float flabs(float data){
    if(data<0){
        return -data;
    }
    else{
        return data;
    }
}

//////////////////////////////////////////////////////
// Definition   : cofactor  
//////////////////////////////////////////////////////
int matrix_cofactor(int row, int col){
    if(!((row+col)%2)){
        return 1;
    }
    else{
        return -1;
    }
}

//////////////////////////////////////////////////////
// Definition   : usual matrix det NOT RECOMENDED
//////////////////////////////////////////////////////
float matrix_det(matrix_t matrix_a){
    int index;
    float det;
      
    if(matrix_a.row_size==2){
        return matrix_get_data(matrix_a,0,0)*matrix_get_data(matrix_a,1,1)
              -matrix_get_data(matrix_a,1,0)*matrix_get_data(matrix_a,0,1);
    }
    
    det=0;
    
    for(index=0;index<matrix_a.col_size;index++){
       det+=matrix_cofactor(0,index)*matrix_get_data(matrix_a,0,index)*matrix_det(matrix_minor(matrix_a,0,index));
    }
    
    return det;
}

//////////////////////////////////////////////////////
// Definition   : calculate matrix det using LU factors  
//////////////////////////////////////////////////////
float matrix_det2(matrix_t matrix_a){
    int row,row_size;
    matrix_t upper_tm,lower_tm;
    float det=-1;
        
    matrix_build_LU2(matrix_a,&lower_tm,&upper_tm);
    
    matrix_remove(&lower_tm);
    
    row_size=upper_tm.row_size;
        
    for(row=0;row<row_size;row++){
        det*=matrix_get_data(upper_tm,row,row);
    }
        
    return det;
}