Logo Search packages:      
Sourcecode: felt version File versions  Download package

basic.c

/*
    This file is part of the FElt finite element analysis package.
    Copyright (C) 1993-2000 Jason I. Gobat and Darren C. Atkinson

    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.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/************************************************************************
 * File:    basic.c            
 *    
 * Description:   
 *          
 ************************************************************************/

# include <time.h>
# include <math.h>
# include <stdio.h>
# include "matrix.h"

#undef srand48
#undef drand48
extern void srand48 ();
extern double drand48 ();

# if (defined DOS || defined __CYGWIN32__)
# define srand48 srand
# define drand48 rand
# endif

# define PRINT_TOL 1.0e-12

int ZeroMatrix (a)            /* a = 0                */
   Matrix   a;          /* matrix to fill with zeros  */
{
   unsigned i,j;

   if (IsCompact(a)) 
      for (i = 1 ; i <= Msize(a) ; i++)
         sdata(a, i, 1) = 0.0;
   else {
      for (i = 1 ; i <= Mrows(a) ; i++)
         for (j = 1 ; j <= Mcols(a) ; j++)
            sdata(a, i, j) = 0;
   }

   return 0;
}

int MirrorMatrix (a)
   Matrix   a;          /* matrix to mirror              */
{
   unsigned i, j;

   if (IsCompact (a))
      return M_COMPACT;

   for (i = 2 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j < i ; j++)
         sdata (a, i, j) = mdata (a, j, i);

   return 0;
}

int CopyMatrix (b, a)         /* b = a                   */
   Matrix   a;          /* source matrix           */
   Matrix   b;          /* destination matrix            */
{
   unsigned i, j;

   if (IsCompact(b) && IsFull(a))
       return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mcols(a) != Mcols(b))
      return M_SIZEMISMATCH;

   if (IsCompact(b)) {
      for (i = 1 ; i <= Msize(a) ; i++) {
         b -> data [i][1] = a -> data [i][1];
         b -> diag [i] = a -> diag [i];
      }
   }
   else {
      for (i = 1 ; i <= Mrows(a) ; i++)
         for (j = 1 ; j <= Mcols(a) ; j++)
            sdata(b,i,j) = mdata(a, i, j);
    }

   return 0;
}

int IdentityMatrix (a)        /* a = [I]                 */
   Matrix   a;          /* destination matrix for identity */
{
   unsigned i;
   int            status;

   if (IsCompact(a))
      return M_COMPACT;

   status = ZeroMatrix (a);
   if (status)
      return status;

   if (Mrows(a) < Mcols(a))
      for (i = 1 ; i <= Mrows(a) ; i++)
         sdata(a, i, i) = 1;
   else
      for (i = 1 ; i <= Mcols(a) ; i++)
         sdata(a, i, i) = 1;

   return 0;
}

int RandomMatrix (a, seed)    /* a(i,j) = rand()            */
   Matrix   a;          /* matrix to randomize        */
   int            seed;       /* optional seed        */
{
   unsigned i, j;
   static int     seeded = 0;

   if (IsCompact(a))
      return M_COMPACT;

   if (!seeded) {
      srand48(time(NULL));
      seeded = 1;
   }

   if (seed)
      srand48 (seed);

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(a,i,j) = drand48();

   return 0;
}

int MultiplyMatrices (c, a, b)      /* c = ab               */
   Matrix   c;          /* destination matrix         */
   Matrix   a;          /* source matrix 1            */
   Matrix   b;          /* source matrix 2            */
{
   unsigned i,j,k;

   if (IsCompact(c))
      return M_COMPACT;

   if (c == a || c == b)
      return M_NOOVERWRITE;

   if (Mrows(a) != Mrows(c) || Mcols(b) != Mcols(c) || Mcols(a) != Mrows(b))
      return M_SIZEMISMATCH;
   
   ZeroMatrix (c);

   for (i = 1 ; i <= Mrows(a) ; i++) 
      for (j = 1 ; j <= Mcols(b) ; j++) 
         for (k = 1 ; k <= Mcols(a) ; k++)
            sdata(c,i,j) += mdata(a,i,k) * mdata(b,k,j);

   return 0;
}

int AddMatrices (c, a, b)     /* c = a + b                  */
   Matrix   c;          /* destination matrix         */
   Matrix   a;          /* source matrix 1            */
   Matrix   b;          /* source matrix 2            */
{
   unsigned i,j;

   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;
  
   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = mdata(a, i, j) + mdata(b, i, j);

   return 0;
}

int SubtractMatrices (c, a, b)      /* c = a - b                  */
   Matrix   c;          /* destination matrix         */
   Matrix   a;          /* source matrix 1            */
   Matrix   b;          /* source matrix 2            */
{
   unsigned i,j;

   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;
  
   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = mdata(a, i, j) - mdata(b, i, j);

   return 0;
}

