/* xglu.h                   freeware                   [[Email Removed]]  */

/* --------------------------------- FUNCTION ------------------------------  */
/* Do    :                                                                    */
/*                                                                            */
/* Call  :                                                                    */
/* Debug :                                                                    */
/* -------------------------------------------------------------------------- */
void LU(
pmatrix m,
pmatrix mL,
pmatrix mU
)
{
fraction pivot;
fraction underpivot;
     int rowpivot;
     int rowunder;
     int col;

    mzero(  mL);
   copymF(m,mU);

         col = 0;
for(rowpivot = 0; rowpivot<(mU->rows); rowpivot++)
{
  pivot.numer = *(mU->pblock+rowpivot *mU->cols+col+1);
  pivot.denom = *(mU->pblock+rowpivot *mU->cols+col  );
  if(pivot.denom<0){pivot.numer*= -1;pivot.denom *= -1;}
  if(!pivot.numer)
    {
     printf("\n Pivot error - A pivot is Zero");
     printf("\n Press return to continue");
     getchar();
     exit(1);
    }

  *(mL->pblock+rowpivot *mL->cols+col  )=*(mU->pblock+rowpivot *mU->cols+col  );
  *(mL->pblock+rowpivot *mL->cols+col+1)=*(mU->pblock+rowpivot *mU->cols+col+1);

    mulrowF(mU,rowpivot,pivot);
  frowminiF(mU,rowpivot);

  for(rowunder = rowpivot+1; rowunder<mU->rows; rowunder++)
     {

*(mL->pblock+rowunder *mL->cols+col  )= (*(mU->pblock+rowunder *mU->cols+col ));
*(mL->pblock+rowunder *mL->cols+col+1)=  *(mU->pblock+rowunder *mU->cols+col+1);

                     underpivot.numer =-(*(mU->pblock+rowunder *mU->cols+col ));
                     underpivot.denom =  *(mU->pblock+rowunder *mU->cols+col+1);

        addrowF(mU,underpivot,rowpivot,rowunder);
      frowminiF(mU,                    rowunder);
     }
  col++;
  col++;
}
}