File: solvers.c

package info (click to toggle)
felt 3.02-4
  • links: PTS
  • area: main
  • in suites: slink
  • size: 16,460 kB
  • ctags: 6,885
  • sloc: ansic: 72,103; fortran: 3,614; yacc: 2,825; lex: 1,172; sh: 311; makefile: 279
file content (96 lines) | stat: -rw-r--r-- 2,411 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/*
    This file is part of the FElt finite element analysis package.
    Copyright (C) 1993-1997 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:	solvers.c
 *	
 * Description:	
 *		
 ************************************************************************/

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

int GaussSeidel(x, A, b)
   Matrix	x;
   Matrix	A;
   Matrix	b;
{
   static int	   maxits = 5000;
   static double   tol = 0.0001;
   int		   iter;
   int		   i, j;
   int		   n;
   double	   sum1, sum2;
   int		   converged;
   double	   new_x;
   double	   diff, base;;

   if (!IsSquare(A))
      return M_NOTSQUARE;

   if (Mrows(b) != Mrows(A) || Mrows(x) != Mrows(A))
      return M_SIZEMISMATCH;

   if (!IsColumnVector(x) || !IsColumnVector(b))
      return M_NOTCOLUMN;

   n = Mrows(A);

   for (i = 1 ; i <= n ; i++) 
      if (mdata(A,i,i) == 0.0)
         return M_SINGULAR;

   converged = 0;

   for (iter = 1 ; iter <= maxits ; iter++) {

      diff = 0.0;
      base = 0.0;
      for (i = 1 ; i <= n ; i++) {

         sum1 = sum2 = 0.0;
         for (j = 1 ; j <= i-1 ; j++)
            sum1 += mdata(A,i,j)*mdata(x,j,1);

         for (j = i+1 ; j <= n ; j++)
            sum2 += mdata(A,i,j)*mdata(x,j,1);

         new_x = (mdata(b,i,1) - sum1 - sum2) / mdata(A,i,i);

         diff += (new_x - mdata(x,i,1))*(new_x - mdata(x,i,1));
         base += new_x*new_x;

         sdata(x, i, 1) = new_x;
      }

      if (base > 0 && diff/base < tol) {
         converged = 1;
         break;
      }
   }

   if (!converged)
      return M_NOTCONVERGED;

   return 0;
}