Breaking News

Matrix operations in C

11:08 AM

this header file contains well structured functions for matrix creation ,and operation.this is the initial version of matrix.h.

/*
 =======================================================================
 Name        : Matrix.h
 Author      : Muhammad
 Version     : 1.0
 Copyright   : Your copyright notice
 Description : Matrix operations in C, Ansi-style
 =======================================================================
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#ifndef __MATRIX_H
#define __MATRIX_H
#define TALLOC(n,typ) (typ *)malloc(n*sizeof(typ))
#define min(x,y) (x > y)?y:x
#define EL_Type double
#ifndef true
#define true  1==1
#define false 1!=1
#endif
typedef unsigned char bool;

typedef struct sMtx {
      int dim_x, dim_y;
      EL_Type *m_stor;
      EL_Type **mtx;
}*Matrix, sMatrix;

typedef struct sRvec {
      int dim_x;
      EL_Type *m_stor;
}*RowVec, sRowVec;

typedef struct LUmtx {
      Matrix LU;
      int m, n, pivsign;
      RowVec piv;
}*LUDecomposition, sLUDecomposition;

#define MtxGet( m, rix, cix ) m->mtx[rix][cix]
#define LUGet( m, rix, cix ) m->LU->mtx[rix][cix]

void LUSet(LUDecomposition m, int rix, int cix, EL_Type val) {
      m->LU->mtx[rix][cix]= val;
}

RowVec NewRowVec(int n) {
      RowVec r;
      r= TALLOC(1,sRowVec);
      r->dim_x = n;
      r->m_stor = TALLOC(n,EL_Type);
      return r;
}

RowVec GetRow(Matrix m, int row) {
      RowVec r = NewRowVec(m->dim_x);
      int i = 0;
      for (; i<m->dim_x; i++)
            r->m_stor[i]=MtxGet(m,row,i);
      return r;
}

Matrix NewMatrix(int x_dim, int y_dim) {
      int n;
      Matrix m;
      m = TALLOC( 1, sMatrix);
      n = x_dim * y_dim;
      m->dim_x = x_dim;
      m->dim_y = y_dim;
      m->m_stor = TALLOC(n, EL_Type);
      m->mtx = TALLOC(m->dim_y, EL_Type *);
      for (n=0; n<y_dim; n++) {
            m->mtx[n] = m->m_stor+n*x_dim;
      }
      return m;
}
Matrix NewMatrixFromVec(int x_dim, int y_dim, EL_Type* data) {
      int n;
      Matrix m;
      m = TALLOC( 1, sMatrix);
      n = x_dim * y_dim;
      m->dim_x = x_dim;
      m->dim_y = y_dim;
      m->m_stor = data;
      m->mtx = TALLOC(m->dim_y, EL_Type *);
      for (n=0; n<y_dim; n++) {
            m->mtx[n] = m->m_stor+n*x_dim;
      }
      return m;
}

Matrix CopyMatrix(Matrix m) {
      Matrix new_matrix = NewMatrix(m->dim_x, m->dim_y);
      int i, j;
      for (i = 0; i < m->dim_y; i++)
            for (j = 0; j<m->dim_x; j++) {
                  new_matrix->mtx[i][j] = MtxGet(m,i,j);
            }
      return new_matrix;
}

void MtxSetRow(Matrix m, int irow, EL_Type *v) {
      int ix;
      EL_Type *mr;
      mr = m->mtx[irow];
      for (ix=0; ix<m->dim_x; ix++)
            mr[ix] = v[ix];
}

Matrix InitMatrix(int x_dim, int y_dim, EL_Type **v) {
      Matrix m;
      int iy;
      m = NewMatrix(x_dim, y_dim);
      for (iy=0; iy<y_dim; iy++)
            MtxSetRow(m, iy, v[iy]);
      return m;
}
/**
 * diplay matrix
 */
void MtxDisplay(Matrix m) {
      int iy, ix;
      const char *sc;
      for (iy=0; iy<m->dim_y; iy++) {
            printf("   ");
            sc = " ";
            for (ix=0; ix<m->dim_x; ix++) {
                  printf("%s %3.3f", sc, m->mtx[iy][ix]);
                  sc = ",";
            }
            printf("\n");
      }
      printf("\n");
}

void MtxMulAndAddRows(Matrix m, int ixrdest, int ixrsrc, EL_Type mplr) {
      int ix;
      EL_Type *drow, *srow;
      drow = m->mtx[ixrdest];
      srow = m->mtx[ixrsrc];
      for (ix=0; ix<m->dim_x; ix++)
            drow[ix] += mplr * srow[ix];
      //    printf("Mul row %d by %d and add to row %d\n", ixrsrc, mplr, ixrdest);
      //    MtxDisplay(m);
}

