矩阵的一些运算

#ifndef MATRIXC
#define MATRIXC

#include <stdio.h>
#include "matrix.h"

#define ISEMPTY(matrix) ((matrix)->line==0||(matrix)->row==0)

void print(double matrix[MAX_MATRIX][MAX_MATRIX],int line,int row)
{
    int i,j;
    for(i=0;i<line;i++)
    {
        for(j=0;j<row;j++)
        {
            printf("%.2lf/t",matrix[i][j]);
        }
        printf("%s","/n/n/n");
    }
}

void read(FILE* fp,double matrix[MAX_MATRIX][MAX_MATRIX],int line,int row)
{
    int i,j;
    for(i=0;i<line;i++)
    {
        for(j=0;j<row;j++)
        {

            fscanf(fp,"%lf",&matrix[i][j]);
        }
    }
}

Matrix* initMatrix(double a[MAX_MATRIX][MAX_MATRIX], size_t m, size_t n)
{
 Matrix *matrix=(Matrix*)malloc(sizeof(Matrix));
 size_t i,j;
 for(i=0;i<m;i++)
 {
  for (j=0;j<n;j++)
  {
   (matrix->el)[i][j]=a[i][j];
  }
 }
 matrix->line=m;
 matrix->row=n;
 return matrix;
}

Matrix* Identity(size_t m)
{
    size_t i,j;
    double mat[MAX_MATRIX][MAX_MATRIX];
    for(i=0;i<m;i++)
    {
        for(j=0;j<m;j++)
        {
            mat[i][j]=0;
        }
    }
    for(i=0;i<m;i++)
    {
        mat[i][i]=1;
    }
    return initMatrix(mat,m,m);
}

size_t isEmptyMatrix(Matrix *a)
{
    return a->line==0||a->row==0;
}

int writeMatrix(Matrix *a, FILE *fp)
{
    int i;
    for(i=0;i<a->line;i++)
    {
        fwrite(a->el+i,sizeof(double),a->row,fp);
    }
    return a->line*a->row;
}

int readMatrix(Matrix *a, size_t m, size_t n, FILE *fp)
{
    int i;
    for(i=0;i<m;i++)
    {
        fread(a->el+i,sizeof(double),n,fp);
    }
    a->line=m;
    a->row=n;
    return m*n;
}

Matrix *addMatrix(Matrix *a, Matrix *b)
{
    int i,j;
    if((a->line!=b->line)||(a->row!=b->row)) return 0;
    for(i=0;i<a->line;i++)
    {
        for(j=0;j<a->row;j++)
        {
            a->el[i][j]+=b->el[i][j];
        }
    }
    return a;
}

int cpymatrix(double matrix1[MAX_MATRIX][MAX_MATRIX],double temp[MAX_MATRIX][MAX_MATRIX],int line,int row)
{
   int i,j;
   for(i=0;i<line;i++)
   {
      for(j=0;j<row;j++)
      {
         matrix1[i][j]=temp[i][j];
      }
   }
   return 1;
}

Matrix *multiMatrix(Matrix *a, Matrix *b)
{
    int i,j,k,line=a->line,mid=a->row,row=b->row;
    double result;
    double temp[MAX_MATRIX][MAX_MATRIX];
    if(a->row!=b->line) return 0;
    for(i=0;i<line;i++)
    {
        for(j=0;j<row;j++)
        {
            result=0;
            for(k=0;k<mid;k++)
            {
                result= result+a->el[i][k]*b->el[k][j];
            }
            temp[i][j]=result;
        }
    }
    cpymatrix(a->el,temp,line,row);
    a->row=b->row;
    return a;
}

int swap(double *a,double *b)
{
    double temp;
    if(a==NULL||b==NULL)
    {
        printf("Error in swap!");
        return -1;
    }
    temp=*a;
    *a=*b;
    *b=temp;
    return 1;
}

Matrix *transMatrix(Matrix *a)
{
    int line=a->line,row=a->row;
    int k=(line<row?line:row);
    int i,j;

    a->line=row;
    a->row=line;

    for(i=0;i<k-1;i++)
    {
        for(j=i+1;j<k;j++)
        {
            swap(&a->el[i][j],&a->el[j][i]);
        }
    }

    if(line>k)
    {
        for(i=k;i<line;i++)
        {
            for(j=0;j<row;j++)
            {
                a->el[j][i]=a->el[i][j];
            }
        }
    }

    if(row>k)
    {
        for(i=0;i<line;i++)
        {
            for(j=k;j<row;j++)
            {
                a->el[j][i]=a->el[i][j];
            }
        }
    }

    return a;
}

