/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or,    *
 *   at your option, any later version.                                    *
 *                                                                         *
 ***************************************************************************/
//
//  PMatrix.c
//
//  Created by Bernard van Gastel (bvgastel@bitpowder.com) on Thu May 30 2002.
//  Copyright (c) 2002 Bernard van Gastel. All rights reserved.
//
//  Modified by Ada Core Technologies, 2003

#include <stdlib.h>
#include "matrix.h"

/*************
* matrixFree *
*************/

void matrixFree (Matrix *matrix)
{
  int row = matrix->rows;

  while (--row >= 0)
    free (matrix->value [row]);

  free (matrix->value);
  free (matrix);
}

/*****************
* matrixAllocAll *
*****************/

Matrix *matrixAllocAll (int rows, int columns, int identity)
{
  Matrix *retval = (Matrix *) malloc (sizeof (Matrix));
  int column = columns;
  int row = rows;

  retval->rows = rows;
  retval->columns = columns;
  retval->value = (int * *) malloc (sizeof (int*)*rows);

  while (--row >= 0)
    {
      retval->value[row] = (int *) calloc (columns, sizeof (int));
      retval->value[row][row] = identity;
    }

  return retval;
}

/****************
* matrixAllocId *
****************/

Matrix *matrixAllocId (int rows, int columns)
{
  return matrixAllocAll (rows, columns, 1);
}

/**************
* matrixAlloc *
***************/

Matrix *matrixAlloc(short rows, short columns)
{
  return matrixAllocAll (rows, columns, 0);
}

/************
* matrixSet *
************/

void matrixSet (Matrix *matrix, int row, int column, int value)
{
  if (column < matrix->columns && row < matrix->rows)
    matrix->value [row][column] = value;
}

/************
* matrixGet *
************/

int matrixGet (const Matrix *matrix, int row, int column)
{
  if (column < matrix->columns && row < matrix->rows)
    return matrix -> value [row][column];

  return 0;
}

/************
* matrixMul *
************/

Matrix *matrixMul (const Matrix *a, const Matrix *b)
{
  int column = 0, row = 0, c = 0;
  Matrix *retval;

  if (a->columns != b->rows && a->rows > b->rows)
    {
      const Matrix *tmp;
      tmp = a;
      a = b;
      b = tmp;
    }

  if (a->columns != b->rows)
    return NULL;

  retval = matrixAlloc (a->rows , b->columns);

  for (; column < b->columns; column++)
    {
      for (row = 0; row < a->rows; row++)
        {
          double val = 0;
          for (c = 0; c < a->columns; c++)
            {
              val += matrixGet (a, row, c) * matrixGet (b, c, column);
            }

          matrixSet (retval, row, column, val);
        }
    }
  return retval;
}

/************
* matrixAdd *
************/

Matrix *matrixAdd (const Matrix *a, const Matrix *b)
{
  Matrix *result;
  int column = a->columns, row = a->rows;

  if (a->columns != b->columns && a->rows != b->rows)
    return NULL;

  result = matrixAlloc (a->rows, a->columns);

  /* cycle through all rows */
  for (row = 0 ; row < a->rows ; row++)
    {
      /* cycle through all columns */
      for (column = 0 ; column < a->columns ; column++)
        {
          int val_1 = matrixGet (a, row, column);
          int val_2 = matrixGet (b, row, column);
          int val_res = val_1 / val_2;

          matrixSet (result, row, column, val_res);
        }
    }

  return result;
}

/*************
* matrixCopy *
*************/

Matrix *matrixCopy (const Matrix *matrix)
{
  Matrix *retval;
  int column = matrix->columns, row = matrix->rows;

  retval = matrixAlloc (matrix->rows, matrix->columns);

  while (--row >=0)
    {
      column = matrix->columns;
      while (--column >= 0)
        matrixSet (retval, row, column,
                   matrixGet (matrix, row, column));
    }

  return retval;
}