math_util.c

multipleshooting to find the roots of the non-linear functions

  long double vec_mul_btc( b, c, n )
     long double *b;
     long double *c;
     int n;
     {
        int i;
        long double a;
  
        a = 0.0;
        for( i = 0; i < n; i++ )
              a += b[i]*c[i];
              
        return( a );
     }

  /* ----------------------------------------------------------------------------
  
  VEC_ADD
  
      This function adds array b to c and puts the result in a. Both
      arrays are long double of length n.
  
                               {a} = {b} + {c}
  
  -----------------------------------------------------------------------------*/
  void vec_add( a, b, c, n )
     long double *a;
     long double *b;
     long double *c;
     int n;
     {
        int i;
  
        for( i = 0; i < n; i++ )
           *(a+i) = *(b+i) + *(c+i);
  
     }
  /* ----------------------------------------------------------------------------
  
  ENORM
  
      This function calculates the Euclidian norm of a vector of length n.
  
                           norm = sqrt( {a} * Transpose{a}  )
  
  -----------------------------------------------------------------------------*/
  long double enorm( a, n )
     long double *a;
     int n;
     {
        long double sum;
        double sqrt( double );
        int i;
  
        sum  = 0.0;
  
        for( i = 0; i < n; i++ )
           sum += a[i] * a[i];
  
        if( sum > 0.0 )
           sum = sqrt((double) sum );
  
        return( sum );
     }
  /* ----------------------------------------------------------------------------
  
  MAT_COPY
  
     This function copies matrix a into matrix b. The dimensions of a and b are
     m by n.
  
                         [a] = [b]
  
  -----------------------------------------------------------------------------*/
  void mat_copy( a, b, m, n )
     long double **a;
     long double **b;
     int m;
     int n;
     {
        int i;
        int j;
  
        for( i = 0; i < m; i++ )
           for( j = 0; j < n; j++ )
              a[i][j] = b[i][j];
  
     }
  /* ----------------------------------------------------------------------------
  
  MAT_SUB
  
      This function substracts matrix b from c and puts the result in a. The
      dimension of a, b and c are m by n.
  
                               [a] = [b] - [c]
  
  -----------------------------------------------------------------------------*/
  void mat_sub( a, b, c, m, n )
     long double **a;
     long double **b;
     long double **c;
     int m;
     int n;
     {
        int i;
        int j;
  
        for( i = 0; i < m; i++ )
           for( j = 0; j < n; j++ )
              a[i][j] = b[i][j] - c[i][j];
  
     }
  /* ----------------------------------------------------------------------------
  
  MAT_ADD
  
      This function adds matrix b to c and puts the result in a. The
      dimension of a, b and c are m by n.
  
                              [a] = [b] + [c]
  
  -----------------------------------------------------------------------------*/
  void mat_add( a, b, c, m, n )
     long double **a;
     long double **b;
     long double **c;
     int m;
     int n;
     {
        int i;
        int j;
  
        for( i = 0; i < m; i++ )
           for( j = 0; j < n; j++ )
              a[i][j] = b[i][j] + c[i][j];
  
     }
  /* ----------------------------------------------------------------------------
  
  MAT_MUL
  
      This function multiplies matrix b with c and puts the result in a. The
      dimension of a is l by n; b is l by m; and c is m by n.
  
                               [a] = [b][c]
  
  -----------------------------------------------------------------------------*/
  void mat_mul( a, b, c, l, m, n )
     long double **a; /* l by n matrix - output */
     long double **b; /* l by m matrix - input  */
     long double **c; /* m by n matrix - input  */
     int l;
     int m;
     int n;
     {
        int i;
        int j;
        int k;
        long double sum;
  
        for( i = 0; i < l; i++ )
           for( j = 0; j < n; j++ ) {
  
              sum = 0.0;
  
              for( k = 0; k < m; k++ )
                 sum += b[i][k] * c[k][j];
  
              a[i][j] = sum;
           }
  
     }
  /* ----------------------------------------------------------------------------
  
  MAT_MULV
  
      This function multiplies matrix b with vector c and puts the result in
      vector a. The length of a is m; b is m by n; and c is n by 1.
  
                               {a} = [b]{c}
  
  -----------------------------------------------------------------------------*/
  void mat_mulv( a, b, c, m, n )
     long double *a;
     long double **b;
     long double *c;
     int m;
     int n;
     {
        int i;
        int j;
  
        for( i = 0; i < m; i++ ) {
           a[i] = 0.0;
           for( j = 0; j < n; j++ )
              a[i] += b[i][j] * c[j];
        }
  
     }
  /* ----------------------------------------------------------------------------
  
  MAT_TRAN_BC
  
      This function multiplies matrix transpose b with c and puts the result in 
      a.
                                        T
                               [a] = [b] [c]
  
  -----------------------------------------------------------------------------*/
  void mat_tran_bc( a, b, c, l, m, n )
     long double **a; /* m by n matrix - output */
     long double **b; /* l by m matrix - input  */
     long double **c; /* l by n matrix - input  */
     int l;
     int m;
     int n;
     {
        int i;
        int j;
        int k;
        long double sum;
  
        for( i = 0; i < m; i++ )
           for( j = 0; j < n; j++ ) {
  
              sum = 0.0;
  
              for( k = 0; k < l; k++ )
                 sum += b[k][i] * c[k][j];
  
              a[i][j] = sum;
           }
     }
  /* ----------------------------------------------------------------------------
  
  MAT_BTRAN_C
  
      This function multiplies matrix b with transpose c and puts the result in 
      a.
                                           T
                               [a] = [b][c]
  
  -----------------------------------------------------------------------------*/
  void mat_btran_c( a, b, c, l, m, n )
     long double **a; /* l by n matrix - output */
     long double **b; /* l by m matrix - input  */
     long double **c; /* n by m matrix - input  */
     int l;
     int m;
     int n;
     {
        int i;
        int j;
        int k;
        long double sum;
  
        for( i = 0; i < l; i++ )
           for( j = 0; j < n; j++ ) {
  
              sum = 0.0;
  
              for( k = 0; k < m; k++ )
                 sum += b[i][k] * c[j][k];
  
              a[i][j] = sum;
           }
     }
  /* ----------------------------------------------------------------------------
  
  MAT_TRAN
  
      This function will transpose matrix b and put the result into a.
  
                                        T
                               [a] = [b] 
  
  -----------------------------------------------------------------------------*/
  void mat_tran( a, b, l, m )
     long double **a; /* m by l matrix - output */
     long double **b; /* l by m matrix - input  */
     int l;
     int m;
     {
        int i;
        int j;
  
        for( i = 0; i < l; i++ )
           for( j = 0; j < m; j++ )
              a[j][i] = b[i][j];
     }
  /* ----------------------------------------------------------------------------
  
  MAT_NORM_A
  
      This function returns a bound on the norm of a square matrix. The bound is
      defined as follows:
      
      
      norm( [a] ) = sum( abs( [a(i,j)] ) )
  
  -----------------------------------------------------------------------------*/
