Back to index

extremetuxracer  0.5beta
gaus.h
Go to the documentation of this file.
00001 /* 
00002  * PPRacer 
00003  * Copyright (C) 2003-2004 The PPRacer Team
00004  *
00005  * Copyright (C) 1999-2001 Jasmin F. Patry
00006  * 
00007  * This program is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU General Public License
00009  * as published by the Free Software Foundation; either version 2
00010  * of the License, or (at your option) any later version.
00011  * 
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  * 
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00020  */
00021 
00022 #include <math.h>
00023 
00024 #ifdef EPS
00025 #   undef EPS
00026 #endif
00027 #define EPS 1.0e-10
00028 
00029 
00030 /*
00031   prototype declarations
00032 */
00033 
00034 static unsigned short order(double *matrix, int n, int pivot);
00035 static void elim(double *matrix, int n, int pivot);
00036 static void backsb(double *matrix, int n, double *soln);
00037 
00038 /*
00039   gauss function -- takes a two dimensional array as its arg,
00040   of dimension n x (n+1)
00041 */
00042 
00043 static int gauss(double *matrix, int n, double *soln)
00044 {
00045     int pivot=0;
00046     unsigned short error=0;
00047 
00048 
00049     while ((pivot<(n-1)) && (!error))
00050     {
00051        if(!(error = order(matrix,n,pivot)))
00052        {
00053            elim(matrix,n,pivot);
00054            pivot++;
00055        }
00056     }
00057 
00058     if (error)
00059     {
00060        return 1;
00061     }
00062     else
00063     {
00064        backsb(matrix, n, soln);
00065     }
00066        
00067     return 0;
00068 }
00069 
00070 
00071 static unsigned short order( double *matrix, int n, int pivot)
00072 {
00073     int row, rmax, k;
00074     double temp;
00075     unsigned short error=0;
00076 
00077     rmax = pivot;
00078 
00079     for (row=pivot+1; row<n; row++)
00080     {
00081        if (fabs(*(matrix+row*(n+1)+pivot)) > fabs(*(matrix+rmax*(n+1)+pivot)))
00082            rmax = row;
00083     }
00084 
00085     if (fabs(*(matrix+rmax*(n+1)+pivot)) < EPS )
00086        error = 1;
00087     else if (rmax != pivot)
00088     {
00089        for (k=0; k<(n+1); k++)
00090        {
00091            temp = *(matrix+rmax*(n+1)+k);
00092            *(matrix+rmax*(n+1)+k) = *(matrix+pivot*(n+1)+k);
00093            *(matrix+pivot*(n+1)+k) = temp;
00094        }
00095     }
00096 
00097     return error;
00098 
00099 }
00100 
00101 static void elim(double *matrix, int n, int pivot)
00102 {
00103     int row, col;
00104     double factor;
00105 
00106     for (row = pivot+1; row < n; row++)
00107     {
00108        factor = (*(matrix+row*(n+1)+pivot))/(*(matrix+pivot*(n+1)+pivot));
00109        *(matrix+row*(n+1)+pivot)=0.0;
00110        for (col=pivot+1l; col<n+1; col++)
00111        {
00112            *(matrix+row*(n+1)+col) = *(matrix+row*(n+1)+col) - 
00113               (*(matrix+pivot*(n+1)+col))*factor;
00114        }
00115     }
00116 }
00117 
00118 
00119 static void backsb(double *matrix, int n, double *soln)
00120 {
00121     int row, col;
00122 
00123 
00124     for (row = n-1; row >=0; row--)
00125     {
00126        for (col = n-1; col >= row+1; col--)
00127        {
00128            *(matrix+row*(n+1)+(n)) = *(matrix+row*(n+1)+n) - 
00129               (*(soln+col))*(*(matrix+row*(n+1)+col));
00130        }
00131        *(soln+row) = (*(matrix+row*(n+1)+n))/(*(matrix+row*(n+1)+row));
00132     }
00133 
00134 }