#include "copyleft.h"

/*
    GEPASI - a simulator of metabolic pathways and other dynamical systems
    Copyright (C) 1989, 1992  Pedro Mendes
*/

/*************************************/
/*                                   */
/*          matrix operations        */
/*                                   */
/*          MICROSOFT C 6.00         */
/*           QuickC/WIN 1.0          */
/*             ULTRIX cc             */
/*              GNU gcc              */
/*                                   */
/*   (include here compilers that    */
/*   compiled GEPASI successfully)   */
/*                                   */
/*************************************/


#include <stdio.h>
#include <string.h>
#include <math.h>
#include "globals.h"
#include "globvar.h"

#define NO_MTX_ERROR 0
#define NO_MULT 1
#define NON_INVERT 2
#define OUT_OF_MEM 3
#define NEARLY_ZERO 1e-10
#define TINY 1.0e-20;

/* reset a vector to zero */
void zerovct( double m1[MAX_MET] )
{
 register int i;

 for(i=0; i<MAX_MET; i++)
   m1[i] = 0.0;
}

/* add two vectors */
void addvct( double m1[MAX_MET], double m2[MAX_MET],
             double m3[MAX_MET], int r               )
{
 register int i;

 for(i=0; i<r; i++)
   m3[i] = m1[i] + m2[i];
}

/* add two matrixes */
void addmtx( double m1[], double m2[], double m3[], int r, int c )
{
 register int i,d;

 d = r * c;
 for(i=0; i<d; i++) m3[i]=m1[i]+m2[i];
}

/* add one matrix with a scalar (double) multiple of another */
void addmtxsm( double m1[MAX_MET][MAX_MET], double m2[MAX_MET][MAX_MET],
          double s,                    double m3[MAX_MET][MAX_MET],
          int r, int c                                              )
{
 register int i,j;

 for(i=0; i<r; i++)
  for(j=0;j<c;j++)
   m3[i][j] = m1[i][j] + s * m2[i][j];
}

/* add one vector with a scalar (double) multiple of another */
void addvctsm( double m1[MAX_MET], double m2[MAX_MET],
               double s,           double m3[MAX_MET],
               int r                                   )
{
 register int i;

 for(i=0; i<r; i++)
   m3[i] = m1[i] + s * m2[i];
}

/* (s1 * v1) + (s2 * v2)  - vectors */
void addvctsm2( double m1[MAX_MET], double s1,
                double m2[MAX_MET], double s2,
                double m3[MAX_MET], int r       )
{
 register int i;

 for(i=0; i<r; i++)
   m3[i] = s1 * m1[i] + s2 * m2[i];
}

/* multiply one matrix by another */
void multmtx( double m1[MAX_STEP][MAX_MET], float m2[MAX_MET][MAX_MET],
             double m3[MAX_STEP][MAX_MET], int  r1, int c1, int r2, int c2 )
{
 register int j,k;
          int i;
 double acum;

/* if (c1!=r2) return(NO_MULT); disabled for speed */
 for(i=0;i<r1;i++)
  for(j=0;j<c2;j++)
  {
   acum=0.0;
   for(k=0;k<c1;k++) acum += m1[i][k] * (double) m2[k][j];
   m3[i][j]=acum;
  }
/*  return(NO_MTX_ERROR); disabled for speed */
}

/* scalar (double) multiple of multiplication of one matrix by a vector */
void multmtxvsm( double m1[MAX_MET][MAX_MET], double m2[MAX_MET],
                 double s,                    double m3[MAX_MET],
                 int r1, int c1, int r2                           )
{
 register int i,k;
 double acum;

/* if (c1!=r2) return(NO_MULT); disabled for speed */
 for(i=0;i<r1;i++)
 {
  acum=0.0;
  for(k=0;k<c1;k++) acum+=m1[i][k] * m2[k];
  m3[i] = s * acum;
 }
/*  return(NO_MTX_ERROR); disabled for speed */
}

