matrix-ops.c

   1 /*
   2 // Copyright (c) 2015 Intel Corporation
   3 //
   4 // Licensed under the Apache License, Version 2.0 (the "License");
   5 // you may not use this file except in compliance with the License.
   6 // You may obtain a copy of the License at
   7 //
   8 //      http://www.apache.org/licenses/LICENSE-2.0
   9 //
  10 // Unless required by applicable law or agreed to in writing, software
  11 // distributed under the License is distributed on an "AS IS" BASIS,
  12 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13 // See the License for the specific language governing permissions and
  14 // limitations under the License.
  15 */
  16
  17 #include "matrix-ops.h"
  18 #include <math.h>
  19
  20
  21 void transpose (int rows, int cols, double m[rows][cols], double m_trans[cols][rows])
  22 {
  23     int i,j;
  24
  25     for (i = 0; i < rows; i++)
  26         for (j = 0; j < cols; j++)
  27             m_trans[j][i] = m[i][j];
  28 }
  29
  30
  31 void multiply (int m, int n, int p, double m1[m][n], double m2[n][p], double result[m][p])
  32 {
  33     int i,j,k;
  34
  35     for (i = 0; i < m; i++)
  36         for (k = 0; k < p; k++) {
  37             result [i][k] = 0;
  38             for (j = 0; j < n; j++)
  39                 result [i][k] += m1[i][j] * m2 [j][k];
  40         }
  41 }
  42
  43
  44 void invert (int s, double m[s][s],  double m_inv[s][s])
  45 {
  46     double t;
  47     int swap,i,j,k;
  48     double tmp[s][s];
  49
  50     for (i = 0; i < s; i++)
  51         for (j = 0; j < s; j++)
  52             m_inv[i][j] = 0;
  53
  54     for (i = 0; i < s; i++)
  55         m_inv[i][i] = 1;
  56
  57     assign(s,s,m,tmp);
  58
  59     for (i = 0; i < s; i++) {
  60         swap = i;
  61         for (j = i+1; j < s; j++) {
  62             if (fabs(tmp[i][j]) > fabs(tmp[i][i]))
  63                 swap = j;
  64         }
  65
  66         if (swap != i) {
  67             /* swap rows */
  68             for (k = 0; k < s; k++) {
  69                 t = tmp[k][i];
  70                 tmp[k][i] = tmp[k][swap];
  71                 tmp[k][swap] = t;
  72
  73                 t = m_inv[k][i];
  74                 m_inv[k][i] = m_inv[k][swap];
  75                 m_inv[k][swap] = t;
  76             }
  77         }
  78
  79         t = 1 / tmp[i][i];
  80
  81         for (k = 0 ; k < s ; k++) {
  82             tmp[k][i] *= t;
  83             m_inv[k][i] *= t;
  84         }
  85
  86         for (j = 0 ; j < s ; j++)
  87             if (j != i) {
  88                 t = tmp[i][j];
  89                 for (k = 0 ; k < s; k++) {
  90                     tmp[k][j] -= tmp[k][i] * t;
  91                     m_inv[k][j] -= m_inv[k][i] * t;
  92                 }
  93             }
  94     }
  95 }
  96
  97
  98 void multiply_scalar_inplace(int rows, int cols, double m[rows][cols], double scalar)
  99 {
 100     int i,j;
 101
 102     for (i = 0; i < rows; i++)
 103         for (j = 0; j < cols; j++)
 104             m[i][j] = m[i][j] * scalar;
 105 }
 106
 107
 108 void assign (int rows, int cols, double m[rows][cols], double m1[rows][cols])
 109 {
 110     int i,j;
 111
 112     for (i = 0; i < rows; i++)
 113         for (j = 0; j < cols; j++)
 114             m1[i][j] = m[i][j];
 115 }
 116
 117 void substract (int rows, int cols, double m1[rows][cols], double m2[rows][cols], double res[rows][cols])
 118 {
 119     int i,j;
 120
 121     for (i = 0; i < rows; i++)
 122         for (j = 0; j < cols; j++)
 123             res[i][j] = m1[i][j] - m2[i][j];
 124 }