matrix-ops.c

   1 /*
   2  * Copyright (C) 2014 Intel Corporation.
   3  */
   4
   5 #include "matrix-ops.h"
   6 #include <math.h>
   7 #include <stdio.h>
   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 void invert (int s, double m[s][s],  double m_inv[s][s])
  31 {
  32     double t;
  33     int swap,i,j,k;
  34     double tmp[s][s];
  35
  36     for (i = 0; i < s; i++)
  37         for (j = 0; j < s; j++)
  38             m_inv[i][j] = 0;
  39     for (i = 0; i < s; i++)
  40         m_inv[i][i] = 1;
  41
  42     assign(s,s,m,tmp);
  43
  44     for (i = 0; i < s; i++) {
  45         swap = i;
  46         for (j = i+1; j < s; j++) {
  47             if (fabs(tmp[i][j]) > fabs(tmp[i][i]))
  48                 swap = j;
  49         }
  50
  51         if (swap != i) {
  52             /* swap rows */
  53             for (k = 0; k < s; k++) {
  54                 t = tmp[k][i];
  55                 tmp[k][i] = tmp[k][swap];
  56                 tmp[k][swap] = t;
  57
  58                 t = m_inv[k][i];
  59                 m_inv[k][i] = m_inv[k][swap];
  60                 m_inv[k][swap] = t;
  61             }
  62         }
  63
  64         t = 1 / tmp[i][i];
  65         for (k = 0 ; k < s ; k++) {
  66             tmp[k][i] *= t;
  67             m_inv[k][i] *= t;
  68         }
  69
  70         for (j = 0 ; j < s ; j++) {
  71             if (j != i) {
  72                 t = tmp[i][j];
  73                 for (k = 0 ; k < s; k++) {
  74                     tmp[k][j] -= tmp[k][i] * t;
  75                     m_inv[k][j] -= m_inv[k][i] * t;
  76                 }
  77             }
  78         }
  79
  80     }
  81 }
  82 void multiply_scalar_inplace(int rows, int cols, double m[rows][cols], double scalar)
  83 {
  84     int i,j;
  85
  86     for (i = 0; i < rows; i++)
  87         for (j = 0; j < cols; j++)
  88             m[i][j] = m[i][j] * scalar;
  89 }
  90 void assign (int rows, int cols, double m[rows][cols], double m1[rows][cols])
  91 {
  92     int i,j;
  93
  94     for (i = 0; i < rows; i++)
  95         for (j = 0; j < cols; j++)
  96             m1[i][j] = m[i][j];
  97 }
  98
  99 void substract (int rows, int cols, double m1[rows][cols], double m2[rows][cols], double res[rows][cols])
 100 {
 101     int i,j;
 102
 103     for (i = 0; i < rows; i++)
 104         for (j = 0; j < cols; j++)
 105             res[i][j] = m1[i][j] - m2[i][j];
 106 }