matrix-ops.c

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