matrix_op.h

Matrix operations solution of AX=B Jordan and newton Methods

//////////////////////////////////////////////////////
// Definition   : exchange given two rows 
//////////////////////////////////////////////////////
void matrix_exchange_rows(matrix_t *matrix, int rowx, int rowy){
    double data;
    int col;
    
    for(col=0;col<matrix->col_size;col++){
        data = matrix_get_data(*matrix,rowx,col);
        matrix_set_data(matrix,rowx,col,matrix_get_data(*matrix,rowy,col));
        matrix_set_data(matrix,rowy,col,data);
        
    }
}

//////////////////////////////////////////////////////
// Definition   : exchange given two rows 
//////////////////////////////////////////////////////
void matrix_exchange_cols(matrix_t *matrix, int colx, int coly){
    double data;
    int row;
    
    for(row=0;row<matrix->row_size;row++){
        data = matrix_get_data(*matrix,row,colx);
        matrix_set_data(matrix,row,colx,matrix_get_data(*matrix,row,coly));
        matrix_set_data(matrix,row,coly,data);
        
    }
}

//////////////////////////////////////////////////////
// Definition   : build permutation vector DONT USE DIRECTLY 
//////////////////////////////////////////////////////
matrix_t matrix_P_decomposition(matrix_t matrix_a){

    int row,col,xcol,mcol,row_size,col_size;
    float piv,data,denum,temp;
    matrix_t perm_v;
                            
    row_size = matrix_a.row_size;
    col_size = matrix_a.col_size;
        
    perm_v=matrix_build_zero(row_size,1);
                    
    for(row=0;row<row_size;row++){
        matrix_set_data(&perm_v,row,0,row);
    }
            
    for(col=0;col<col_size;col++){
        
        piv=0.0;
        for(row=col;row<row_size;row++){
            
            data=flabs(matrix_get_data(matrix_a,row,col));
            if(data>piv){
               piv=data;
               mcol=row;
            } //end if
        } //end row
            
        if(!piv){
           error_num=3;
           matrix_error_print();
           exit(1);
        }
        
        /**************************************/
        data=matrix_get_data(perm_v,col,0);
        temp=matrix_get_data(perm_v,mcol,0);
        
        matrix_set_data(&perm_v,col,0,temp);
        matrix_set_data(&perm_v,mcol,0,data);
        /**************************************/

        matrix_exchange_rows(&matrix_a,col,mcol);
        
        for(row=col+1;row<row_size;row++){
            
            denum=matrix_get_data(matrix_a,row,col)/matrix_get_data(matrix_a,col,col);
            matrix_set_data(&matrix_a,row,col,denum);
                  
            for(xcol=col+1;xcol<col_size;xcol++){
                                                
                data  = matrix_get_data(matrix_a,row,col);
                denum = matrix_get_data(matrix_a,col,xcol);
                data  = data * denum;
                
                denum = matrix_get_data(matrix_a,row,xcol);
                data = denum-data;
                matrix_set_data(&matrix_a,row,xcol,data);
                
            } //end xcol
        } //end row
    } //end col
    return perm_v;
}

//////////////////////////////////////////////////////
// Definition   : build permutation vector USE IT
//////////////////////////////////////////////////////
matrix_t matrix_build_Pv(matrix_t matrix_a){
    int row_size,col_size;
    matrix_t matrix;
    
    
    row_size=matrix_a.row_size;
    col_size=matrix_a.col_size;
    
    matrix_copy(matrix_a,&matrix);               
    return matrix_P_decomposition(matrix);
}

//////////////////////////////////////////////////////
// Definition   : build permutation matrix 
//////////////////////////////////////////////////////
matrix_t matrix_build_P(matrix_t matrix_a){
    int row,col,row_size,col_size;
    matrix_t perm,perm_v;
    
    
    row_size=matrix_a.row_size;
    col_size=matrix_a.col_size;
    
    perm_v=matrix_build_Pv(matrix_a);        
    perm=matrix_build_zero(row_size,col_size);
    
     for(row=0;row<row_size;row++){
          col=matrix_get_data(perm_v,row,0);
          matrix_set_data(&perm,row,col,1);
     }
     
     return perm;
}

//////////////////////////////////////////////////////
// Definition   : build LU matricies A=LU
//////////////////////////////////////////////////////
void matrix_build_LU(matrix_t matrix,matrix_t *lower_tm,matrix_t *upper_tm){
    int row,col,ind,row_size,col_size;
    float temp;
    matrix_t matrix_a;
    
    matrix_copy(matrix,&matrix_a);
   
    row_size = matrix_a.row_size;
    col_size = matrix_a.col_size;
    
	*upper_tm=matrix_build_zero(row_size,col_size);
	*lower_tm=matrix_build_zero(row_size,col_size);
	
    for(col=0;col<col_size;col++){
          matrix_set_data(upper_tm,col,col,matrix_get_data(matrix_a,col,col));
          
          for(row=col;row<row_size;row++){
              matrix_set_data(lower_tm,row,col,matrix_get_data(matrix_a,row,col)/matrix_get_data(*upper_tm,col,col));
              matrix_set_data(upper_tm,col,row,matrix_get_data(matrix_a,col,row));
          }
          for(row=col;row<row_size;row++){
              for(ind=col;ind<row_size;ind++){
                  temp = matrix_get_data(matrix_a,row,ind)- matrix_get_data(*lower_tm,row,col) * matrix_get_data(*upper_tm,col,ind);
                     
				  matrix_set_data(&matrix_a,row,ind,temp);
              }//ind
          }//row
     }//col 
}       