void MtxSwapRows(Matrix m, int rix1, int rix2) {
      EL_Type *r1, *r2, temp;
      int ix;
      if (rix1 == rix2)
            return;
      r1 = m->mtx[rix1];
      r2 = m->mtx[rix2];
      for (ix=0; ix<m->dim_x; ix++)
            temp = r1[ix];
      r1[ix]=r2[ix];
      r2[ix]=temp;
      //    printf("Swap rows %d and %d\n", rix1, rix2);
      //    MtxDisplay(m);
}

void MtxNormalizeRow(Matrix m, int rix, int lead) {
      int ix;
      EL_Type *drow;
      EL_Type lv;
      drow = m->mtx[rix];
      lv = drow[lead];
      for (ix=0; ix<m->dim_x; ix++)
            drow[ix] /= lv;
      //    printf("Normalize row %d\n", rix);
      //    MtxDisplay(m);
}

void MtxToReducedREForm(Matrix m) {
      int lead;
      int rix, iix;
      EL_Type lv;
      int rowCount = m->dim_y;

      lead = 0;
      for (rix=0; rix<rowCount; rix++) {
            if (lead >= m->dim_x)
                  return;
            iix = rix;
            while (0 == MtxGet(m, iix,lead)) {
                  iix++;
                  if (iix == rowCount) {
                        iix = rix;
                        lead++;
                        if (lead == m->dim_x)
                              return;
                  }
            }
            MtxSwapRows(m, iix, rix);
            MtxNormalizeRow(m, rix, lead);
            for (iix=0; iix<rowCount; iix++) {
                  if (iix != rix) {
                        lv = MtxGet(m, iix, lead );
                        MtxMulAndAddRows(m, iix, rix, -lv) ;
                  }
            }
            lead++;
      }
}
/**
 * Compute LU such that
 * A = L*U
 *
 */
LUDecomposition NewLUDecomposition(Matrix A) {

      // Use a "left-looking", dot-product, Crout/Doolittle algorithm.

      LUDecomposition LU= TALLOC( 1, sLUDecomposition);
      LU->LU = CopyMatrix(A);
      LU->m = A->dim_y;
      LU->n = A->dim_x;
      int m = LU->m;
      int n = LU->n;
      LU->piv = NewRowVec(LU->m);
      int i, j;
      for (i = 0; i < m; i++) {
            LU->piv->m_stor[i] = i;
      }
      LU->pivsign = 1;
      EL_Type* LUrowi;
      RowVec LUcolj = NewRowVec(m);

      // Outer loop.

      for (j = 0; j < n; j++) {

            // Make a copy of the j-th column to localize references.

            for (i = 0; i < m; i++) {
                  LUcolj->m_stor[i] = LUGet(LU,i,j);
            }

            // Apply previous transformations.

            for (i = 0; i < m; i++) {
                  LUrowi = LU->LU->mtx[i];
                  //               
                  // Most of the time is spent in the following dot product.

                  int kmax= min(i,j);
                  double s = 0.0;
                  int k;
                  for (k = 0; k < kmax; k++) {
                        //printf("LUrowi->m_stor[%d]=%f LUcolj->m_stor[%d]=%f\n",k,LUrowi[k],k,LUcolj->m_stor[k]);
                        s += LUrowi[k] * LUcolj->m_stor[k];
                  }

                  LUrowi[j] = LUcolj->m_stor[i] -= s;
                  //printf("s=%f and LUrowi->m_stor[%d]=%f LUcolj->m_stor[%d]=%f\n",s,j,LUrowi->m_stor[j],i,LUcolj->m_stor[i]);
            }

            // Find pivot and exchange if necessary.

            int p = j;
            for (i = j+1; i < m; i++) {
                  if (abs(LUcolj->m_stor[i]) > abs(LUcolj->m_stor[p])) {
                        p = i;

                  }
            }
            if (p != j) {
                  int k;
                  for (k = 0; k < n; k++) {
                        double t= LUGet(LU,p,k);
                        LUSet(LU, p, k, LUGet(LU,j,k));
                        LUSet(LU, j, k, t);
                  }
                  EL_Type val = LU->piv->m_stor[p];
                  LU->piv->m_stor[p] = LU->piv->m_stor[j];
                  LU->piv->m_stor[j] = val;
                  LU->pivsign = -LU->pivsign;
            }

            // Compute multipliers.

            if ((j < m) & (LUGet(LU,j,j) != 0)) {
                  for (i = j+1; i < m; i++) {
                        EL_Type v= LUGet(LU,i,j);
                        //LU[i][j] /= LU[j][j];
                        LUSet(LU, i, j, v/LUGet(LU,j,j));
                  }
            }
      }
      //printf("\nmy LU\n");
      //MtxDisplay(LU->LU);
      return LU;
}
/**
 * check if the diagonal is 0
 */
