int simeq(double *a,double *b,int n)
/*----------------------------------------------------------------------
 * Solves a linear inhomogeneous equation system with precision "double".
 *
 * Input parameters:
 *    n number of equations and variables.
 *
 * Input and output parameters:
 *    a left-hand-side matrix stored columnwise. It is
 *      destroyed in the computation. The size of matrix a is n by n.
 *    b right-hand-side vector (length n). It is
 *      replaced by the solution vector.
 *
 * Return value:
 *    0 for a normal solution
 *    1 for a singular set of equations (i.e. det(a)=0.)
 *
 * Transformed from FORTRAN to C, enhanced by trivial cases (n=1...3)
 * and by scaling by: Henning Weede.
 * Reference: SSP (scientific subroutine package) by IBM
 *----------------------------------------------------------------------
 */{
   int ij,i1,i2,ia,ib,ic,ix,ixj,ixjx,jjx,jx;
   long i,imax,iqs,it,j,jj,jy,k,ny;
   double biga,save2,tol,det,b0,b1,b2;
   static double *p=(double *)NULL;
   static double *q=(double *)NULL;
   static int nBak=-1;

   /* solve trivial cases analytically */
   switch(n)
      {
      case 1:
         if(1.e-10*fabs(b[0])>1.e10*fabs(a[0]))
            return(1);
         b[0] = b[0]/a[0];
         return(0);
      case 2:
         det = a[0]*a[3]-a[1]*a[2];
         b0 = a[3]*b[0]-a[2]*b[1];
         b1 = a[0]*b[1]-a[1]*b[0];
         if((1.e-10*fabs(b0)>1.e10*fabs(det))||(1.e-10*fabs(b1)>1.e10*fabs(det)))
            return(1);
         b[0] = b0/det;
         b[1] = b1/det;
         return(0);
      case 3:
         det=a[0]*a[4]*a[8]+a[1]*a[5]*a[6]+a[2]*a[3]*a[7]-a[2]*a[4]*a[6]-a[0]*a[5]*a[7]-a[1]*a[3]*a[8];
         b0=(a[4]*a[8]-a[5]*a[7])*b[0]+(a[5]*a[6]-a[3]*a[8])*b[1]+(a[3]*a[7]-a[4]*a[6])*b[2];
         b1=(a[2]*a[7]-a[1]*a[8])*b[0]+(a[0]*a[8]-a[2]*a[6])*b[1]+(a[1]*a[6]-a[0]*a[7])*b[2];
         b2=(a[1]*a[5]-a[2]*a[4])*b[0]+(a[2]*a[3]-a[0]*a[5])*b[1]+(a[0]*a[4]-a[1]*a[3])*b[2];
         if(
            (1.e-10*fabs(b0)>1.e10*fabs(det))||
            (1.e-10*fabs(b1)>1.e10*fabs(det))||
            (1.e-10*fabs(b2)>1.e10*fabs(det))
         ) return(1);
         b[0]=b0/det;
         b[1]=b1/det;
         b[2]=b2/det;
         return(0);
      }

   /* reallocation if necessary */
   if((n!=nBak)||(p==(double *)NULL))
      {
      if(
         ((p=(double *)realloc((void *)p,n*sizeof(double)))==(double *)NULL)
         ||
         ((q=(double *)realloc((void *)q,n*sizeof(double)))==(double *)NULL)
      )  {
         fprintf(stderr,"lsimeq(): Cannot (re)allocate memory.\n");
         exit(1);
         }
      }

   /* pre-scaling */
   for(i1=0; i1<n; ++i1)
      {
      for(i2=0,p[i1]=q[i1]=0.; i2<n; ++i2)
         {
         p[i1]+=a[i1+i2*n]*a[i1+i2*n];
         q[i1]+=a[i2+i1*n]*a[i2+i1*n];
         }
      p[i1]=1./sqrt(sqrt(p[i1]));
      q[i1]=1./sqrt(sqrt(q[i1]));
      }
   for(i1=0; i1<n; ++i1)
      for(i2=0,b[i1]=b[i1]*p[i1]; i2<n; ++i2)
         a[i1+i2*n]=a[i1+i2*n]*p[i1]*q[i2];

   /* forward solution */
   tol=0.0;
   jj=-(long)n;
   for(j=1L; j<=(long)n; ++j)
      {
      jy=j+1L;
      jj += (long)n+1L;
      biga=0.0;
      it=jj-j;
      for(i=j; i<=(long)n; ++i)
         {
         ij=it+i-1L;
         if(((biga<0.0)?-biga:biga)<((a[ij]<0.0)?-a[ij]:a[ij]))
            {
            biga=a[ij];
            imax=i;
            }
         }
      if(((biga<0.0)?-biga:biga)<=tol)
         return(1);
      i1=j+(long)n*(j-2L)-1L;
      it=imax-j;
      for(k=j; k<=(long)n; ++k)
         {
         i1+=n;
         i2=i1+(int)it;
         save2=a[i1];
         a[i1]=a[i2];
         a[i2]=save2;
         a[i1]=a[i1]/biga;
         }
      i1=imax-1L;
      i2=j-1L;
      save2=b[i1];
      b[i1]=b[i2];
      b[i2]=save2/biga;
      if(j==(long)n)
         break;
      iqs=(long)n*(j-1);
      for(ix=(int)jy-1; ix<n; ++ix)
         {
         ixj=(int)iqs+ix;
         it=(int)j-ix-1;
         for(jx=(int)jy-1; jx<n; ++jx)
            {
            ixjx=n*jx+ix;
            jjx=ixjx+(int)it;
            a[ixjx] -= (a[ixj]*a[jjx]);
            }
         b[ix] -= (b[(int)j-1]*a[ixj]);
         }
      }

   /* back solution */
   ny=(long)n-1L;
   it=((long)n)*((long)n);
   for(j=1L; j<=ny; ++j)
      {
      ia=(int)(it-j)-1;
      ib=n-(int)j-1;
      ic=n-1;
      for(k=1L; k<=j; ++k)
         {
         b[ib] -= a[ia]*b[ic];
         ia -= n;
         ic -= 1;
         }
      }

   /* post-scaling */
   for(i1=0; i1<n; ++i1)
      b[i1]=b[i1]*q[i1];

   return(0);

   } /* simeq() */