/* multiplication of one matrix by a vector */
void multmtxv( double m1[MAX_MET][MAX_MET], double m2[MAX_MET],
               double m3[MAX_MET], int r1, int c1, int r2       )
{
 register int i,k;

/* if (c1!=r2) return(NO_MULT); disabled for speed */
 for(i=0;i<r1;i++)
 {
  m3[i] = 0.0;
  for(k=0;k<c1;k++) m3[i] += m1[i][k] * m2[k];
 }
/*  return(NO_MTX_ERROR); disabled for speed */
}

/* mutiply a matrix by a scalar (double) */
void multmtxs( double m1[MAX_MET][MAX_MET], double s,
               double m2[MAX_MET][MAX_MET], int r, int c )
{
 register int i,j;

 for(i=0;i<r;i++)
  for(j=0;j<c;j++)
   m2[i][j] = m1[i][j] * s;
}

void back_sub(double (*a)[MAX_MET][MAX_MET], int n, int *indx, double *b)
{
 int i, ii, ip, j;
 double sum;

 for( i=0, ii=-1; i<n; i++ )
 {
  ip = indx[i];
  sum = b[ip];
  b[ip] = b[i];
  if( ii != -1 )
   for( j=ii; j<i; j++ ) sum -= (*a)[i][j] * b[j];
  else if( sum ) ii = i;
  b[i] = sum;
 }
 for( i=n-1; i>=0; i-- )
 {
  sum = b[i];
  for( j=i+1; j<n; j++ )
   sum -= (*a)[i][j] * b[j];
  b[i] = sum/(*a)[i][i];
 }
}


int lu_decomp(double (*a)[MAX_MET][MAX_MET], int n, int *indx, double *d)
{
 int i, imax, j, k;
 double big, dum, sum, temp;
 double vv[MAX_MET];

 *d = (double) 1;
 for( i=0; i<n; i++ )
 {
  big = (double) 0;
  for( j=0; j<n; j++ )
   if( ( temp = fabs( (*a)[i][j] ) ) > big ) big = temp;
  if( big == (double) 0 )
  {
   return NON_INVERT;
  }
  vv[i] = (double) 1/big;
 }
 for( j=0; j<n; j++ )
 {
  for( i=0; i<j; i++ )
  {
   sum = (*a)[i][j];
   for( k=0; k<i; k++ ) sum -= (*a)[i][k] * (*a)[k][j];
   (*a)[i][j] = sum;
  }
  big = (double) 0;
  for( i=j; i<n; i++ )
  {
   sum = (*a)[i][j];
   for( k=0; k<j; k++ )
    sum -= (*a)[i][k] * (*a)[k][j];
   (*a)[i][j] = sum;
   if( ( dum = vv[i]*fabs( sum ) ) >= big )
   {
    big = dum;
    imax = i;
   }
  }
  if( j != imax )
  {
   for( k=0; k<n; k++ )
   {
    dum = (*a)[imax][k];
    (*a)[imax][k] = (*a)[j][k];
    (*a)[j][k] = dum;
   }
   *d = -(*d);
   vv[imax] = vv[j];
  }
  indx[j] = imax;
  if( (*a)[j][j] == (double) 0 ) (*a)[j][j] = TINY;
  if( j != n-1 )
  {
   dum = (double) 1/((*a)[j][j]);
   for( i=j+1; i<n; i++ ) (*a)[i][j] *= dum;
  }
 }
}


int lu_inverse( double (*m1)[MAX_MET][MAX_MET], double (*m2)[MAX_MET][MAX_MET], int n )
{
 double dt;
 register int i,j;
 double tvv[MAX_MET];
 int emptyrow, mat_indx[MAX_MET];				/* indexes for matrix row permut.	*/

 for( i=0; i<n; i++)
 {
  emptyrow = 1;
  for( j=0; j<n; j++)
   if ( (*m1)[i][j] != (double) 0 ) emptyrow = 0;
  if ( emptyrow ) return NON_INVERT;			/* if one row is empty m1 singular	*/
 }

 if ( !( j = lu_decomp(m1, n, mat_indx, &dt) ) ) return j;
 for( j=0; j<n; j++ )
 {
  for( i=0; i<n; i++ )
   tvv[i] = i==j ? (double) 1 : (double) 0;
  back_sub(m1, n, mat_indx, tvv);
  for( i=0; i<n; i++ ) (*m2)[i][j] = tvv[i];
 }
 return NO_MTX_ERROR;
}