int ModMatrices (c, a, b)     /* c = a % b                  */
   Matrix   c;          /* destination matrix         */
   Matrix   a;          /* source matrix 1            */
   Matrix   b;          /* source matrix 2            */
{
   unsigned i,j;

   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;
  
   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = fmod (mdata(a, i, j), mdata(b, i, j));

   return 0;
}

int Saxpy (c, a, b, alpha)    /* c = b + alpha*a            */    
   Matrix   c;          /* destination vector         */
   Matrix   a;          /* source 1 vector            */
   Matrix   b;          /* source 2 vector            */
   double   alpha;            /* scale factor               */
{
   unsigned i;

   if (IsCompact(c))
      return M_COMPACT;

   if (Mcols(a) != 1 || Mcols(b) != 1 || Mcols(c) != 1)
      return M_SIZEMISMATCH;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      sdata(c,i,1) = alpha*mdata(a, i, 1) + mdata(b, i, 1);

   return 0;
}

int Gaxpy (c, a, b, A)        /* c = b + Aa                 */
   Matrix   c;          /* destination vector         */
   Matrix   a;          /* source 1 vector            */
   Matrix   b;          /* source 2 vector            */
   Matrix   A;          /* source matrix        */
{
   unsigned i,j;

   if (IsCompact(c))
      return M_COMPACT;

   if (a == c)
      return M_NOOVERWRITE;

   if (Mcols(a) != 1 || Mcols(b) != 1 || Mcols(c) != 1)
      return M_SIZEMISMATCH;

   if (Mrows(a) != Mcols(A) || Mrows(b) != Mrows(A) || Mrows(c) != Mrows(b))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(b) ; i++)
      sdata(c,i,1) = mdata(b, i, 1);

   for (i = 1 ; i <= Mcols(A) ; i++)
      for (j = 1 ; j <= Mrows(A) ; j++)
         sdata(c,j,1) += mdata(a, i, 1) * mdata(A, j, i); 

   return 0;
}

int ScaleMatrix(b, a, factor, offset)     /* b(i,j) = factor*a(i,j) + offset  */
   Matrix   b;                /* destination matrix             */
   Matrix   a;                /* source matrix            */
   double   factor;                 /* multiplicative scale factor          */
   double   offset;                 /* additive offset                */
{
   unsigned i,j;

   if (IsCompact(b))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mcols(a) != Mcols(b))
      return M_SIZEMISMATCH;
 
   if (factor == 1) {
      for (i = 1 ; i <= Mrows(a) ; i++)
         for (j = 1 ; j <= Mcols(a) ; j++)
            sdata(b,i,j) = mdata(a, i, j) + offset;
   }
   else if (offset == 0) {
      for (i = 1 ; i <= Mrows(a) ; i++)
         for (j = 1 ; j <= Mcols(a) ; j++)
            sdata(b,i,j) = mdata(a, i, j)*factor;
   }
   else { 
      for (i = 1 ; i <= Mrows(a) ; i++)
         for (j = 1 ; j <= Mcols(a) ; j++)
            sdata(b,i,j) = mdata(a, i, j)*factor + offset;
   }

   return 0;
}

int SqrtMatrix(b, a)          /* b(i,j) = sqrt(a(i,j))        */
   Matrix   b;          /* destination matrix         */
   Matrix   a;          /* source matrix          */
{
   unsigned i,j;

   if (IsCompact(b))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mcols(a) != Mcols(b))
      return M_SIZEMISMATCH;
 
    for (i = 1 ; i <= Mrows(a) ; i++)
       for (j = 1 ; j <= Mcols(a) ; j++)
          sdata(b,i,j) = sqrt (fabs (mdata(a, i, j)));

   return 0;
}

int DotBProduct(x, a, b)      /* x = aTb              */
   double   *x;         /* pointer to result location */
   Matrix   a;          /* source vector (row) 1      */
   Matrix   b;          /* source vector 2            */
{
   unsigned i;
   double   result;

   if (Mrows(a) != 1 || Mcols(b) != 1 || Mcols(a) != Mrows(b))
      return M_SIZEMISMATCH;

   result = 0;
   for (i = 1 ; i <= Mcols(a) ; i++)
      result += mdata(a,1,i)*mdata(b,i,1);

   *x = result;
   return 0;
}

int TransposeMatrix(b, a)     /* b = aT               */
   Matrix   b;          /* destination matrix         */
   Matrix   a;          /* source matrix        */
{
   unsigned i, j;

   if (IsCompact(b))
      return M_COMPACT;

   if (b == a)
      return M_NOOVERWRITE;

   if (Mrows(a) != Mcols(b) || Mcols(a) != Mrows(b))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(b,j,i) = mdata(a, i, j);

   return 0;
}

