matrix-ops.c

   1 #include "matrix-ops.h"
   2 #include <math.h>
   3 #include <stdio.h>
   4
   5 void transpose (int rows, int cols, double m[rows][cols], double m_trans[cols][rows])
   6 {
   7     int i,j;
   8
   9     for (i = 0; i < rows; i++)
  10         for (j = 0; j < cols; j++)
  11             m_trans[j][i] = m[i][j];
  12 }
  13
  14
  15 void multiply (int m, int n, int p, double m1[m][n], double m2[n][p], double result[m][p])
  16 {
  17     int i,j,k;
  18
  19     for (i = 0; i < m; i++)
  20         for (k = 0; k < p; k++) {
  21             result [i][k] = 0;
  22             for (j = 0; j < n; j++)
  23                 result [i][k] += m1[i][j] * m2 [j][k];
  24         }
  25 }
  26 void invert (int s, double m[s][s],  double m_inv[s][s])
  27 {
  28     double t;
  29     int swap,i,j,k;
  30     double tmp[s][s];
  31
  32     for (i = 0; i < s; i++)
  33         for (j = 0; j < s; j++)
  34             m_inv[i][j] = 0;
  35     for (i = 0; i < s; i++)
  36         m_inv[i][i] = 1;
  37
  38     assign(s,s,m,tmp);
  39
  40     for (i = 0; i < s; i++) {
  41         swap = i;
  42         for (j = i+1; j < s; j++) {
  43             if (fabs(tmp[i][j]) > fabs(tmp[i][i]))
  44                 swap = j;
  45         }
  46
  47         if (swap != i) {
  48             /* swap rows */
  49             for (k = 0; k < s; k++) {
  50                 t = tmp[k][i];
  51                 tmp[k][i] = tmp[k][swap];
  52                 tmp[k][swap] = t;
  53
  54                 t = m_inv[k][i];
  55                 m_inv[k][i] = m_inv[k][swap];
  56                 m_inv[k][swap] = t;
  57             }
  58         }
  59
  60         t = 1 / tmp[i][i];
  61         for (k = 0 ; k < s ; k++) {
  62             tmp[k][i] *= t;
  63             m_inv[k][i] *= t;
  64         }
  65
  66         for (j = 0 ; j < s ; j++) {
  67             if (j != i) {
  68                 t = tmp[i][j];
  69                 for (k = 0 ; k < s; k++) {
  70                     tmp[k][j] -= tmp[k][i] * t;
  71                     m_inv[k][j] -= m_inv[k][i] * t;
  72                 }
  73             }
  74         }
  75
  76     }
  77 }
  78 void multiply_scalar_inplace(int rows, int cols, double m[rows][cols], double scalar)
  79 {
  80     int i,j;
  81
  82     for (i = 0; i < rows; i++)
  83         for (j = 0; j < cols; j++)
  84             m[i][j] = m[i][j] * scalar;
  85 }
  86 void assign (int rows, int cols, double m[rows][cols], double m1[rows][cols])
  87 {
  88     int i,j;
  89
  90     for (i = 0; i < rows; i++)
  91         for (j = 0; j < cols; j++)
  92             m1[i][j] = m[i][j];
  93 }
  94
  95 void substract (int rows, int cols, double m1[rows][cols], double m2[rows][cols], double res[rows][cols])
  96 {
  97     int i,j;
  98
  99     for (i = 0; i < rows; i++)
 100         for (j = 0; j < cols; j++)
 101             res[i][j] = m1[i][j] - m2[i][j];
 102 }