int isNonsingular(LUDecomposition LU) {
      int j;
      for (j = 0; j < LU->n; j++) {
            if (LUGet(LU,j,j) == 0)
                  return false;
      }
      return true;
}
/**
 * matrix multiplication
 * A*B = X
 */
Matrix times(Matrix A, Matrix B) {
      if (B->dim_y != A->dim_x) {
            printf("Matrix inner dimensions must agree function times.");
            return 0;
      }
      Matrix X = NewMatrix(B->dim_x, A->dim_y);

      RowVec Bcolj = NewRowVec(A->dim_x);
      int j, k, i;
      for (j = 0; j < B->dim_x; j++) {
            for (k = 0; k < B->dim_x; k++) {
                  Bcolj->m_stor[k] = B->mtx[k][j];
            }
            for (i = 0; i < A->dim_y; i++) {
                  RowVec Arowi = GetRow(A, i);
                  double s = 0;
                  for (k = 0; k < A->dim_x; k++) {
                        s += Arowi->m_stor[k]*Bcolj->m_stor[k];
                  }
                  X->mtx[i][j] = s;
            }
      }
      return X;
}
Matrix identity(int m, int n) {
      Matrix A = NewMatrix(m, n);
      int i, j;
      for (i = 0; i < m; i++) {
            for (j = 0; j < n; j++) {
                  A->mtx[i][j] = (i == j ? 1 : 0);
            }
      }
      return A;
}
/** Get a submatrix.
 @param r    Array of row indices.
 @param i0   Initial column index
 @param i1   Final column index
 @return     A(r(:),j0:j1)  
 */
Matrix getMatrix(Matrix A, RowVec r, int j0, int j1) {
      Matrix X = NewMatrix(j1-j0+1, r->dim_x);
      int i, j;

      for (i = 0; i < r->dim_x; i++) {
            for (j = j0; j <= j1; j++) {
                  X->mtx[i][j-j0] = A->mtx[(int)r->m_stor[i]][j];
            }
      }
      return X;
}
/**
 * A*X = B
 * we would like to find X
 * in case of B is identity
 * X is the inverse of A
 */
Matrix solve(LUDecomposition LU, Matrix B) {
      if (B->dim_y != LU->m) {
            printf("Error in dimentions in solve func");
            return 0;
      }

      if (!isNonsingular(LU)) {
            printf("the matrix is singular");
            return 0;
      }

      // Copy right hand side with pivoting
      int nx = B->dim_x;
      Matrix Xmat = getMatrix(B, LU->piv, 0, nx-1);

      // Solve L*Y = B(piv,:)
      int k, i, j;
      for (k = 0; k < LU->n; k++) {
            for (i = k+1; i < LU->n; i++) {
                  for (j = 0; j < nx; j++) {
                        Xmat->mtx[i][j] -= Xmat->mtx[k][j]* LUGet(LU,i,k);
                  }
            }
      }
      // Solve U*X = Y;
      for (k = LU->n-1; k >= 0; k--) {
            for (j = 0; j < nx; j++) {
                  Xmat->mtx[k][j] = MtxGet(Xmat,k,j)/ LUGet(LU,k,k);
            }
            for (i = 0; i < k; i++) {
                  for (j = 0; j < nx; j++) {
                        Xmat->mtx[i][j] -= MtxGet(Xmat,k,j)* LUGet(LU,i,k);
                  }
            }
      }
      return Xmat;
}
/**
 *
 */
Matrix inverse(Matrix A) {
      return solve(NewLUDecomposition(A), identity(A->dim_y, A->dim_y));
}

/**
 * returns new N*M matrix with random values
 */
Matrix random(int m, int n) {
      Matrix A = NewMatrix(m, n);
      int i, j;
      for (i = 0; i < m; i++) {
            for (j = 0; j < n; j++) {

                  A->mtx[i][j] = (EL_Type)rand();
            }
      }
      return A;
}
#endif

Thanks Alot For My Friend Muhammed Hamed 

Yours Moataz Muhammed

No comments

Subscribe to: Post Comments ( Atom )