/* gjcpeg.c: Gauss-Jordan column pivot for a group of linear system.*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef EPSILON
#define EPSILON 0.000000001
#endif
#define DIM  451
#ifndef GROUP
#define GROUP DIM         /* for computing invert matrix of A */
#endif

int gjcpegd3(int process, double A[DIM][DIM], double xx[DIM][GROUP] )
{


int k,i,j,i0;
double pelement;

for(k=0;k<DIM;k++)
  {
  pelement=fabs(A[k][k]);  i0=k;
  for(i=k+1;i<DIM;i++)
    {
    if( fabs(A[i][k]) > pelement )
       { pelement=fabs(A[i][k]); i0=i; }
    }
  if( i0 != k )
    {
    for(j=k;j<DIM;j++)
      { pelement=A[k][j]; A[k][j]=A[i0][j]; A[i0][j]=pelement; }
    for(j=0;j<GROUP;j++)
      { pelement=xx[k][j]; xx[k][j]=xx[i0][j]; xx[i0][j]=pelement;}
    }
  if( fabs(A[k][k]) < EPSILON ) return(1);
  for(j=k+1;j<DIM; j++ )  A[k][j]= A[k][j]/A[k][k];
  for(j=0;j<GROUP;j++) xx[k][j]=xx[k][j]/A[k][k];
  A[k][k]=1.0;
  for(i=0;i<DIM;i++)
    {
    if( i!=k )
     {
     for(j=0;j<GROUP;j++) xx[i][j]=xx[i][j]-A[i][k]*xx[k][j];
     for(j=k+1;j<DIM;j++) A[i][j]=A[i][j]-A[i][k]*A[k][j];
     A[i][k]=0.0;
     }
    }
  
  }
return(0);
}