/* @(#)mo_gaussj.c	17.1.1.1 (ES0-DMD) 01/25/02 17:49:06 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  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 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* @(#)mo_gaussj.c	17.1.1.1 (ESO-DAG) 01/25/02 17:49:06 */

#include <math.h>

#define SWAP(a,b) {float temp=(a);(a)=(b);(b)=temp;}

void MO_GAUSSJ(a,n,b,m)
float **a,**b;
int n,m;
{
	int   *indxc,*indxr,*ipiv;
	int   i,icol,irow,j,k,l,ll;
	float big,dum,pivinv;

	indxc = (int *) osmmget(n * sizeof(int));
	indxr = (int *) osmmget(n * sizeof(int));
	ipiv  = (int *) osmmget(n * sizeof(int));

	for (j = 1; j<= n;j++) ipiv[j]=0;
	for (i = 1; i <= n; i++) {
	   big=0.0;
	   for (j = 1; j <= n; j++)
	      if (ipiv[j] != 1)
	         for (k = 1; k <= n; k++) {
		    if (ipiv[k] == 0) {
		       if (fabs(a[j][k]) >= big) {
		          big  = fabs(a[j][k]);
		          irow = j;
		          icol = k;
		       }
		    } 
              else if (ipiv[k] > 1) 
                 SCETER(2, "*** FATAL: GAUSSJ: Singular Matrix-1");
				}
	   ++(ipiv[icol]);
	   if (irow != icol) {
	      for (l = 1; l <= n; l++) SWAP(a[irow][l],a[icol][l])
	      for (l = 1; l <= m; l++) SWAP(b[irow][l],b[icol][l])
	    }

	   indxr[i]       = irow;
	   indxc[i]       = icol;
	   if (a[icol][icol] == 0.0) 
              SCETER(2, "*** FATAL: GAUSSJ - Singular Matrix-2");
	   pivinv         = 1.0/a[icol][icol];
	   a[icol][icol]  = 1.0;
	   for (l = 1; l <= n; l++) a[icol][l] *= pivinv;
	   for (l = 1; l <= m; l++) b[icol][l] *= pivinv;
	   for (ll = 1; ll <= n; ll++)
	   if (ll != icol) {
	      dum         = a[ll][icol];
	      a[ll][icol] = 0.0;
	      for (l = 1; l <= n; l++) a[ll][l] -= a[icol][l]*dum;
	      for (l = 1; l <= m; l++) b[ll][l] -= b[icol][l]*dum;
	   }
	}

	for (l = n; l >= 1; l--) {
	   if (indxr[l] != indxc[l])
	      for (k = 1; k <= n; k++)
		 SWAP(a[k][indxr[l]],a[k][indxc[l]]);
	}

	osmmfree((char *) ipiv);
	osmmfree((char *) indxr);
	osmmfree((char *) indxc);
}

#undef SWAP