int contrary(double matrix[MAX_MATRIX][MAX_MATRIX],int line)
{
 int i,j,k,m;
    double a[MAX_MATRIX][MAX_MATRIX],temp;

 /*生成一个单位矩阵a[line][line]*/
 for(i=0;i<line;i++)
 {
  for(j=0;j<line;j++)
  {
   a[i][j]=0.0;
  }
 }
 for(i=0;i<line;i++) a[i][i]=1.0;

    for(i=0;i<line;i++)
    {
        /*找到非元*/
        for(j=i;j<line;j++)
            if(matrix[j][i]!=0) break;
        if(j!=i)
        {
            /*把matrix[j][i]移到matrix[i][i]处*/
            for(k=i;k<line;k++)
            {
                swap(&matrix[j][k],&matrix[i][k]);
            }
   for(k=0;k<line;k++)
   {
    swap(&a[j][k],&a[i][k]);
   }
        }
        for(j=i+1;j<line;j++)
        {
            temp=matrix[j][i]/matrix[i][i];
   /*把matrix[i][i]下的元素全部化为0*/
            for(k=i;k<line;k++)
            {
                matrix[j][k]=matrix[j][k]-temp*matrix[i][k];
            }
    /*对a[LINE][ROW]作同样的变换*/
   for(k=0;k<line;k++)
            {
                a[j][k]=a[j][k]-temp*a[i][k];
            }
        }
    }

    for(m=line-1;m>=0;m--)
    {
        /*把matrix[LINE][ROW]化成对角形*/
        for(j=m-1;j>=0;j--)
        {
            temp=matrix[j][m]/matrix[m][m];
            matrix[j][m]=0;
            for(k=0;k<line;k++)
            {
                a[j][k] = a[j][k]- temp*a[m][k];
            }

        }
    }

    /*把matrix[LINE][ROW]化成一个单位矩阵*/
    for(i=0;i<line;i++)
    {
        for(j=0;j<line;j++)
        {
            a[i][j]=a[i][j]/matrix[i][i];
        }
    }
 cpymatrix(matrix,a,line,line);
    return 1 ;
}

Matrix *inverseMatrix(Matrix *a)
{
    contrary(a->el,a->line);
    return a;
}

int order(double matrix[MAX_MATRIX][MAX_MATRIX],int line,int row)
{

    int i,j,k,m;
    double temp;
    for(i=0;i<line;i++)
    {
        /*找到非元素*/
        for(m=i;m<line;m++)
        {
            k=0;
            for(j=i;j<row;j++)
            {
                if(matrix[m][j]!=0)
                {
                    k=1;
                    break;
                }
            }
            if(k==1) break;
        }
        if(m==line) break;
        /*把非0元素matrix[m][j]移到matrix[i][i]处*/
        for(k=i; k < row; k++)
  {
   swap(&matrix[m][k],&matrix[i][k]);
  }
        for(k=0;k<line;k++)
  {
   swap(&matrix[k][j],&matrix[k][i]);
  }

        /*把matrix[i][i]下的元素为0*/
        for(k=i+1;k<line;k++)
        {
            temp=matrix[k][i];
            for(j=i;j<row;j++)
            {
                matrix[k][j] -= (temp/matrix[i][i])*matrix[i][j];
            }

        }
    }
    for(i=0;i<line;i++)
    if(matrix[i][i]==0) break;
    return i;
}

double result(double matrix[MAX_MATRIX][MAX_MATRIX],int line)
{
    int i;
    double result=1;
    order(matrix,line,line);
    for(i=0;i<line;i++)
    result*=matrix[i][i];
    return result;
}

double resultMatrix(Matrix *a)
{
    return result(a->el,a->line);
}

Matrix *lineexchangeMatrix(Matrix *a,Matrix*b)
{
    multiMatrix(a,b);
    return a;
}

Matrix *rowexchangeMatrix(Matrix *a,Matrix*b)
{
    multiMatrix(a,b);
    return a;
}

Matrix *RMatrix(Matrix *b)
{
    orderMatrix(b);
    return b;
}

Matrix* LMatrix(Matrix *a)
{
    int i,j;
    orderMatrix(a);

    for(i=0;i<a->line/2;i++)
    {
        for(j=0;j<a->row;j++)
        {
            swap(&(a->el[i][j]),&(a->el[a->line-1-i][j]));
        }
    }

    for(i=0;i<a->row/2;i++)
    {
        for(j=0;j<a->line;j++)
        {
            swap(&(a->el[j][i]),&(a->el[j][a->row-1-i]));
        }
    }
    return a;
}

#endif