gaussjordan.c 3.88 KB
 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 ``````/** * taken from: http://coastal.er.usgs.gov/rvm/toolkit/rvmlibv2.c * unknown license */ /*------------------------------------------------------*/ /* gauss.c */ /* */ /* 20 jan 89 */ /* */ /* Now does full pivoting--row and column swapping. */ /* Requires keeping track of which variable corresponds */ /* to each vector position. */ /* */ /* 18 jan 89 */ /* paul o'neill */ /* */ /* from rob holman's pascal program--geo.pas */ /* */ /* Gauss-Jordan procedure to solve even-determined */ /* square matrix equation Ax = vec,where A is N x N */ /* and vec is N x 1. This is taken roughly from */ /* Menke's book (Geophysical Data Analysis: Discrete */ /* Inverse Theory--1984 p210 ), but performs actual */ /* row switching to simplify the programming. */ /* Partial pivoting is used. */ /* */ /* A[][] is a square matrix, N x N */ /* vec[] is N x 1 of the matrix */ /* nsize is the size of the equation system */ /* */ /* returns 0 if successful */ /* returns -1 if ill-conditioned matrix */ /*------------------------------------------------------*/ `````` Michael Beck committed Jun 06, 2006 35 36 ``````#include "config.h" `````` 37 38 ``````#include #include `````` Michael Beck committed Jun 06, 2006 39 ``````#include "xmalloc.h" `````` 40 `````` `````` Matthias Braun committed Mar 09, 2010 41 42 ``````#include "gaussjordan.h" `````` 43 44 45 46 ``````#define SMALL 0.00001 int firm_gaussjordansolve(double *A, double *vec, int nsize) { `````` Michael Beck committed Apr 26, 2007 47 48 `````` int i, j, row, col, col2, biggest_r = 0, biggest_c = 0, t; double big, temp, sum; `````` 49 50 `````` double *scramvec = XMALLOCN(double, nsize); int *x = XMALLOCN(int, nsize); `````` Michael Beck committed Apr 26, 2007 51 `````` int res = 0; `````` 52 53 `````` #define _A(row,col) A[(row)*nsize + (col)] `````` Michael Beck committed Apr 26, 2007 54 55 56 `````` /* init x[] */ for (i = 0; i < nsize; ++i) x[i] = i; `````` 57 `````` `````` Michael Beck committed Apr 26, 2007 58 `````` /* triangularize A */ `````` Christoph Mallon committed Mar 09, 2011 59 `````` /* ie A has zeros below its diagonal */ `````` Michael Beck committed Apr 26, 2007 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 `````` for (col = 0; col < nsize - 1; ++col) { big = 0; /* find the largest left in LRH box */ for (row = col; row < nsize; ++row) { for (col2 = col; col2 < nsize; ++col2) { temp = fabs(_A(row,col2)); if (temp > big) { biggest_r = row; biggest_c = col2; big = temp; /* largest element left */ } } } if (big < SMALL) { res = -1; `````` 75 76 `````` goto end; } `````` 77 `````` `````` 78 `````` /* swap rows */ `````` Christoph Mallon committed Feb 13, 2010 79 `````` for (i=0;i=0;i--) { `````` 121 `````` sum = 0; `````` Christoph Mallon committed Feb 13, 2010 122 `````` for (j=i+1;j