long double mat_norm_a( a, n )
long double **a;
int n;
   {
      int i, j;
      long double sum;
      
      sum = 0.0;
      for( i = 0; i < n; i++ )
         for( j = 0; j < n; j++ )
            sum += fabs((double) a[i][j]);
            
      return( sum );
   }
  /* ----------------------------------------------------------------------------
  
  MAT_NORM_B
  
      This function returns a bound on the norm of a square matrix. The bound is
      defined as follows:
      
      
      norm( [a] ) = sqrt( sum( [a(i,j)]**2 ) )
  
  -----------------------------------------------------------------------------*/
long double mat_norm_b( a, n )
long double **a;
int n;
   {
      int i, j;
      long double sum;
      
      sum = 0.0;
      for( i = 0; i < n; i++ )
         for( j = 0; j < n; j++ )
            sum += a[i][j]*a[i][j];
      return( (long double) sqrt( (double) sum ) );
   }
  /* ----------------------------------------------------------------------------
  
  MAT_NORM_C
  
      This function returns a bound on the norm of a square matrix. The bound is
      defined as follows:
      
      
      norm( [a] ) = max( sum( abs( [a(i,j)] ) ) )
                     i    j
  -----------------------------------------------------------------------------*/
long double mat_norm_c( a, n )
long double **a;
int n;
   {
      int i, j;
      long double sum, max;
      
      max = 0.0;
      for( i = 0; i < n; i++ ) {
         sum = 0.0;
         for( j = 0; j < n; j++ )
            sum += fabs((double) a[i][j]);
         if( sum > max )
            max = sum;
      }
      return( max );     
   }
 /* ----------------------------------------------------------------------------
  
  MAT_MULS
  
      This function multiplies matrix b by a scalar c and puts the result in
      matrix a. The dimension of a is n by m; b is n by m.
  
                               [a] = [b]*c
  
  -----------------------------------------------------------------------------*/
  void mat_muls( a, b, c, n, m )
     long double **a;
     long double **b;
     long double c;
     int m;
     int n;
     {
        int i;
        int j;
  
        for( i = 0; i < n; i++ ) 
           for( j = 0; j < m; j++ )
              a[i][j] = b[i][j] * c;
        }
/* ----------------------------------------------------------------------------
   MAT_SOLVE
   
   This function returns the solution to the matrix equation
   
   AX = B
      
   INPUT DATA:
               A - an n by n matrix of long double.
               B - an n by m matrix of long double.
               n - matrix dimension, int.
               m - matrix dimension, int.
   	
   OUTPUT DATA:
               X - an n by m matrix of long double. The solution to the equation.
   				
   RETURN CODES:
   	       0 - normal completion.
           1 - n <= 0 or m <= 0.
           2 - Unable to allocate storage.
	       3 - The A matrix is singular.
 ---------------------------------------------------------------------------- */ 
int mat_solve( x, a, b, n, m )
long double **x, **a, **b;
int n, m;
   {
      int i, j, err;      
      int *ipivot;
      long double *bvec, *xvec;
           
      void a1d_free_int( int * );
      int *a1d_allo_int( int );      
      int factor( long double **, int, int * );
      void solve( long double **, int *, long double *, int, long double * ); 
      void vec_column( long double *, long double **, int, int );
      void a1d_free_dbl( long double * );
      long double *a1d_allo_dbl( int );
           
/*
   Do some error checks
*/
      if( n <= 0 || m <= 0 )
         return( 1 );

      ipivot = a1d_allo_int( n );
      xvec = a1d_allo_dbl( n );
      bvec = a1d_allo_dbl( n );
       
      if( ipivot == NULL || bvec == NULL || xvec == NULL ) {
         printf(" Unable to allocate storage for ipivot, bvec and xvec in mat_solve\n");
         return( 2 );
      }
	 
      err = 0;
      	 
      err = factor( a, n, ipivot );
      
      if( err != 1 && err != -1 ) {
         printf(" Error -> %d, from factor in mat_solve\n",err);
         if( err == 0 )
	    return( 3 );
	    
	 return( 2 );
      }