//////////////////////////////////////////////////////
// Definition   : LUP decomposition PA=LU
//////////////////////////////////////////////////////
void matrix_build_LU2(matrix_t matrix_a,matrix_t *lower_tm,matrix_t *upper_tm){
    int row,col,row_size,col_size;
    float data;
    matrix_t matrix;
       
    row_size = matrix_a.row_size;
    col_size = matrix_a.col_size;
       
    *upper_tm=matrix_build_zero(row_size,col_size);
    *lower_tm=matrix_build_zero(row_size,col_size); 
    
    matrix_copy(matrix_a,&matrix);
    
    matrix_P_decomposition(matrix);
    
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
                data=matrix_get_data(matrix,row,col);
                if(row<col){
                   matrix_set_data(upper_tm,row,col,data);
                   matrix_set_data(lower_tm,row,col,0);
                }
                else if(row==col){
                        matrix_set_data(upper_tm,row,col,data);
                        matrix_set_data(lower_tm,row,col,1);
                     }
                     else{
                        matrix_set_data(upper_tm,row,col,0);
                        matrix_set_data(lower_tm,row,col,data);
                     }
        }
    }        
}
    
//////////////////////////////////////////////////////
// Definition   : build a zero matrix at the given dims  
//////////////////////////////////////////////////////
matrix_t matrix_build_id(int row_size,int col_size){
    int row,col;
    matrix_t matrix_id;
    
    matrix_id=matrix_build_zero(row_size,col_size);
    
    for(row=0;row<row_size;row++){
        for(col=0;col<col_size;col++){
            if(row==col){
				matrix_set_data(&matrix_id,row,col,1);
			}
			else{
				matrix_set_data(&matrix_id,row,col,0);
			}
        }
    }
    
    return matrix_id;
}
//////////////////////////////////////////////////////
// Definition   : calculate inverse of the given L matrix  
//////////////////////////////////////////////////////
matrix_t matrix_inverse_L(matrix_t matrix_lm){
    matrix_t id,inv;
    int row,col,xcol,row_size,col_size;
    float sum,data;
    
    row_size=matrix_lm.row_size;
    col_size=matrix_lm.col_size;
        
    inv=matrix_build_zero(row_size,col_size);
    id=matrix_build_id(row_size,col_size);
            
    for(col=0;col<col_size;col++){
        for(row=0;row<row_size;row++){
            sum=0;
            for(xcol=col;xcol<row;xcol++){
                sum+=matrix_get_data(matrix_lm,row,xcol)*matrix_get_data(inv,xcol,col);
            }
            data=(matrix_get_data(id,row,col)-sum)/matrix_get_data(matrix_lm,row,row);
            matrix_set_data(&inv,row,col,data);
        }
    }
        
    return inv;
}

//////////////////////////////////////////////////////
// Definition   : calculate inverse of the given U matrix  
//////////////////////////////////////////////////////
matrix_t matrix_inverse_U(matrix_t matrix_um){
    matrix_t id,inv;
    int row,col,xcol,row_size,col_size;
    float sum,data;
    
    row_size=matrix_um.row_size;
    col_size=matrix_um.col_size;
    
    
    inv=matrix_build_zero(row_size,col_size);
    id=matrix_build_id(row_size,col_size);
            
    for(col=col_size-1;col>=0;col--){
        for(row=col;row>=0;row--){
            sum=0;
            for(xcol=row;xcol<col+1;xcol++){
                sum+=matrix_get_data(matrix_um,row,xcol)*matrix_get_data(inv,xcol,col);
            }
            data=(matrix_get_data(id,row,col)-sum)/matrix_get_data(matrix_um,row,row);
            matrix_set_data(&inv,row,col,data);
        }
    }
        
    return inv;
}

//////////////////////////////////////////////////////
// Definition   : matrix inverse using choleski method  
//////////////////////////////////////////////////////
matrix_t matrix_inverse2(matrix_t matrix){
    matrix_t lower_tm,upper_tm;
    
    matrix_build_LU(matrix,&lower_tm,&upper_tm);
    
    return matrix_mul(matrix_inverse_U(upper_tm),matrix_inverse_L(lower_tm));
}

//////////////////////////////////////////////////////
// Definition   : adjoint matrix using inverse matrix
//////////////////////////////////////////////////////
matrix_t matrix_adj2(matrix_t matrix){
        
    float det=-1;
    int row,row_size;
    matrix_t lower_tm,upper_tm;
    
    row_size=matrix.row_size;
    matrix_build_LU(matrix,&lower_tm,&upper_tm);
    
    for(row=0;row<row_size;row++){
        det*=matrix_get_data(upper_tm,row,row);
    }
    
    return matrix_smul(matrix_mul(matrix_inverse_U(upper_tm),matrix_inverse_L(lower_tm)),1/det);
    
}
//////////////////////////////////////////////////////
//  M A T R I X   O P E R A T I O N S   ---END---  
//////////////////////////////////////////////////////