int CompareEQMatrices (c, a, b)     /* c = (a == b)         */
   Matrix   c;                /* destination matrix   */
   Matrix   a;                /* first RHS matrix     */
   Matrix   b;                /* second RHS matrix    */
{
   unsigned i, j;
  
   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = (mdata(a, i, j) == mdata(b, i, j));

   return 0;
}

int CompareNEQMatrices (c, a, b)    /* c = (a != b)         */
   Matrix   c;                /* destination matrix   */
   Matrix   a;                /* first RHS matrix     */
   Matrix   b;                /* second RHS matrix    */
{
   unsigned i, j;
  
   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = (mdata(a, i, j) != mdata(b, i, j));

   return 0;
}

int CompareGTMatrices (c, a, b)     /* c = (a > b)          */
   Matrix   c;                /* destination matrix   */
   Matrix   a;                /* first RHS matrix     */
   Matrix   b;                /* second RHS matrix    */
{
   unsigned i, j;
  
   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = (mdata(a, i, j) > mdata(b, i, j));

   return 0;
}

int CompareLTMatrices (c, a, b)     /* c = (a < b)          */
   Matrix   c;                /* destination matrix   */
   Matrix   a;                /* first RHS matrix     */
   Matrix   b;                /* second RHS matrix    */
{
   unsigned i, j;
  
   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = (mdata(a, i, j) < mdata(b, i, j));

   return 0;
}

int CompareLTEMatrices (c, a, b)    /* c = (a <= b)         */
   Matrix   c;                /* destination matrix   */
   Matrix   a;                /* first RHS matrix     */
   Matrix   b;                /* second RHS matrix    */
{
   unsigned i, j;
  
   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = (mdata(a, i, j) <= mdata(b, i, j));

   return 0;
}

int CompareGTEMatrices (c, a, b)    /* c = (a >= b)         */
   Matrix   c;                /* destination matrix   */
   Matrix   a;                /* first RHS matrix     */
   Matrix   b;                /* second RHS matrix    */
{
   unsigned i, j;
  
   if (IsCompact(c))
      return M_COMPACT;

   if (Mrows(a) != Mrows(b) || Mrows(b) != Mrows(c))
      return M_SIZEMISMATCH;

   if (Mcols(a) != Mcols(b) || Mcols(b) != Mcols(c))
      return M_SIZEMISMATCH;

   for (i = 1 ; i <= Mrows(a) ; i++)
      for (j = 1 ; j <= Mcols(a) ; j++)
         sdata(c,i,j) = (mdata(a, i, j) >= mdata(b, i, j));

   return 0;
}

int PrintMatrix (m, fp)       /* print matrix m to file fp  */
   Matrix   m;          /* matrix to print            */
   FILE           *fp;        /* file pointer for output    */
{
   double       val;
   unsigned start, end;
   unsigned     i, j;
  
   if (fp == NULL)
      fp = stdout; 

   for (start = 1 ; start <= Mcols(m) ; start += 7) {
      end = start + 6 > Mcols(m) ? Mcols(m) : start + 6;
   
      if (Mcols(m) > 7)
         fprintf (fp, "\nColumns %d through %d\n\n", start, end);
      else
         fprintf (fp, "\n");

      for (j = 1 ; j <= Mrows(m) ; j++) {
         for (i = start ; i <= end ; i++) {
            val = mdata(m,j,i);
            if (fabs(val) < PRINT_TOL)
               val = 0.0;

            fprintf (fp, "%10.4g ", val);
         }
         fprintf (fp, "\n");
      }
   }
     
   fprintf (fp, "\n");
   
   return 0;
}

int PrintMatrixSubsection (m, sr, sc, er, ec, fp) /* print matrix m to fp    */
   Matrix   m;                        /* matrix to print         */
   unsigned sr;                       /* starting row      */
   unsigned sc;                       /* starting column         */
   unsigned er;                       /* ending row              */
   unsigned ec;                       /* ending column           */
   FILE           *fp;                      /* file pointer for output */
{
   double       val;
   unsigned     i, j;

   if (fp == NULL)
      fp = stdout;

   if (sr == 0)
      sr = 1;
   if (sc == 1)
      sc = 1;
   if (er == 0)
      er = Mrows(m);
   if (ec == 0)
      ec = Mcols(m);

   for (i = sr ; i <= er ; i ++) {
      for (j = sc ; j <= ec ; j ++) {
         val = mdata(m, i, j);
         if (fabs(val) < PRINT_TOL)
            val = 0.0;

         fprintf (fp, "%9.4g ", val);
         if (j % 8 == 0)
            fprintf (fp, "\n> ");
      }
      fprintf (fp, "\n");
   }

   fprintf (fp, "\n");
   
   return 0;
}

Generated by  Doxygen 1.6.0   Back to index