home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / MATRIX.C < prev    next >
C/C++ Source or Header  |  1993-08-26  |  6KB  |  272 lines
  1. #include "copyleft.h"
  2.  
  3. /*
  4.     GEPASI - a simulator of metabolic pathways and other dynamical systems
  5.     Copyright (C) 1989, 1992  Pedro Mendes
  6. */
  7.  
  8. /*************************************/
  9. /*                                   */
  10. /*          matrix operations        */
  11. /*                                   */
  12. /*          MICROSOFT C 6.00         */
  13. /*           QuickC/WIN 1.0          */
  14. /*             ULTRIX cc             */
  15. /*              GNU gcc              */
  16. /*                                   */
  17. /*   (include here compilers that    */
  18. /*   compiled GEPASI successfully)   */
  19. /*                                   */
  20. /*************************************/
  21.  
  22.  
  23. #include <stdio.h>
  24. #include <string.h>
  25. #include <math.h>
  26. #include "globals.h"
  27. #include "globvar.h"
  28.  
  29. #define NO_MTX_ERROR 0
  30. #define NO_MULT 1
  31. #define NON_INVERT 2
  32. #define OUT_OF_MEM 3
  33. #define NEARLY_ZERO 1e-10
  34. #define TINY 1.0e-20;
  35.  
  36. /* reset a vector to zero */
  37. void zerovct( double m1[MAX_MET] )
  38. {
  39.  register int i;
  40.  
  41.  for(i=0; i<MAX_MET; i++)
  42.    m1[i] = 0.0;
  43. }
  44.  
  45. /* add two vectors */
  46. void addvct( double m1[MAX_MET], double m2[MAX_MET],
  47.              double m3[MAX_MET], int r               )
  48. {
  49.  register int i;
  50.  
  51.  for(i=0; i<r; i++)
  52.    m3[i] = m1[i] + m2[i];
  53. }
  54.  
  55. /* add two matrixes */
  56. void addmtx( double m1[], double m2[], double m3[], int r, int c )
  57. {
  58.  register int i,d;
  59.  
  60.  d = r * c;
  61.  for(i=0; i<d; i++) m3[i]=m1[i]+m2[i];
  62. }
  63.  
  64. /* add one matrix with a scalar (double) multiple of another */
  65. void addmtxsm( double m1[MAX_MET][MAX_MET], double m2[MAX_MET][MAX_MET],
  66.           double s,                    double m3[MAX_MET][MAX_MET],
  67.           int r, int c                                              )
  68. {
  69.  register int i,j;
  70.  
  71.  for(i=0; i<r; i++)
  72.   for(j=0;j<c;j++)
  73.    m3[i][j] = m1[i][j] + s * m2[i][j];
  74. }
  75.  
  76. /* add one vector with a scalar (double) multiple of another */
  77. void addvctsm( double m1[MAX_MET], double m2[MAX_MET],
  78.                double s,           double m3[MAX_MET],
  79.                int r                                   )
  80. {
  81.  register int i;
  82.  
  83.  for(i=0; i<r; i++)
  84.    m3[i] = m1[i] + s * m2[i];
  85. }
  86.  
  87. /* (s1 * v1) + (s2 * v2)  - vectors */
  88. void addvctsm2( double m1[MAX_MET], double s1,
  89.                 double m2[MAX_MET], double s2,
  90.                 double m3[MAX_MET], int r       )
  91. {
  92.  register int i;
  93.  
  94.  for(i=0; i<r; i++)
  95.    m3[i] = s1 * m1[i] + s2 * m2[i];
  96. }
  97.  
  98. /* multiply one matrix by another */
  99. void multmtx( double m1[MAX_STEP][MAX_MET], float m2[MAX_MET][MAX_MET],
  100.              double m3[MAX_STEP][MAX_MET], int  r1, int c1, int r2, int c2 )
  101. {
  102.  register int j,k;
  103.           int i;
  104.  double acum;
  105.  
  106. /* if (c1!=r2) return(NO_MULT); disabled for speed */
  107.  for(i=0;i<r1;i++)
  108.   for(j=0;j<c2;j++)
  109.   {
  110.    acum=0.0;
  111.    for(k=0;k<c1;k++) acum += m1[i][k] * (double) m2[k][j];
  112.    m3[i][j]=acum;
  113.   }
  114. /*  return(NO_MTX_ERROR); disabled for speed */
  115. }
  116.  
  117. /* scalar (double) multiple of multiplication of one matrix by a vector */
  118. void multmtxvsm( double m1[MAX_MET][MAX_MET], double m2[MAX_MET],
  119.                  double s,                    double m3[MAX_MET],
  120.                  int r1, int c1, int r2                           )
  121. {
  122.  register int i,k;
  123.  double acum;
  124.  
  125. /* if (c1!=r2) return(NO_MULT); disabled for speed */
  126.  for(i=0;i<r1;i++)
  127.  {
  128.   acum=0.0;
  129.   for(k=0;k<c1;k++) acum+=m1[i][k] * m2[k];
  130.   m3[i] = s * acum;
  131.  }
  132. /*  return(NO_MTX_ERROR); disabled for speed */
  133. }
  134.  
  135. /* multiplication of one matrix by a vector */
  136. void multmtxv( double m1[MAX_MET][MAX_MET], double m2[MAX_MET],
  137.                double m3[MAX_MET], int r1, int c1, int r2       )
  138. {
  139.  register int i,k;
  140.  
  141. /* if (c1!=r2) return(NO_MULT); disabled for speed */
  142.  for(i=0;i<r1;i++)
  143.  {
  144.   m3[i] = 0.0;
  145.   for(k=0;k<c1;k++) m3[i] += m1[i][k] * m2[k];
  146.  }
  147. /*  return(NO_MTX_ERROR); disabled for speed */
  148. }
  149.  
  150. /* mutiply a matrix by a scalar (double) */
  151. void multmtxs( double m1[MAX_MET][MAX_MET], double s,
  152.                double m2[MAX_MET][MAX_MET], int r, int c )
  153. {
  154.  register int i,j;
  155.  
  156.  for(i=0;i<r;i++)
  157.   for(j=0;j<c;j++)
  158.    m2[i][j] = m1[i][j] * s;
  159. }
  160.  
  161. void back_sub(double (*a)[MAX_MET][MAX_MET], int n, int *indx, double *b)
  162. {
  163.  int i, ii, ip, j;
  164.  double sum;
  165.  
  166.  for( i=0, ii=-1; i<n; i++ )
  167.  {
  168.   ip = indx[i];
  169.   sum = b[ip];
  170.   b[ip] = b[i];
  171.   if( ii != -1 )
  172.    for( j=ii; j<i; j++ ) sum -= (*a)[i][j] * b[j];
  173.   else if( sum ) ii = i;
  174.   b[i] = sum;
  175.  }
  176.  for( i=n-1; i>=0; i-- )
  177.  {
  178.   sum = b[i];
  179.   for( j=i+1; j<n; j++ )
  180.    sum -= (*a)[i][j] * b[j];
  181.   b[i] = sum/(*a)[i][i];
  182.  }
  183. }
  184.  
  185.  
  186. int lu_decomp(double (*a)[MAX_MET][MAX_MET], int n, int *indx, double *d)
  187. {
  188.  int i, imax, j, k;
  189.  double big, dum, sum, temp;
  190.  double vv[MAX_MET];
  191.  
  192.  *d = (double) 1;
  193.  for( i=0; i<n; i++ )
  194.  {
  195.   big = (double) 0;
  196.   for( j=0; j<n; j++ )
  197.    if( ( temp = fabs( (*a)[i][j] ) ) > big ) big = temp;
  198.   if( big == (double) 0 )
  199.   {
  200.    return NON_INVERT;
  201.   }
  202.   vv[i] = (double) 1/big;
  203.  }
  204.  for( j=0; j<n; j++ )
  205.  {
  206.   for( i=0; i<j; i++ )
  207.   {
  208.    sum = (*a)[i][j];
  209.    for( k=0; k<i; k++ ) sum -= (*a)[i][k] * (*a)[k][j];
  210.    (*a)[i][j] = sum;
  211.   }
  212.   big = (double) 0;
  213.   for( i=j; i<n; i++ )
  214.   {
  215.    sum = (*a)[i][j];
  216.    for( k=0; k<j; k++ )
  217.     sum -= (*a)[i][k] * (*a)[k][j];
  218.    (*a)[i][j] = sum;
  219.    if( ( dum = vv[i]*fabs( sum ) ) >= big )
  220.    {
  221.     big = dum;
  222.     imax = i;
  223.    }
  224.   }
  225.   if( j != imax )
  226.   {
  227.    for( k=0; k<n; k++ )
  228.    {
  229.     dum = (*a)[imax][k];
  230.     (*a)[imax][k] = (*a)[j][k];
  231.     (*a)[j][k] = dum;
  232.    }
  233.    *d = -(*d);
  234.    vv[imax] = vv[j];
  235.   }
  236.   indx[j] = imax;
  237.   if( (*a)[j][j] == (double) 0 ) (*a)[j][j] = TINY;
  238.   if( j != n-1 )
  239.   {
  240.    dum = (double) 1/((*a)[j][j]);
  241.    for( i=j+1; i<n; i++ ) (*a)[i][j] *= dum;
  242.   }
  243.  }
  244. }
  245.  
  246.  
  247. int lu_inverse( double (*m1)[MAX_MET][MAX_MET], double (*m2)[MAX_MET][MAX_MET], int n )
  248. {
  249.  double dt;
  250.  register int i,j;
  251.  double tvv[MAX_MET];
  252.  int emptyrow, mat_indx[MAX_MET];                /* indexes for matrix row permut.    */
  253.  
  254.  for( i=0; i<n; i++)
  255.  {
  256.   emptyrow = 1;
  257.   for( j=0; j<n; j++)
  258.    if ( (*m1)[i][j] != (double) 0 ) emptyrow = 0;
  259.   if ( emptyrow ) return NON_INVERT;            /* if one row is empty m1 singular    */
  260.  }
  261.  
  262.  if ( !( j = lu_decomp(m1, n, mat_indx, &dt) ) ) return j;
  263.  for( j=0; j<n; j++ )
  264.  {
  265.   for( i=0; i<n; i++ )
  266.    tvv[i] = i==j ? (double) 1 : (double) 0;
  267.   back_sub(m1, n, mat_indx, tvv);
  268.   for( i=0; i<n; i++ ) (*m2)[i][j] = tvv[i];
  269.  }
  270.  return NO_MTX_ERROR;
  271. }
  272.