[I2 logo] [RWTH logo] MOVES: Software Modeling and Verification
(Informatik 2)
Computer Science / RWTH / I2 / Research / Skil / Examples / Gauss
Printer-friendly
Skil Sample Applications

Gaussian Elimination



#include "Par/par.c"
#include "Array/array.h"


typedef struct _elemrec
        {
           float val ;
           int   row ;
           int   col ;
        }
        elemrec ;


#define rabs(x) (x >= 0 ? x : -x)               /* abs (x) for x float */


static void pck_f ($t v, Buff buff)
{
   return ;
}


static void upck_f (Buff buff, $t *vptr)
{
   return ;
}


static $t dup_f ($t v)
{
   return (v) ;
}


static int *lowerbd1 (int v_blocksize, int procNr)
{
   int *res = (int *) VMalloc (2 * sizeof (int), "lowerbd1") ;

   res[0] = procNr * v_blocksize ;
   res[1] = 0                    ;

   return (res) ;
}


static int *lowerbd2 (int procNr)
{
   int *res = (int *) VMalloc (2 * sizeof (int), "lowerbd2") ;

   res[0] = procNr ;
   res[1] = 0      ;

   return (res) ;
}


static $t init_f (int size, Index ix)
{
   return (ix[0] == ix[1] ? 1 : ix[1] == size ? 5 * (2 * size - 1) : 2) ;
}


static $t zero (Index ix)
{
   return (0) ;
}


static void destroy_f ($t v)
{
   return ;
}


static $t copy_pivot (array <$t> a, int k, $t v, Index ix)
{
   Bounds  bds   = array_part_bounds (a) ;
   $t      akk                           ;
   char   *FName = "copy_pivot"          ;
   Index   elem_index                    ;

   if (bds->lowerBd[0] <= k && k < bds->upperBd[0])
   {
      elem_index[0] = k ;
      elem_index[1] = k ;

      akk = array_get_elem (a, elem_index) ;

      if (akk == 0.0)
         LogError (EC_ERROR, FName, "pivot element a[%d,%d] = 0.0\n", k, k) ;

      elem_index[1] = ix[1] ;

      return (array_get_elem (a, elem_index) / akk) ;
   }
   else
      return (v) ;
}


static $t eliminate (int k, array <$t> a, array <$t> piv, $t v, Index ix)
{
   Index elem_index1, elem_index2 ;

   if (ix[0] == k || ix[1] < k)
      return (v) ;
   else
   {
      elem_index1[0] = ix[0]  ;
      elem_index1[1] = k      ;
      elem_index2[0] = procId ;
      elem_index2[1] = ix[1]  ;

      return (v - array_get_elem (a, elem_index1) *
	          array_get_elem (piv, elem_index2)) ;
   }
}


static $t normalize (array <$t> a, int refcol, $t v, Index ix)
{
   $t     aii                 ;
   char  *FName = "normalize" ;
   Index  elem_index          ;

   if (ix[1] != refcol)
      return (v) ;

   elem_index[0] = ix[0] ;
   elem_index[1] = ix[0] ;

   aii = array_get_elem (a, elem_index) ;

   if (aii == 0.0)
      LogError (EC_ERROR, FName, "pivot element a[%d,%d] = 0.0\n", ix[0],ix[0]);

   return (v / aii) ;
}


static elemrec make_elemrec ($t v, Index ix)
{
   elemrec res ;

   res.val = v     ;
   res.row = ix[0] ;
   res.col = ix[1] ;

   return (res) ;
}


static elemrec max_abs_in_col (int col, elemrec ref, elemrec new)
{
   elemrec res ;

   if (ref.col != col || ref.row < col)
      if (new.col != col || new.row < col)
      {
         res.val = 0.0 ;
         res.row = -1  ;
         res.col = -1  ;
      }
      else
         res = new ;
   else
      if (new.col != col || new.row < col)
         res = ref ;
      else
         if (rabs (ref.val) > rabs (new.val))
            res = ref ;
         else if (rabs (ref.val) < rabs (new.val))
            res = new ;
         else /* rabs (ref.val) == rabs (new.val) */
            if (ref.row == col)
               res = ref ;
            else if (new.row == col)
               res = new ;
            else /* ref.row != col && new.row != col */
               res = ref ;

   return (res) ;
}


static int switch_rows (int row1, int row2, int crt_row)
{
   int ret ;

   if (row1 == crt_row)
      ret = row2 ;
   else if (row2 == crt_row)
      ret = row1 ;
   else
      ret = crt_row ;

   return (ret) ;
}


int main (int argc, char **argv)
{
   array <float>  a, b, piv                ;
   elemrec        e                        ;
   int            k, size, v_blocksize     ;
   char          *FName = "main"           ;
   Index          totalsize, blocksize, ix ;

   size = atoi (argv[1]) ;

   par_init   () ;
   array_init () ;

   totalsize[0] = size           ;
   totalsize[1] = size + 1       ;

   v_blocksize  = size / netSize ;

   blocksize[0] = v_blocksize    ;
   blocksize[1] = size + 1       ;

   a = array_create (2, totalsize, blocksize, lowerbd1 (v_blocksize),
		     my_rand (size)) ;
   b = array_create (2, totalsize, blocksize, lowerbd1 (v_blocksize),
		     zero) ;

   totalsize[0] = yNetDim ;
   blocksize[0] = 1       ;

   piv = array_create (2, totalsize, blocksize, lowerbd2, zero) ;

   ix[1] = 0 ;

   for (k = 0 ; k < size ; k++)
   {
      e = array_fold (make_elemrec, max_abs_in_col (k), a, pck_f, upck_f,
                      dup_f) ;

      if (e.val == 0.0)
      {
         if (procId == 0)
            printf ("Matrix is singular.\n") ;

         return (0) ;
      }

      if (e.row != k)
         array_permute_rows (a, switch_rows (e.row, k), pck_f, upck_f) ;

      ix[0] = k / v_blocksize ;

      array_map (copy_pivot (a, k), piv, piv)        ;
      array_broadcast_block (piv, ix, pck_f, upck_f) ;
      array_map (eliminate (k, a, piv), a, b)        ;
      array_copy (b, a) ;
   }

   array_map (normalize (a, size), b, a) ;

   array_destroy (piv, destroy_f) ;
   array_destroy (b,   destroy_f) ;
   array_destroy (a,   destroy_f) ;

   array_close () ;
   par_close   () ;

   return (0) ;
}


Back to the Skil home page

Valid HTML 4.01 Strict! Valid CSS!