#define CS	1	/* check subscripts on */

#include "common.h"
#include "corefunc.h"

#include <values.h>

#include "_linalg.h"

#define VAR	FloatHi
#define VARPTR	FloatHiPtr

#define MOIE	1	/* Manually Optimize Inductive Expresions */

/* Pointer to standard linea algebra's functions */
Local(StdLinAlgDataPtr) stdalg;

/* phase dim, numbers of mesh intervals, collocation points, boundary and integral conditions*/
Local(Int2) nphase,ntst,ncol,nbc,nnint;

/* global constants */
Local(Int2) nfree,nrc,nrow,nclm;
Local(Int2) ndecomp,nback;

/*
  Address evaluation constants:
  Loc(A(i,j,k))=Loc(mp->A)+Ak*k+Aj*j+i
*/
Local(int) Ak,Aj,Bk,Cj,A1k,A2k,BBk,CCk,S1k,S2k;
#define Bj nfree
#define A1j nphase
#define A2j nphase
#define BBj nphase
#define CCj nphase
#define S1j nphase
#define S2j nphase
Local(int) ICFj;

/* i - column, j - row, k - block */
#define MATRIXELEMA(mp,i,j,k) (*((mp)->A+Ak*(k)+Aj*(j)+(i)))
#define MATRIXELEMB(mp,i,j,k) (*((mp)->B+Bk*(k)+Bj*(j)+(i)))
#define MATRIXELEMC(mp,i,j,k) (*((mp)->C+(mp)->Ck*(k)+Cj*(j)+(i)))
#define MATRIXELEMD(mp,i,j)   (*((mp)->D+nfree*(j)+(i)))
#define MATRIXELEMA1(mp,i,j,k) (*((mp)->A1+A1k*(k)+A1j*(j)+(i)))
#define MATRIXELEMA2(mp,i,j,k) (*((mp)->A2+A2k*(k)+A2j*(j)+(i)))
#define MATRIXELEMBB(mp,i,j,k) (*((mp)->BB+BBk*(k)+BBj*(j)+(i)))
#define MATRIXELEMCC(mp,i,j,k) (*((mp)->CC+CCk*(k)+CCj*(j)+(i)))
#define MATRIXELEMS1(mp,i,j,k) (*((mp)->S1+S1k*(k)+S1j*(j)+(i)))
#define MATRIXELEMS2(mp,i,j,k) (*((mp)->S2+S2k*(k)+S2j*(j)+(i)))
#define MATRIXELEMFAA(mp,i,j)  (*((mp)->FAA+nphase*(j)+(i)))
#define MATRIXELEMFCC(mp,i)    (*((mp)->FCC+(i)))
#define MATRIXELEMBUF(mp,i)    (*((mp)->BUF+(i)))
#define MATRIXELEMSOL1(mp,i,j)  (*((mp)->SOL1+nphase*(j)+(i)))
#define MATRIXELEMSOL2(mp,i,j)  (*((mp)->SOL2+nphase*(j)+(i)))
#define MATRIXELEMSOL3(mp,i,j)  (*((mp)->SOL3+nphase*(j)+(i)))
#define MATRIXELEME(mp,i,j)     (*((mp)->E+(nbc+nphase+nfree)*(j)+(i)))
#define MATRIXELEMX(mp,i)       (*((mp)->X+(i)))
#define MATRIXELEMP0(mp,i,j)    (*((mp)->P0+nphase*(j)+(i)))
#define MATRIXELEMP1(mp,i,j)    (*((mp)->P1+nphase*(j)+(i)))
#define MATRIXELEMIRF(mp,i,j)   (*((mp)->IRF+nphase*ncol*(j)+(i)))
#define MATRIXELEMICF(mp,i,j)   (*((mp)->ICF+ICFj*(j)+(i)))
#define MATRIXELEMIPR(mp,i,j)   (*((mp)->IPR+nphase*(j)+(i)))
#define MATRIXELEMICF1(mp,i,j)  (*((mp)->ICF1+nphase*(j)+(i)))
#define MATRIXELEMICF2(mp,i,j)  (*((mp)->ICF2+nphase*(j)+(i)))
#define MATRIXELEMIR(mp,i)      (*((mp)->IR+(i)))
#define MATRIXELEMIC(mp,i)      (*((mp)->IC+(i)))


#ifdef Nint
#undef Nint	/* ncbimath.h defines this */
#endif

/* Initialization function */
Entry(void) bvp_Init(StdLinAlgDataPtr Stdalg, Int2 Nphase, 
					      Int2 Ntst, 
					      Int2 Ncol, 
					      Int2 Nbc, 
					      Int2 Nint, 
					      Int2 Nfree) {
  stdalg=Stdalg;

  ndecomp=nback=0;

  nphase=Nphase;
  ntst=Ntst;
  ncol=Ncol;
  nbc=Nbc;
  nnint=Nint;
  nfree=Nfree;
  nrc=nbc+nfree;
  nrow=nphase*ncol;
  nclm=nphase*(ncol+1);

  Ak=(ncol+1)*nphase*nphase*ncol;
  Aj=nphase*(ncol+1);
  Bk=ncol*nphase*nfree;
  Cj=nphase*(ncol+1);
  A1k=nphase*nphase;
  A2k=nphase*nphase;
  BBk=nphase*nfree;
  CCk=nphase*(nphase+nfree);
  S1k=nphase*nphase;
  S2k=nphase*nphase;
  ICFj=nphase*(ncol+1);
}

#ifdef WIN_MSWIN
  static FILE *deb=NULL;
  static void DebOpen(void) {
    if (!deb) deb=fopen("debbvp","wt");
  }
  #define Pr fprintf(deb,
#else
  #define Pr printf(
#endif

#if CS
Boolean Check(Int2 block, Int2 nblock, Int2 sub, Int2 lim, Int2 line) {
  if (block<0 || block>=nblock || sub<0 || sub>=lim) {
    if (Message(MSG_OKC,"Something is out of range:\n"
                   " Block %i not in [0,%i), or\n"
                   " index %i not in [0,%i)\n"
		   "_bvpalg.c, line %i",
		   (int)block,(int)nblock,(int)sub,(int)lim,(int)line) == ANS_CANCEL) {
#ifdef WIN_MSWIN
      div(1,0);
#else
      abort();
#endif
    } else return TRUE;
  }
  return FALSE;
}
#endif



/***********/
/* Vectors */

/* k - block number, j - elem number */

#define VECTNUM(k,j)      	(nphase*ncol*(k)+(j))

#define VECTORELEMFA(vp,k,j) 	(*(vp+VECTNUM((k),(j))))
#define VECTORELEMFC(vp,j) 	(*(vp+VECTNUM(ntst,(j))))

Entry(void) bvp_SetVectorElementFA(Vector vp, Int2 k, Int2 j, VAR e) {
#if CS
  Check(k,ntst,j,nphase*ncol,__LINE__);
#endif
  vp[VECTNUM(k,j)]=e;
}

Entry(VAR) bvp_GetVectorElementFA(Vector vp, Int2 k, Int2 j) {
#if CS
  Check(k,ntst,j,nphase*ncol,__LINE__);
#endif
  return vp[VECTNUM(k,j)];
}

Entry(void) bvp_SetVectorElementFC(Vector vp, Int2 j, VAR e) {
#if CS
  Check(0,1,j,stdalg->GetVectorDim(vp)-nphase*ncol*ntst,__LINE__);
#endif
  vp[VECTNUM(ntst,j)]=e;
}

Entry(VAR) bvp_GetVectorElementFC(Vector vp, Int2 j) {
#if CS
  Check(0,1,j,stdalg->GetVectorDim(vp)-nphase*ncol*ntst,__LINE__);
#endif
  return vp[VECTNUM(ntst,j)];
}

Entry(void) bvp_PrintVector(Vector vp, CharPtr title) {
  int n=stdalg->GetVectorDim(vp),j,k;
  #ifdef WIN_MSWIN
  DebOpen();
  #endif
  Pr "\n%s\n",title);
  Pr "  nphase=%i, ncol=%i, ntst=%i, nfree=%i\n",(int)nphase,(int)ncol,(int)ntst,(int)nfree);
  Pr "  n=%i\n",n);
  Pr "  FA:\n");
  for (k=0; k<ntst; k++) {
    Pr "  Block %i: ",k);
    for (j=0; j<nphase*ncol; j++)
      Pr " %.15lg",bvp_GetVectorElementFA(vp,k,j));
    Pr "\n");
  }
  Pr "  FC: ");
  for (j=0; j<(n-ntst*ncol*nphase); j++)
    Pr " %.15lg",bvp_GetVectorElementFC(vp,j));
  Pr "\n");
}


/************/
/* Matrixes */

/* A matrix */
typedef struct {
  Int2 n,m;		/* dimensions */
  int Ck;
  VAR * A;		/* ntst regular matrices nphase*ncol x nphase*(ncol+1)*/
  VAR * B;		/* ntst parameter matrices nfpar X nphase*ncol */
  VAR * C;		/* ntst condition matrices nphase*(ncol+1) x (nfree+1)*nphase */
  VAR * D;		/* 1 bottom matrix nphase*nfree x (nfree+1)*nphase */
  VAR *A1,*A2,*BB,*CC,*S1,*S2,*FAA,*FCC,*BUF,*SOL1,*SOL2,*SOL3,*E,*X,*P0,*P1;    
                       /* working double arrays  used in Decomp and Solve */
  Int2Ptr IRF,ICF,IPR,ICF1,ICF2,IR,IC;  
                       /* working integer arrays  used in Decomp and Solve */
  /* MUST be the last field */
  VAR align_to_multiple_of_8;
} bvpMatrix, PNTR bvpMatrixPtr, bvpMatrix1[1];


#define MATRIXNUMA  	ntst*(ncol+1)*nphase*nphase*ncol
#define MATRIXNUMB 	ntst*ncol*nphase*nfree
#define MATRIXNUMC(n) 	ntst*((n)-ntst*nphase*ncol)*(ncol+1)*nphase
#define MATRIXNUMD(n) 	nfree*((n)-ntst*nphase*ncol)
#define MATRIXNUMA1 	(ntst+1)*nphase*nphase
#define MATRIXNUMA2 	(ntst+1)*nphase*nphase
#define MATRIXNUMBB 	(ntst+1)*nphase*nfree
#define MATRIXNUMCC 	(ntst+2)*nphase*(nphase+nfree)
#define MATRIXNUMS1 	(ntst+1)*nphase*nphase
#define MATRIXNUMS2 	(ntst+1)*nphase*nphase
#define MATRIXNUMFAA  	ntst*nphase
#define MATRIXNUMFCC  	(nbc+nphase+nfree+2*nphase*nphase+1)  /* nfree->nnint+1*/
#define MATRIXNUMBUF  	(nbc+nphase+nfree)       
#define MATRIXNUMSOL1 	nphase*ntst
#define MATRIXNUMSOL2 	nphase*ntst
#define MATRIXNUMSOL3 	nphase*ntst
#define MATRIXNUME   	(nbc+nphase+nfree)*(nbc+nphase+nfree)
#define MATRIXNUMX      nphase*ncol
#define MATRIXNUMP0     nphase*nphase
#define MATRIXNUMP1     nphase*nphase
#define MATRIXNUMIRF 	(ntst+1)*nphase*ncol
#define MATRIXNUMICF 	(ntst+1)*nphase*(ncol+1)
#define MATRIXNUMIPR 	(ntst+1)*nphase
#define MATRIXNUMICF1 	(ntst+1)*nphase
#define MATRIXNUMICF2 	(ntst+1)*nphase
#define MATRIXNUMIR   	(nbc+nphase+nfree)
#define MATRIXNUMIC   	(nbc+nphase+nfree)


#define MATRSIZE(n,m) (sizeof(bvpMatrix1)+\
                       (MATRIXNUMA+MATRIXNUMB+MATRIXNUMC(n)+MATRIXNUMD(n))*sizeof(VAR) +\
                       (MATRIXNUMA1+MATRIXNUMA2+MATRIXNUMBB+MATRIXNUMCC+MATRIXNUMS1+MATRIXNUMS2+\
                        MATRIXNUMFAA+MATRIXNUMFCC+MATRIXNUMBUF+MATRIXNUMSOL1+MATRIXNUMSOL2+\
                        MATRIXNUMSOL3+MATRIXNUME+MATRIXNUMX+MATRIXNUMP0+MATRIXNUMP1)*sizeof(VAR) +\
                       (MATRIXNUMIRF+MATRIXNUMICF+MATRIXNUMIPR+MATRIXNUMICF1+MATRIXNUMICF2+\
                        MATRIXNUMIR+MATRIXNUMIC)*sizeof(Int2) \
                      )


Entry(void) bvp_PrintMatrix(bvpMatrixPtr mp, CharPtr title, int Level) {
  int n=mp->n,m=mp->m,i,j,k;
  #ifdef WIN_MSWIN
  DebOpen();
  #endif
  Pr "\n%s\n",title);
  Pr "  nphase=%i, ncol=%i, ntst=%i, nfree=%i\n",(int)nphase,(int)ncol,(int)ntst,(int)nfree);
  Pr "  n=%i, m=%i\n",n,m);
if (Level==0) {
  Pr "  A,B\n");
  for (k=0; k<ntst; k++) {
    Pr "  Block %i\n",k);
    for (j=0; j<nphase*ncol; j++) {
      for (i=0; i<nphase*(ncol+1); i++)
        Pr " %.10lg",MATRIXELEMA(mp,i,j,k));
      Pr " | ");
      for (i=0; i<nfree; i++)
        Pr " %.10lg",MATRIXELEMB(mp,i,j,k));
      Pr "\n");
    }
    Pr "\n");
  }
  Pr "  C\n");
  for (k=0; k<ntst; k++) {
    Pr "\n");
    Pr "  Block %i\n",k);
    for (j=0; j<mp->n-ntst*nphase*ncol; j++) {
      for (i=0; i<nphase*(ncol+1); i++)
        Pr " %.10lg",MATRIXELEMC(mp,i,j,k));
      Pr "\n");
    }
  }
  Pr "  D\n");
  for (j=0; j<mp->n-ntst*nphase*ncol; j++) {
   for (i=0; i<nfree; i++)
      Pr " %.10lg",MATRIXELEMD(mp,i,j));
    Pr "\n");
  }
}
if (Level==1) {
  Pr "  A1,A2,BB\n");
  for (k=0; k<ntst; k++) {
    Pr "  Block %i\n",k);
    for (j=0; j<nphase; j++) {
      for (i=0; i<nphase; i++)
        Pr " %.10lg",MATRIXELEMA1(mp,j,i,k));       	/* a transposed matrix is printed !!! */
      Pr " | ");
      for (i=0; i<nphase; i++)
        Pr " %.10lg",MATRIXELEMA2(mp,j,i,k));		/* a transposed matrix is printed !!! */
      Pr " | ");
      for (i=0; i<nfree; i++)
        Pr " %.10lg",MATRIXELEMBB(mp,j,i,k));		/* a transposed matrix is printed !!! */
      Pr "\n");
    }
  }
  for (i=0 ; i<(nphase+nfree)*11; i++) Pr "-");
  Pr "\n");
  Pr "  CC\n");
  for (k=0; k<ntst; k++) {
    Pr "  Block %i\n",k);
    for (j=0; j<nphase+nnint+1; j++) {
      for (i=0; i<nphase; i++)
        Pr " %.10lg",MATRIXELEMCC(mp,i,j,k));		
      Pr "\n");
    }
  }
  Pr "\n");
  Pr "  S1,S2\n");
  for (k=0; k<ntst; k++) {
    Pr "  Block %i\n",k);
    for (j=0; j<nphase; j++) {
      for (i=0; i<nphase; i++)
        Pr " %.10lg",MATRIXELEMS1(mp,j,i,k));       	/* a transposed matrix is printed !!! */
      Pr " | ");
      for (i=0; i<nphase; i++)
        Pr " %.10lg",MATRIXELEMS2(mp,j,i,k));		/* a transposed matrix is printed !!! */
      Pr "\n");
    }
  }
}
if (Level==3) {
  Pr "  FAA\n");
  for (i=0; i<nphase; i++) {
    for (j=0; j<ntst; j++) {
      Pr " %.10lg",MATRIXELEMFAA(mp,i,j));
    }
    Pr "\n");
  }
}
if (Level==4) {
  Pr "  E,XE\n");
  for (i=0; i<(nbc+nphase+nfree); i++) {
    for (j=0; j<(nbc+nphase+nfree); j++) {
      Pr " %.10lg",MATRIXELEME(mp,i,j));     		
    }
    Pr " | ");
    Pr " %.10lg",MATRIXELEMBUF(mp,i));
    Pr "\n");
  }
}
if (Level==5) {
  Pr "  P0\n");
  for (i=0; i<nphase; i++) {
    for (j=0; j<nphase; j++) {
      Pr " %.10lg",MATRIXELEMP0(mp,i,j));
    }
    Pr "\n");
  }
  Pr "\n");
  Pr "  P1\n");
  for (i=0; i<nphase; i++) {
    for (j=0; j<nphase; j++) {
      Pr " %.10lg",MATRIXELEMP1(mp,i,j));
    }
    Pr "\n");
  }
}
if (Level==6) {
  Pr "  FCC\n");
  for (i=0; i<(nbc+nphase+nfree+2*nphase*nphase+1); i++) {
    Pr " %.10lg",MATRIXELEMFCC(mp,i));
  }
  Pr "\n");
}
if (Level==7) {
  Pr "  SOL1\n");
  for (i=0; i<nphase; i++) {
    for (j=0; j<ntst; j++) {
      Pr " %.10lg",MATRIXELEMSOL1(mp,i,j));
    }
    Pr "\n");
  }
  Pr "\n");
  Pr "  SOL2\n");
  for (i=0; i<nphase; i++) {
    for (j=0; j<ntst; j++) {
      Pr " %.10lg",MATRIXELEMSOL2(mp,i,j));
    }
    Pr "\n");
  }
  Pr "\n");
  Pr "  SOL3\n");
  for (i=0; i<nphase; i++) {
    for (j=0; j<ntst; j++) {
      Pr " %.10lg",MATRIXELEMSOL3(mp,i,j));
    }
    Pr "\n");
  }
}


#ifdef WIN_MSWIN
  fflush(deb);
#endif
}

Local(void) MatrixSetPtrs(bvpMatrixPtr mp) {
  CharPtr vp;
  vp=(CharPtr)mp+sizeof(bvpMatrix1);
  mp->A=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMA;
  mp->B=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMB;
  mp->C=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMC(mp->n);
  mp->D=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMD(mp->n);
  mp->A1=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMA1;
  mp->A2=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMA2;
  mp->BB=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMBB;
  mp->CC=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMCC;
  mp->S1=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMS1;
  mp->S2=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMS2;
  mp->FAA=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMFAA;
  mp->FCC=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMFCC;
  mp->BUF=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMBUF;
  mp->SOL1=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMSOL1;
  mp->SOL2=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMSOL2;
  mp->SOL3=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMSOL3;
  mp->E=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUME;
  mp->X=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMX;
  mp->P0=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMP0;
  mp->P1=(VAR *)vp;
  vp+=sizeof(VAR)*MATRIXNUMP1;
  mp->IRF=(Int2Ptr)vp;
  vp+=sizeof(Int2)*MATRIXNUMIRF;
  mp->ICF=(Int2Ptr)vp;
  vp+=sizeof(Int2)*MATRIXNUMICF;
  mp->IPR=(Int2Ptr)vp;
  vp+=sizeof(Int2)*MATRIXNUMIPR;
  mp->ICF1=(Int2Ptr)vp;
  vp+=sizeof(Int2)*MATRIXNUMICF1;
  mp->ICF2=(Int2Ptr)vp;
  vp+=sizeof(Int2)*MATRIXNUMICF2;
  mp->IR=(Int2Ptr)vp;
  vp+=sizeof(Int2)*MATRIXNUMIR;
  mp->IC=(Int2Ptr)vp;
}

/* Creates a Matrix of dimensions n,m */
Entry(bvpMatrixPtr) bvp_CreateMatrix(int n, int m DEBPAR) {
  bvpMatrixPtr mp;
  mp = MemNewDeb(MATRSIZE(n,m) DEBARG);
  mp->n=(Int2)n;
  mp->m=(Int2)m;
  mp->Ck=(mp->n-ntst*ncol*nphase)*(ncol+1)*nphase;
  MatrixSetPtrs(mp);
  return mp;
}

/* Destroys a Matrix and free the memory it occupies */
Entry(bvpMatrixPtr) bvp_DestroyMatrix(bvpMatrixPtr mp) {
  return MemFree(mp);
}

/* Copies a Matrix */
Entry(bvpMatrixPtr) bvp_CopyMatrix(bvpMatrixPtr from, bvpMatrixPtr to DEBPAR) {
  if (!to) to=bvp_CreateMatrix(from->n,from->m DEBARG);
  MemCopy(to,from,MATRSIZE(from->n,from->m));
  MatrixSetPtrs(to);
  return to;
}

/* Clears a Matrix */
Entry(void) bvp_ClearMatrix(bvpMatrixPtr mp) {
  if (mp) {
    MemFill((CharPtr)mp+sizeof(bvpMatrix1),0,MATRSIZE(mp->n,mp->m)-sizeof(bvpMatrix1));
  }
}

/* Returns number of rows */
Entry(Int2) bvp_GetRowNum(bvpMatrixPtr mp) {
  return mp->n;
}

/* Returns number of columns */
Entry(Int2) bvp_GetColNum(bvpMatrixPtr mp) {
  return mp->m;
}

Entry(VARPTR) bvp_GetMatrixA(bvpMatrixPtr mp, Int2 k, Int2Ptr mAj) {
  *mAj=Aj;
  return &MATRIXELEMA(mp,0,0,k);
}

Entry(void) bvp_SetMatrixAElement(bvpMatrixPtr mp, Int2 i, Int2 j, Int2 k, VAR e) {
  MATRIXELEMA(mp,i,j,k)=e;
}

Entry(VAR) bvp_GetMatrixAElement(bvpMatrixPtr mp, Int2 i, Int2 j, Int2 k) {
  return MATRIXELEMA(mp,i,j,k);
}

Entry(void) bvp_SetMatrixBElement(bvpMatrixPtr mp, Int2 i, Int2 j, Int2 k, VAR e) {
  MATRIXELEMB(mp,i,j,k)=e;
}

Entry(VAR) bvp_GetMatrixBElement(bvpMatrixPtr mp, Int2 i, Int2 j, Int2 k) {
  return MATRIXELEMB(mp,i,j,k);
}

Entry(void) bvp_SetMatrixCElement(bvpMatrixPtr mp, Int2 i, Int2 j, Int2 k, VAR e) {
  MATRIXELEMC(mp,i,j,k)=e;
}

Entry(VAR) bvp_GetMatrixCElement(bvpMatrixPtr mp, Int2 i, Int2 j, Int2 k) {
  return MATRIXELEMC(mp,i,j,k);
}

Entry(void) bvp_SetMatrixDElement(bvpMatrixPtr mp, Int2 i, Int2 j, VAR e) {
  MATRIXELEMD(mp,i,j)=e;
}

Entry(VAR) bvp_GetMatrixDElement (bvpMatrixPtr mp, Int2 i, Int2 j) {
  return MATRIXELEMD(mp,i,j);
}

Entry(VAR) bvp_GetMatrixEElement (bvpMatrixPtr mp, Int2 i, Int2 j) {
  return MATRIXELEME(mp,i,j);
}

Entry(VAR) bvp_GetMatrixP0Element (bvpMatrixPtr mp, Int2 i, Int2 j) {
  return MATRIXELEMP0(mp,i,j);
}

Entry(VAR) bvp_GetMatrixP1Element (bvpMatrixPtr mp, Int2 i, Int2 j) {
  return MATRIXELEMP1(mp,i,j);
}

#define A(i,j,k) MATRIXELEMA(mp,i,j,k)
#define B(i,j,k) MATRIXELEMB(mp,i,j,k)
#define C(i,j,k) MATRIXELEMC(mp,i,j,k)
#define D(i,j)   MATRIXELEMD(mp,i,j)
#define FA(k,j)  VECTORELEMFA(vp,k,j)     /* k counts BLOCKS */
#define FC(j)    VECTORELEMFC(vp,j)

#define NDIM    nphase
#define NTST    ntst
#define NCOL    ncol
#define NFPR    nfree


/* Multilpy Matrix by vector */
Entry(Vector) bvp_MultiplyMatrixVector(bvpMatrixPtr mp, Vector vp, Vector vpres DEBPAR) {

#define FASV(k,j)  VECTORELEMFA(vpres,k,j)       /* k counts BLOCKS */
#define FCSV(j)    VECTORELEMFC(vpres,j)


  VAR SUM;
  Int2 I,J,K,nrow;
#if MOIE
  VARPTR pAi,pAj,pAk,pBi,pBj,pBk,pCi,pCk;
  int iCj;
#endif
  nrow=bvp_GetRowNum(mp);
  
  if (!vpres) vpres=stdalg->CreateVector(nrow);

#if MOIE
  pAi=mp->A;
  pBi=mp->B;
#endif
  for (I=0; I<NTST; I++) {
#if MOIE
    pAj=pAi;
    pBj=pBi;
#endif
    for (J=0; J<NDIM*NCOL; J++) { 
#if MOIE
      pAk=pAj;
      for (SUM=0.0,K=0; K<NDIM*NCOL; K++) {
	SUM+=*pAk*FA(I,K); 
	pAk++;
      }
#else
      for (SUM=0.0,K=0; K<NDIM*NCOL; K++) SUM+=A(K,J,I)*FA(I,K); 
#endif
#if MOIE
      pAk=pAj+NDIM*NCOL;
#endif
      if (I<NTST-1) {
#if MOIE
        for (K=0; K<NDIM; K++) {
	  SUM+=*pAk*FA(I+1,K);
	  pAk++;
	}
#else
        for (K=0; K<NDIM; K++) SUM+=A(NDIM*NCOL+K,J,I)*FA(I+1,K);
#endif
      } else {
#if MOIE
        for (K=0; K<NDIM; K++) {
	  SUM+=*pAk*FC(K);
	  pAk++;
	}
#else
        for (K=0; K<NDIM; K++) SUM+=A(NDIM*NCOL+K,J,I)*FC(K);
#endif
      }
#if MOIE
      pBk=pBj;
      for (K=0; K<NFPR; K++) {
	SUM+=*pBk*FC(NDIM+K);
	pBk++;
      }
#else
      for (K=0; K<NFPR; K++) SUM+=B(K,J,I)*FC(NDIM+K);
#endif
     
      FASV(I,J)=SUM;
#if MOIE
      pAj+=Aj;
      pBj+=Bj;
#endif
    }
#if MOIE
    pAi+=Ak;
    pBi+=Bk;
#endif
  }
#if MOIE
  iCj=0;
#endif
  for (J=0; J<(nrow-NTST*NDIM*NCOL); J++) {
#if MOIE
    pCi=mp->C+iCj;
#endif
    for (SUM=0.0,I=0; I<NTST; I++) {
#if MOIE
      pCk=pCi;
      for (K=0; K<NDIM*NCOL; K++) {
	SUM+=*pCk*FA(I,K);
	pCk++;
      }
#else
      for (K=0; K<NDIM*NCOL; K++) SUM+=C(K,J,I)*FA(I,K);
#endif
#if MOIE
      pCk=pCi+NDIM*NCOL;
#endif
      if (I<NTST-1) {
#if MOIE
        for (K=0; K<NDIM; K++) {
	  SUM+=*pCk*FA(I+1,K);
	  pCk++;
	}
#else
        for (K=0; K<NDIM; K++) SUM+=C(NDIM*NCOL+K,J,I)*FA(I+1,K);
#endif
      } else {
#if MOIE
        for (K=0; K<NDIM; K++) {
	  SUM+=*pCk*FC(K);
	  pCk++;
	}
#else
        for (K=0; K<NDIM; K++) SUM+=C(NDIM*NCOL+K,J,I)*FC(K);
#endif
      } 
#if MOIE
      pCi+=mp->Ck;
#endif
    }
    for (K=0; K<NFPR; K++) SUM+=D(K,J)*FC(NDIM+K);
    FCSV(J)=SUM;
#if MOIE
    iCj+=Cj;
#endif
  }
  return vpres;

#undef FASV
#undef FCSV

}

/* Appends a vector to a Matrix as the last row */
Entry(bvpMatrixPtr) bvp_AppendMatrixVector(bvpMatrixPtr mpsrc, Vector vp, bvpMatrixPtr mp DEBPAR) {
  Int2 i,j,k;
  if (!mp) mp=bvp_CreateMatrix(mpsrc->n+1,mpsrc->m DEBARG);
  MemCopy(mp->A,mpsrc->A,MATRIXNUMA*sizeof(VAR));
  MemCopy(mp->B,mpsrc->B,MATRIXNUMB*sizeof(VAR));
  for (k=0; k<ntst; k++)
    MemCopy(mp->C+k*(MATRIXNUMC(mpsrc->n+1)/ntst),
            mpsrc->C+k*(MATRIXNUMC(mpsrc->n)/ntst),
            (MATRIXNUMC(mpsrc->n)/ntst)*sizeof(VAR));
  MemCopy(mp->D,mpsrc->D,MATRIXNUMD(mpsrc->n)*sizeof(VAR));
  j=mpsrc->n-ntst*nphase*ncol;
  k=0;
  for (i=0; i<ncol*nphase; i++) C(i,j,k)=FA(0,i);
  for (k=1; k<ntst-1; k++) {
    for (i=0; i<nphase; i++)
      C(i,j,k)=C(nphase*ncol+i,j,k-1)=0.5*FA(k,i);
    for (i=nphase; i<ncol*nphase; i++)
      C(i,j,k)=FA(k,i);
    for (i=ncol*nphase; i<(ncol+1)*nphase; i++)
      C(i,j,k)=C(i-nphase*ncol,j,k+1)=0.5*FA(k+1,i-nphase*ncol);
  }
  for (i=0; i<nphase; i++)
    C(i,j,ntst-1)=C(nphase*ncol+i,j,ntst-2)=0.5*FA(ntst-1,i);
  for (i=nphase; i<ncol*nphase; i++)
    C(i,j,ntst-1)=FA(ntst-1,i);
  for (i=ncol*nphase; i<(ncol+1)*nphase; i++)
    C(i,j,ntst-1)=FC(i-nphase*ncol);
  for (i=0; i<nfree; i++)
    D(i,j)=FC(i+nphase);
  return mp;
}

#undef MATRSIZE


/******************************/
/* BVP decomposer and solver  */

#define IRF(i,j)	MATRIXELEMIRF(mp,i,j)
#define ICF(i,j)	MATRIXELEMICF(mp,i,j)

#define NOV	nphase
#define NA	ntst
#define NRA	nrow
#define NCA	nclm
#define NCB	nfree
#define NBC	nbc
#define NRC  	nrc

Entry(Int2) bvp_conpar(bvpMatrixPtr mp) {
/**** CONDENSATION OF PARAMETERS (ELIMINATION OF LOCAL VARIABLES) ******/

int NEX,M1,M2,IR1,IRP,ICP1,IPIV,JPIV,ITMP;
int I,J,IC,K1,K2,L,IR;
VAR ZERO,PIV,TPIV,RM;

#if MOIE
  VARPTR pAi,pAj,pAl,pA1l,pBi,pBl,pB1l,pCi,pCir,pCj,pCl,pAw;
  Int2Ptr pICF_i,pICFk2i,pICFli;
  VAR wA;
  int icf,irf_ir,irf_irp;
#endif
  
  NEX = NCA-2*NOV;
  if (NEX == 0) return 0; 

  M1    = NOV+1;
  M2    = NOV+NEX;

/* INITIALIZATION  */

  for (I=0; I<NA; I++) {
    for (J=0; J<NRA; J++) IRF(J,I)=J;
    for (J=0; J<NCA; J++) ICF(J,I)=J;
  }

  ZERO = 0.0;
#if MOIE
  pAi=mp->A;
  pBi=mp->B;
  pCi=mp->C;
  pICF_i=mp->ICF;
#endif
  for (I=0; I<NA; I++) {
    for (IC=M1-1; IC<M2; IC++) {
      IR1=IC-NOV+1;
      IRP=IR1-1;
      ICP1=IC+1;

/* SEARCH FOR PIVOT (COMPLETE PIVOTING) */

      PIV = ZERO;
      for (K1=IRP; K1<NRA; K1++) {
#if MOIE
	pAj=pAi+Aj*IRF(K1,I);
	pICFk2i=pICF_i+IC;
#endif
        for (K2=IC; K2<M2; K2++) { 
#if MOIE
          TPIV = *(pAj+*pICFk2i);
	  pICFk2i++;
#else
          TPIV = A(ICF(K2,I),IRF(K1,I),I);
#endif
          if (TPIV < ZERO) TPIV=-TPIV;
          if (PIV < TPIV) {
            PIV = TPIV;
            IPIV = K1;
            JPIV = K2;
          }
        }
      }

/* INDICES ARE CHANGED */
       
      ITMP        = ICF(IC,I);
      ICF(IC,I)   = ICF(JPIV,I);
      ICF(JPIV,I) = ITMP;

      ITMP        = IRF(IRP,I);
      IRF(IRP,I)  = IRF(IPIV,I);
      IRF(IPIV,I) = ITMP;

/* ELIMINATION */

#if MOIE
      wA=A(ICF(IC,I),IRF(IRP,I),I);
      pAj=pAi+ICF(IC,I);
      irf_irp=IRF(IRP,I);
      pAw=pAi+Aj*irf_irp;
#endif
      for (IR=IR1; IR<NRA; IR++) {
#if MOIE
	irf_ir=IRF(IR,I);
        RM=*(pAj+Aj*irf_ir)/wA;
        *(pAj+Aj*irf_ir)=RM;
#else
        RM=A(ICF(IC,I),IRF(IR,I),I)/A(ICF(IC,I),IRF(IRP,I),I);
        A(ICF(IC,I),IRF(IR,I),I)=RM;
#endif
        if (RM != 0.0) {
#if MOIE
	  pAl=pAi+Aj*irf_ir;
	  pA1l=pAw;
#endif
          for (L=0; L<NOV; L++) {
#if MOIE
	    *pAl-=RM**pA1l;
	    pAl++;
	    pA1l++;
#else
	    A(L,IRF(IR,I),I)-=RM*A(L,IRF(IRP,I),I);
#endif
	  }	 
#if MOIE
	  pAl=pAi+Aj*irf_ir;
	  pA1l=pAw;
	  pICFli=pICF_i+ICP1;
#endif
          for (L=ICP1; L<NCA; L++) {
#if MOIE
	    icf=*pICFli++;
	    *(pAl+icf)-=RM**(pA1l+icf);
#else
            A(ICF(L,I),IRF(IR,I),I)-=RM*A(ICF(L,I),IRF(IRP,I),I);
#endif
	  }
#if MOIE
	  pBl=pBi+Bj*irf_ir;
	  pB1l=pBi+Bj*irf_irp;
#endif
          for (L=0; L<NCB; L++) {
#if MOIE
	    *pBl-=RM**pB1l;
	    pBl++;
	    pB1l++;
#else
	    B(L,IRF(IR,I),I)-=RM*B(L,IRF(IRP,I),I);
#endif
	  }
        }
      }
#if MOIE
      pCj=pCi+Cj*NBC;
      pCir=pCj+ICF(IC,I);
      pAj=pAw;
#endif
      for (IR=NBC; IR<NRC; IR++) {
#if MOIE
        RM=*pCir/wA;
        *pCir=RM;
#else
        RM=C(ICF(IC,I),IR,I)/A(ICF(IC,I),IRF(IRP,I),I);
        C(ICF(IC,I),IR,I)=RM;
#endif
        if (RM != 0.0) {
#if MOIE
	  pCl=pCj;
	  pAl=pAj;
#endif
          for (L=0; L<NOV; L++) {
#if MOIE
            *pCl-=RM**pAl;
	    pCl++;
	    pAl++;
#else
            C(L,IR,I)-=RM*A(L,IRF(IRP,I),I);
#endif
	  }
#if MOIE
	  pCl=pCj;
	  pAl=pAj;
	  pICFli=pICF_i+ICP1;
#endif
          for (L=ICP1; L<NCA; L++) {
#if MOIE
	    icf=*pICFli++;
            *(pCl+icf)-=RM**(pAl+icf);
#else
            C(ICF(L,I),IR,I)-=RM*A(ICF(L,I),IRF(IRP,I),I);
#endif
	  }
#if MOIE
	  pB1l=pBi+Bj*irf_irp;
#endif
          for (L=0; L<NCB; L++) {
#if MOIE
            D(L,IR)-=RM**pB1l;
	    pB1l++;
#else
            D(L,IR)-=RM*B(L,IRF(IRP,I),I);
#endif
	  }
        }	/* if (RM != 0.0) */
#if MOIE
	pCj+=Cj;
	pCir+=Cj;
#endif
      }		/* for (IR=NBC; IR<NRC; IR++) */
    }		/* for (IC=M1-1; IC<M2; IC++) */
#if MOIE
    pAi+=Ak;
    pBi+=Bk;
    pCi+=mp->Ck;
    pICF_i+=ICFj;
#endif
  }		/* for (I=0; I<NA; I++) */
  return 0;
}

#define A1(i,j,k)	MATRIXELEMA1(mp,i,j,k)
#define A2(i,j,k)	MATRIXELEMA2(mp,i,j,k)
#define BB(i,j,k)       MATRIXELEMBB(mp,i,j,k)
#define CC(i,j,k)	MATRIXELEMCC(mp,i,j,k)

Entry(Int2) bvp_copycp(bvpMatrixPtr mp) {
/****** COPIES THE CONDENSED SYSTEM GENERATED BY CONPAR INTO THE WORKSPACE *******/
Int2 NAP1,IC1;
Int2 I,IR,IC,IRFIR;

#if MOIE
  VARPTR pAi,pAic,pAic1;
  VARPTR pA1i,pA1ir,pA1ic;
  VARPTR pA2i,pA2ir,pA2ic;
  VARPTR pBi,pBic;
  VARPTR pBBi,pBBir,pBBic;
  VARPTR pCCi,pCCir,pCCic;
  VARPTR pCi,pCir,pCic;
  VARPTR pC1ir,pC1ic;
#endif

  NAP1=NA+1;
#if MOIE
  pAi=mp->A;
  pA1i=mp->A1;
  pA2i=mp->A2;
  pBi=mp->B;
  pBBi=mp->BB;
#endif
  for (I=0; I<NA; I++) {
#if MOIE
    pA1ir=pA1i;
    pA2ir=pA2i;
    pBBir=pBBi;
#endif
    for (IR=0; IR<NOV; IR++) {
      IRFIR=IRF(NRA-NOV+IR,I);
#if MOIE
      pAic=pAi+Aj*IRFIR;
      pAic1=pAic+NCA-NOV+0;	/* IC1=NCA-NOV+IC, see below */
      pA1ic=pA1ir;
      pA2ic=pA2ir;
#endif
      for (IC=0; IC<NOV; IC++) {
#if MOIE
        *pA1ic=*pAic;
        *pA2ic=*pAic1;
	pAic++;
	pAic1++;
	pA1ic+=A1j;
	pA2ic+=A2j;
#else
        IC1=NCA-NOV+IC;
        A1(IR,IC,I)=A(IC,IRFIR,I);
        A2(IR,IC,I)=A(IC1,IRFIR,I);
#endif
      } 
#if MOIE
      pBic=pBi+Bj*IRFIR;
      pBBic=pBBir;
      for (IC=0; IC<NCB; IC++) {
	*pBBic=*pBic;
	pBic++;
	pBBic+=BBj;
      }
#else
      for (IC=0; IC<NCB; IC++) BB(IR,IC,I)=B(IC,IRFIR,I);
#endif
#if MOIE
    pA1ir++;
    pA2ir++;
    pBBir++;
#endif
    }
#if MOIE
  pAi+=Ak;
  pA1i+=A1k;
  pA2i+=A2k;
  pBi+=Bk;
  pBBi+=BBk;
#endif
  }
#if MOIE
  pCCi=mp->CC;
  pCi=mp->C;
#endif
  for (I=0; I<NAP1; I++) {
#if MOIE
    pCCir=pCCi;
    pCir=pCi;
    pC1ir=pCi-mp->Ck;	/* C(...,...,I-1) */
#endif
    for (IR=0; IR<NRC; IR++) {
#if MOIE
      pCCic=pCCir;
      pCic=pCir;
      pC1ic=pC1ir+NRA;	/* C(NRA+IC...) */
#endif
      for (IC=0; IC<NOV; IC++) {
        if (I == 0) 
#if MOIE
          *pCCic=*pCic;
#else
          CC(IC,IR,I)=C(IC,IR,I);
#endif
        else {
          if (I == NA) 
#if MOIE
            *pCCic=*pC1ic;
#else
            CC(IC,IR,I)=C(NRA+IC,IR,I-1);
#endif
          else 
#if MOIE
            *pCCic=*pCic+*pC1ic;
#else
            CC(IC,IR,I)=C(IC,IR,I)+C(NRA+IC,IR,I-1);
#endif
        }
#if MOIE
	pCCic++;
	pCic++;
	pC1ic++;
#endif
      }
#if MOIE
      pCCir+=CCj;
      pCir+=Cj;
      pC1ir+=Cj;
#endif
    }
#if MOIE
  pCCi+=CCk;
  pCi+=mp->Ck;
#endif
  }
  return 0;
}

#define S1(i,j,k)	MATRIXELEMS1(mp,i,j,k)
#define S2(i,j,k)	MATRIXELEMS2(mp,i,j,k)
#define E(i,j)		MATRIXELEME(mp,i,j)
#define P0(i,j)  	MATRIXELEMP0(mp,i,j)
#define P1(i,j)  	MATRIXELEMP1(mp,i,j)
#define ICF1(i,j)	MATRIXELEMICF1(mp,i,j)
#define ICF2(i,j)	MATRIXELEMICF2(mp,i,j)
#define IPR(i,j)	MATRIXELEMIPR(mp,i,j)             

Entry(Int2) bvp_reduce(bvpMatrixPtr mp) {

Int2 NAM1,IPIV1,JPIV1,IPIV2,JPIV2,ITMP;
Int2 NOVPI,NOVPJ,NOVPJ2;
Int2 I,I1,I2,I3,J,K1,K2,IR,IC,ICP1,L;
VAR ZERO,PIV1,PIV2,TPIV,TMP,RM;

#if MOIE
  VARPTR pS1i,pS1k2,pS1k1;
  VARPTR pS2i,pS2k2,pS2k1;
  VARPTR pA2__i1,pA2k1_i1,pA2ic_i1,pA2ipiv_i1,pA2ir_i1;
  VARPTR pA2__i2,pA2ipivli2;
  VARPTR pA1__i2,pA1k1_i2,pA1ipiv_i2;
  VARPTR pS1__i1,pS1ic_i1,pS1icli1,pS1ipivli1;
  VARPTR pS2__i1,pS2ic_i1,pS2icli1,pS2ipivli1;
  VAR wA2;
  int iw,jw;
#endif

  ZERO    = 0.0;
  NAM1    = NA-1;

/*  INITIALIZATION */

#if MOIE
  pS1i=mp->S1;
  pS2i=mp->S2;
#endif
  for (I=0; I<NA; I++) {
#if MOIE
    pS1k2=pS1i;
    pS2k2=pS2i;
#endif
    for (K1=0; K1<NOV; K1++) {
      ICF1(K1,I)     = K1;
      ICF2(K1,I)     = K1;
      IPR(K1,I)      = K1;
#if MOIE
      pS1k1=pS1k2;
      pS2k1=pS2k2;
#endif
      for (K2=0; K2<NOV; K2++) {
#if MOIE
        *pS1k1 = 0.0;
        *pS2k1 = 0.0;
	pS1k1++;
	pS2k1++;
#else
        S2(K1,K2,I) = 0.0;
        S1(K1,K2,I) = 0.0;
#endif
      }
#if MOIE
    pS1k2+=S1j;
    pS2k2+=S2j;
#endif
    }
#if MOIE
  pS1i+=S1k;
  pS2i+=S2k;
#endif
  }
  for (IR=0; IR<NOV; IR++) { 
    for (IC=0; IC<NOV; IC++) S1(IR,IC,0)=A1(IR,IC,0);
  }

/*    REDUCE PROCESS   */

#if MOIE
  pA2__i1=mp->A2;
  pA2__i2=mp->A2+A2k;	/* I2=I1+1 */
  pA1__i2=mp->A1+A2k;	/* ---'--- */
  pS1__i1=mp->S1;
  pS2__i1=mp->S2;
#endif
  for (I1=0; I1<NAM1; I1++) {
    I2=I1+1;
    I3=I2+1;
#if MOIE
    pS1ic_i1=pS1__i1;
    pS2ic_i1=pS2__i1;
    pA2ic_i1=pA2__i1;
#endif
    for (IC=0; IC<NOV; IC++) {
      ICP1=IC+1;

/*    PIVOTING (COMPLETE PIVOTING)  */
/*    ROWS ARE SWAPPED PHYSICALLY, COLUMNS ONLY SWAP INDEXES  */

      PIV1 = ZERO;
#if MOIE
      pA2k1_i1=pA2ic_i1;
#endif
      for (K1=IC; K1<NOV; K1++) {
        for (K2=IC; K2<NOV; K2++) {
#if MOIE
          TPIV      = *(pA2k1_i1+A2j*ICF2(K2,I1));
#else
          TPIV      = A2(K1,ICF2(K2,I1),I1);
#endif
          if (TPIV < ZERO) TPIV=-TPIV;
          if (PIV1 < TPIV) {
            PIV1   = TPIV;
            IPIV1  = K1;
            JPIV1  = K2;
          }
        }
#if MOIE
	pA2k1_i1++;
#endif
      }
      PIV2 = ZERO;
#if MOIE
      pA1k1_i2=pA1__i2;
#endif
      for (K1=0; K1<NOV; K1++) {
        for (K2=IC; K2<NOV; K2++) {
#if MOIE
          TPIV      = *(pA1k1_i2+A1j*ICF1(K2,I2));
#else
          TPIV      = A1(K1,ICF1(K2,I2),I2);
#endif
          if (TPIV < ZERO) TPIV=-TPIV;
          if (PIV2 < TPIV) {
            PIV2   = TPIV;
            IPIV2  = K1;
            JPIV2  = K2;
          }
        }
#if MOIE
	pA1k1_i2++;
#endif
      }
      if (PIV1 >= PIV2) {
        IPR(IC,I1)        = IPIV1;
        ITMP              = ICF2(IC,I1);
        ICF2(IC,I1)       = ICF2(JPIV1,I1);
        ICF2(JPIV1,I1)    = ITMP;
        ITMP              = ICF1(IC,I2);
        ICF1(IC,I2)       = ICF1(JPIV1,I2);
        ICF1(JPIV1,I2)    = ITMP;

/*  SWAPPING  */
#if MOIE
	pS1icli1=pS1ic_i1;
	pS2icli1=pS2ic_i1;
	pS1ipivli1=pS1__i1+IPIV1;
	pS2ipivli1=pS2__i1+IPIV1;
	pA2ipiv_i1=pA2__i1+IPIV1;
#endif
        for (L=0; L<NOV; L++) {
#if MOIE
          TMP            = *pS1icli1;
          *pS1icli1      = *pS1ipivli1;
          *pS1ipivli1   = TMP;
#else
          TMP            = S1(IC,L,I1);
          S1(IC,L,I1)    = S1(IPIV1,L,I1);
          S1(IPIV1,L,I1) = TMP;
#endif
          if (L >= IC) {
#if MOIE
	    iw=A2j*ICF2(L,I1);
            TMP=*(pA2ic_i1+iw);
            *(pA2ic_i1+iw)=*(pA2ipiv_i1+iw);
            *(pA2ipiv_i1+iw)= TMP;
#else
            TMP=A2(IC,ICF2(L,I1),I1);
            A2(IC,ICF2(L,I1),I1)=A2(IPIV1,ICF2(L,I1),I1);
            A2(IPIV1,ICF2(L,I1),I1)= TMP;
#endif
          }
#if MOIE
          TMP            = *pS2icli1;
          *pS2icli1      = *pS2ipivli1;
          *pS2ipivli1 = TMP;
#else
          TMP            = S2(IC,L,I1);
          S2(IC,L,I1)    = S2(IPIV1,L,I1);
          S2(IPIV1,L,I1) = TMP;
#endif
#if MOIE
	  pS1icli1+=S1j;
	  pS2icli1+=S2j;
	  pS1ipivli1+=S1j;
	  pS2ipivli1+=S2j;
#endif
        }
	/* Don't optimize - NCB is small (~2) */
        for (L=0; L<NCB; L++) {
          TMP            = BB(IC,L,I1);
          BB(IC,L,I1)    = BB(IPIV1,L,I1);
          BB(IPIV1,L,I1) = TMP;
        }
      } else {
        IPR(IC,I1)        = NOV+IPIV2;
        ITMP              = ICF2(IC,I1);
        ICF2(IC,I1)       = ICF2(JPIV2,I1);
        ICF2(JPIV2,I1)    = ITMP;
        ITMP              = ICF1(IC,I2);
        ICF1(IC,I2)       = ICF1(JPIV2,I2);
        ICF1(JPIV2,I2)    = ITMP;               

/*  SWAPPING  */

#if MOIE
	pA1ipiv_i2=pA1__i2+IPIV2;
	pS1icli1=pS1ic_i1;
	pS2icli1=pS2ic_i1;
	pA2ipivli2=pA2__i2+IPIV2;
#endif
        for (L=0; L<NOV; L++) {
          if (L >= IC) {
#if MOIE
	    iw=A2j*ICF2(L,I1);
	    jw=A1j*ICF2(L,I1);
            TMP  = *(pA2ic_i1+iw);
            *(pA2ic_i1+iw)=*(pA1ipiv_i2+jw);
            *(pA1ipiv_i2+jw) = TMP;
#else
            TMP  = A2(IC,ICF2(L,I1),I1);
            A2(IC,ICF2(L,I1),I1)=A1(IPIV2,ICF2(L,I1),I2);
            A1(IPIV2,ICF2(L,I1),I2) = TMP;
#endif
          }
#if MOIE
          TMP            = *pS2icli1;
          *pS2icli1      = *pA2ipivli2;
          *pA2ipivli2    = TMP;
	  pS2icli1+=S2j;
	  pA2ipivli2+=A2j;
#else
          TMP            = S2(IC,L,I1);
          S2(IC,L,I1)    = A2(IPIV2,L,I2);
          A2(IPIV2,L,I2) = TMP;
#endif
#if MOIE
          TMP            = *pS1icli1;
          *pS1icli1      = S1(IPIV2,L,I2);
          S1(IPIV2,L,I2) = TMP;                 
	  pS1icli1+=S1j;
#else
          TMP            = S1(IC,L,I1);
          S1(IC,L,I1)    = S1(IPIV2,L,I2);
          S1(IPIV2,L,I2) = TMP;                 
#endif
        }
	/* Don't optimize - NCB is small (~2) */
        for (L=0; L<NCB; L++) {
          TMP            = BB(IC,L,I1);
          BB(IC,L,I1)    = BB(IPIV2,L,I2);
          BB(IPIV2,L,I2) = TMP;
        }
      }
/*  END OF PIVOTING  */

/*  ELIMINATION   */
#if MOIE
      wA2=*(pA2ic_i1+A2j*ICF2(IC,I1));
      pA2ir_i1=pA2__i1+A2j*ICF2(IC,I1)+ICP1;
#endif
      for (IR=ICP1; IR<NOV; IR++) {
#if MOIE
        RM = *pA2ir_i1/wA2;
        *pA2ir_i1   = RM;
#else
        RM = A2(IR,ICF2(IC,I1),I1)/A2(IC,ICF2(IC,I1),I1);
        A2(IR,ICF2(IC,I1),I1)   = RM;
#endif
        if (RM != 0.0) {
          for (L=ICP1; L<NOV; L++) {
            A2(IR,ICF2(L,I1),I1) -=RM*A2(IC,ICF2(L,I1),I1);
          }
#if MOIE
	  pS1icli1=pS1ic_i1;
	  pS2icli1=pS2ic_i1;
#endif
          for (L=0; L<NOV; L++) {
#if MOIE
            S1(IR,L,I1) -= RM**pS1icli1;
            S2(IR,L,I1) -= RM**pS2icli1;
	    pS1icli1+=S1j;
	    pS2icli1+=S2j;
#else
            S1(IR,L,I1) -= RM*S1(IC,L,I1);
            S2(IR,L,I1) -= RM*S2(IC,L,I1);
#endif
          }
          for (L=0; L<NCB; L++) {
            BB(IR,L,I1) = BB(IR,L,I1)-RM*BB(IC,L,I1);
          }
        }
#if MOIE
	pA2ir_i1++;
#endif
      }
      for (IR=0; IR<NOV; IR++) {
#if MOIE
        RM = A1(IR,ICF1(IC,I2),I2)/wA2;
        A1(IR,ICF1(IC,I2),I2)   = RM;
#else
        RM = A1(IR,ICF1(IC,I2),I2)/A2(IC,ICF2(IC,I1),I1);
        A1(IR,ICF1(IC,I2),I2)   = RM;
#endif
        if (RM != 0.0) {
          for (L=ICP1; L < NOV; L++) {
            A1(IR,ICF1(L,I2),I2) -= RM*A2(IC,ICF2(L,I1),I1);
          }
          for (L=0; L<NOV; L++) {
            S1(IR,L,I2) -= RM*S1(IC,L,I1);
            A2(IR,L,I2) -= RM*S2(IC,L,I1);
          }
          for (L=0; L<NCB; L++) BB(IR,L,I2) -= RM*BB(IC,L,I1);
        }
      }
      for (IR=NBC; IR<NRC; IR++) {
#if MOIE
        RM = CC(ICF2(IC,I1),IR,I2)/wA2;
        CC(ICF2(IC,I1),IR,I2)   = RM;                  
#else
        RM = CC(ICF2(IC,I1),IR,I2)/A2(IC,ICF2(IC,I1),I1);
        CC(ICF2(IC,I1),IR,I2)   = RM;                  
#endif
        if (RM != 0.0) {
          for (L=ICP1; L<NOV; L++)
            CC(ICF2(L,I1),IR,I2) -= RM*A2(IC,ICF2(L,I1),I1);
          for (L=0; L<NOV; L++) {
            CC(L,IR,0)  -= RM*S1(IC,L,I1);
            CC(L,IR,I3) -= RM*S2(IC,L,I1);
          }  
          for (L=0; L<NCB; L++) D(L,IR) -= RM*BB(IC,L,I1);
        }
      }
#if MOIE
      pS1ic_i1++;
      pS2ic_i1++;
      pA2ic_i1++;
#endif
    }           
#if MOIE
    pA2__i1+=A2k;
    pA2__i2+=A2k;
    pA1__i2+=A1k;
    pS1__i1+=S1k;
    pS2__i1+=S2k;
#endif
  }            

/*   COPYING TO P0,P1 (here, not in dimrge) and E (also in dimrge) */
  for (I=0; I<NOV; I++) {
    for (J=0; J<NOV; J++) {
      P0(I,J)    = S1(I,J,NA-1);
      P1(I,J)    = A2(I,J,NA-1);
    }
  }
  for (I=0; I<NOV; I++) {
    for (J=0; J<NOV; J++) {
      NOVPJ      = NOV+J;
      E(I,J)     = S1(I,J,NA-1);
      E(I,NOVPJ) = A2(I,J,NA-1);
    }
    for (J=0; J<NCB; J++) {
      NOVPJ2      = 2*NOV+J;
      E(I,NOVPJ2) = BB(I,J,NA-1);
    }
  }

  for (I=0; I<NRC; I++) {     
    NOVPI=NOV+I;
    for (J=0; J<NOV; J++) {
      NOVPJ          = NOV+J;
      E(NOVPI,J)     = CC(J,I,0);
      E(NOVPI,NOVPJ) = CC(J,I,NA);
    }
    for (J=0; J<NCB; J++) {
      NOVPJ2          = 2*NOV+J;
      E(NOVPI,NOVPJ2) = D(J,I);
    }
  }
/* bvp_PrintMatrix(mp,"P-matrices:",5); */

  for (I=0; I<NOV; I++) ICF2(I,NA)=I;
  return 0;
}

Entry(Int2) bvp_conrhs(bvpMatrixPtr mp, Vector vp) {
/******* CONDENSATION OF RIGHT HAND SIDE  *******/

Int2 NEX,NBCP1,M1,M2;
Int2 I,IC,IR1,IRP,IRFIRP,ICFIC,IRFIR,IR;

  NEX=NCA-2*NOV;
  if (NEX == 0) return 0;

  NBCP1 = NBC+1;
  M1    = NOV+1;
  M2    = NOV+NEX;

  for (I=0; I<NA; I++) {
    for (IC=M1-1; IC<M2; IC++) {
      IR1=IC-NOV+1;
      IRP=IR1-1;
      IRFIRP = IRF(IRP,I);
      ICFIC  = ICF(IC,I);
      for (IR=IR1; IR<NRA; IR++) {
        IRFIR=IRF(IR,I);
        if (A(ICFIC,IRFIR,I) != 0.0) {
	  FA(I,IRFIR)-=A(ICFIC,IRFIR,I)*FA(I,IRFIRP);
	}
      }
      for (IR=NBCP1-1; IR<NRC; IR++) {
        if (C(ICFIC,IR,I) != 0.0) {
	  FC(IR)-=C(ICFIC,IR,I)*FA(I,IRFIRP);
	}
      }
    }
  }
  return 0;
}

#define FAA(i,j)	MATRIXELEMFAA(mp,i,j)

Entry(Int2) bvp_cpyrhs(bvpMatrixPtr mp, Vector vp) {
/****   COPIES THE RHS ***/
Int2 I,IR,IRFIR;
      
  for (I=0; I<NA; I++) {
    for (IR=0; IR<NOV; IR++) {
      IRFIR=IRF(NRA-NOV+IR,I);
      FAA(IR,I)=FA(I,IRFIR);
    }        
  }
  return 0;
}

Entry(Int2) bvp_redrhs(bvpMatrixPtr mp, Vector vp) {
Int2 NAM1;
Int2 I1,I2,IC,IR,ICP1,IPIV1,L1;
VAR  TMP,RM;

  NAM1    = NA-1;

  for (I1=0; I1<NAM1; I1++) { 
    I2=I1+1;
    for (IC=0; IC<NOV; IC++) {
      ICP1  = IC+1;
      IPIV1 = IPR(IC,I1);
      if (IPIV1 < NOV) {
        TMP           = FAA(IC,I1);
        FAA(IC,I1)    = FAA(IPIV1,I1);
        FAA(IPIV1,I1) = TMP;
      } else {
        L1         = IPIV1-NOV;
        TMP        = FAA(IC,I1);
        FAA(IC,I1) = FAA(L1,I2);
        FAA(L1,I2) = TMP;
      }
      for (IR=ICP1; IR<NOV; IR++) {
        L1=ICF2(IC,I1);
        RM=A2(IR,L1,I1);
        FAA(IR,I1)=FAA(IR,I1)-RM*FAA(IC,I1);
      }
      for (IR=0; IR<NOV; IR++) {
        L1=ICF1(IC,I2);
        RM=A1(IR,L1,I2);
        FAA(IR,I2) = FAA(IR,I2)-RM*FAA(IC,I1);
      }
      for (IR=NBC; IR<NRC; IR++) {
        L1=ICF2(IC,I1);
        RM=CC(L1,IR,I2);
        FC(IR)= FC(IR)-RM*FAA(IC,I1);
      }
    }
  }
  return 0;
}            

#define N       (NRC+NOV)
#define U(i)    MATRIXELEMFCC(mp,i)
#define F(i)    MATRIXELEMBUF(mp,i)
#define IC(i)	MATRIXELEMIC(mp,i)
#define IR(i)	MATRIXELEMIR(mp,i)

Local(Int2) ge(bvpMatrixPtr mp, VAR *DET) {
Int2 I,J,JJ,IPIV,JPIV,K,JJP1,IP1,I1,NM1,L;
VAR P,PIV,RM,SM;

/* Initialization */

       for (I=0; I<N; I++) {
         IC(I)=I;
         IR(I)=I;
       }

/*   Elimination  */

       *DET=1.0;
       NM1=N-1;

       for (JJ=0; JJ<NM1; JJ++) {
         IPIV=JJ;
         JPIV=JJ;
         PIV=0.0;
	 for (I=JJ; I<N; I++) {
	   for (J=JJ; J<N; J++) {
             P=fabs(E(IR(I),IC(J)));
	     if (P > PIV) {
               PIV=P;
               IPIV=I;
               JPIV=J;
	     }
           }
         }

         *DET=*DET*E(IR(IPIV),IC(JPIV));
	 if (IPIV != JJ) *DET=-*DET;
	 if (JPIV != JJ) *DET=-*DET;
 
         K=IR(JJ);
         IR(JJ)=IR(IPIV);
         IR(IPIV)=K;

         K=IC(JJ);
         IC(JJ)=IC(JPIV);
         IC(JPIV)=K;

         JJP1=JJ+1;
	 for (L=JJP1; L<N; L++) {
           RM=E(IR(L),IC(JJ))/E(IR(JJ),IC(JJ));
	   for (I=JJP1; I<N; I++) 
	     E(IR(L),IC(I))-=RM*E(IR(JJ),IC(I));
           F(IR(L))-=RM*F(IR(JJ));
         }
       }
       *DET=*DET*E(IR(N-1),IC(N-1));

/*   Backsubstitution:  */

         U(IC(N-1))=F(IR(N-1))/E(IR(N-1),IC(N-1));
	 for (I1=0; I1<NM1; I1++) {
           I=N-2-I1;
           IP1=I+1;
	   for (J=IP1, SM=0.0; J<N; J++) SM+=E(IR(I),IC(J))*U(IC(J));
           U(IC(I))=(F(IR(I))-SM)/E(IR(I),IC(I));
         }
  return 0;
}


Local(Int2) nlvc(bvpMatrixPtr mp, Int2 K) {
Int2 I,J,JJ,IPIV,JPIV,JJP1,IP1,I1,NMK,L,KK;
VAR P,PIV,RM,SM;
/*
 Finds a null-vector of a singular (N x N)-matrix E.
 The null space of E is assumed to be K-dimensional.
 On exit U contains the null vector
*/

/* Initialization */

       for (I=0; I<N; I++) {
         IC(I)=I;
         IR(I)=I;
       }

/*   Elimination  */

       NMK=N-K;

       for (JJ=0; JJ<NMK; JJ++) {
         IPIV=JJ;
         JPIV=JJ;
         PIV=0.0;
	 for (I=JJ; I<N; I++) {
	   for (J=JJ; J<N; J++) {
             P=fabs(E(IR(I),IC(J)));
	     if (P > PIV) {
               PIV=P;
               IPIV=I;
               JPIV=J;
	     }
           }
         }

         KK=IR(JJ);
         IR(JJ)=IR(IPIV);
         IR(IPIV)=KK;

         KK=IC(JJ);
         IC(JJ)=IC(JPIV);
         IC(JPIV)=KK;

         JJP1=JJ+1;
	 for (L=JJP1; L<N; L++) {
           RM=E(IR(L),IC(JJ))/E(IR(JJ),IC(JJ));
	   for (I=JJP1; I<N; I++) 
	     E(IR(L),IC(I))-=RM*E(IR(JJ),IC(I));
         }
       }
              
/*   Backsubstitution:  */
         for (I=0; I<K; I++) U(IC(N-1-I))=1.0;
	 for (I1=0; I1<NMK; I1++) {
           I=NMK-1-I1;
           IP1=I+1;
	   for (J=IP1, SM=0.0; J<N; J++) SM+=E(IR(I),IC(J))*U(IC(J));
           U(IC(I))=-SM/E(IR(I),IC(I));
         }
  return 0;
}

#define XE(i)	 	MATRIXELEMBUF(mp,i)
#define FCC(i)		MATRIXELEMFCC(mp,i)

Entry(Int2) bvp_dimrge(bvpMatrixPtr mp, Vector vp, Int2 NILV) {
Int2 NCR,I,J,K,NOVPJ,NOVPJ2,NOVPI,K1,K2,KR,KC;
VAR DET;

  NCR     = NRC+NOV;

/*  COPYING  */
  for (I=0; I<NOV; I++) {
    for (J=0; J<NOV; J++) {
      NOVPJ      = NOV+J;
      E(I,J)     = S1(I,J,NA-1);
      E(I,NOVPJ) = A2(I,J,NA-1);
    }
    for (J=0; J<NCB; J++) {
      NOVPJ2      = 2*NOV+J;
      E(I,NOVPJ2) = BB(I,J,NA-1);
    }
  }

  for (I=0; I<NRC; I++) {     
    NOVPI=NOV+I;
    for (J=0; J<NOV; J++) {
      NOVPJ          = NOV+J;
      E(NOVPI,J)     = CC(J,I,0);
      E(NOVPI,NOVPJ) = CC(J,I,NA);
    }
    for (J=0; J<NCB; J++) {
      NOVPJ2          = 2*NOV+J;
      E(NOVPI,NOVPJ2) = D(J,I);
    }
  }

  for (I=0; I<NOV; I++) XE(I)=FAA(I,NA-1);
  for (I=0; I<NRC; I++) XE(NOV+I) = FC(I);


/* bvp_PrintMatrix(mp,"REDUCED SYSTEM:",4); */

/*  SOLVE FOR FCC  */
  if (NILV == 0) {
    ge (mp, &DET);
  } else {
    nlvc(mp, NILV);
    DET=0.0;
  }


/* ????? */
  K1=NCR;
  K2=K1+NOV*NOV;
  for (KR=0; KR<NOV; KR++) {
    for (KC=0; KC<NOV; KC++) {
      K=NOV*KR+KC;
      FCC(K1+K)=P0(KR,KC);
      FCC(K2+K)=P1(KR,KC);
    }
  }
  FCC(NCR+2*NOV*NOV)=DET;
/* bvp_PrintMatrix(mp,"SOLUTION VECTOR :",6); */
  
  for (I=0; I<NRC; I++) FC(I)=FCC(NOV+I);
  return 0;
}

#define SOL1(i,j)	MATRIXELEMSOL1(mp,i,j)
#define SOL2(i,j)	MATRIXELEMSOL2(mp,i,j)
#define SOL3(i,j)	MATRIXELEMSOL3(mp,i,j)

Entry(Int2) bvp_bcksub(bvpMatrixPtr mp, Vector vp) {

Int2 L,I,K,L1,L2,NAM1;
VAR SM;
      
/*  BACKSUBSTITUTE PROCESS FOR REDUCE IS A RECURSIVE PROCEDURE  */

  for (L=0; L<NOV; L++) {
    SOL1(L,NA-1) = FCC(L);
    SOL2(L,NA-1) = FC(L);
  }
  for (L=0; L<NOV; L++) SOL2(L,NA-1)=FC(L);
  if (NA > 1) {
    for (L=0; L<NOV; L++) {
      SOL1(L,NA-2)=SOL1(L,NA-1);
      SOL3(L,NA-2)=SOL2(L,NA-1);
    }
  }
  
/* BACKSUBSITUTE PROCESS  */

  NAM1=NA-1;
  for (I=NAM1-1; I>=0; I--) {
    for (K=NOV-1; K>=0; K--) {
       for (L=0, SM=0.0; L<NOV; L++) {
          SM+=SOL1(L,I)*S1(K,L,I);
          SM+=SOL3(L,I)*S2(K,L,I);
       }
       for (L=0; L<NCB; L++) SM+=FC(NOV+L)*BB(K,L,I);
       for (L=K+1; L<NOV; L++) {
         L1=ICF2(L,I);
         SM+=SOL2(L1,I)*A2(K,L1,I);
       }
       L2=ICF2(K,I);
       SOL2(L2,I)=(FAA(K,I)-SM)/A2(K,L2,I);
    }
    for (L=0; L<NOV; L++) {
      SOL1(L,I+1)=SOL2(L,I);
      if (I > 0) {
        SOL3(L,I-1)=SOL2(L,I);
        SOL1(L,I-1)=SOL1(L,I);            
      }
    }
  }
  return 0;
}

#define X(i)	MATRIXELEMX(mp,i)

Entry(Int2) bvp_infpar(bvpMatrixPtr mp, Vector vp) {
Int2  I,IR,J,NRAM,NRAPJ,IRP1,IRFIR,NOVPJ,NOVPIR,J1,ICFNOVPIR,ICFJ1;
VAR   SM;
/*   DETERMINE THE "LOCAL" VARIABLES BY BACKSUBSTITUTION   */

  NRAM=NRA-NOV;

/*   BACKSUBSTITUTE PROCESS IN THE CONDENSATION OF PARAMETERS   */

  for (I=0; I<NA; I++) {
    for (IR=NRAM-1; IR>=0; IR--) {
      IRP1=IR+1;
      IRFIR=IRF(IR,I);
      for (J=0, SM=0.0; J<NOV; J++) {
        NRAPJ=NRA+J;
        SM+=A(J,IRFIR,I)*SOL1(J,I);
        SM+=A(NRAPJ,IRFIR,I)*SOL2(J,I);
      }
      for (J=0; J<NCB; J++) {
        NOVPJ=NOV+J;
        SM+=B(J,IRFIR,I)*FC(NOVPJ);
      }
      for (J=IRP1; J<NRAM; J++) {
        J1=J+NOV;
        ICFJ1=ICF(J1,I);
        SM+=A(ICFJ1,IRFIR,I)*X(ICFJ1);
      }
      NOVPIR=NOV+IR;
      ICFNOVPIR=ICF(NOVPIR,I);
      X(ICFNOVPIR)=(FA(I,IRFIR)-SM)/A(ICFNOVPIR,IRFIR,I);
    }
	 
/*  COPY THE SOLUTION STORED IN SOL1 AND X INTO FA  */

    for (J=0; J<NOV; J++) FA(I,J)=SOL1(J,I);
    for (J=NOV; J<NRA; J++) FA(I,J)=X(J);
  }
  return 0;
}


Entry(Int2) bvp_Decomp(bvpMatrixPtr mp DEBPAR) {
   ndecomp++;
/* printf("ndecomp=%i, nback=%i\n",ndecomp, nback); */
/* translation of BRBDD from AUTO */
  bvp_conpar(mp);  
/* bvp_PrintMatrix(mp,"After conpar:",0); */ 
  bvp_copycp(mp);
/* bvp_PrintMatrix(mp,"Before reduce:",1); */
  bvp_reduce(mp);
/* bvp_PrintMatrix(mp,"After reduce:",1); */
  return 0;
}

Entry(void) bvp_Solve(bvpMatrixPtr mp, Vector vp) {
   nback++;
/*   printf("ndecomp=%i, nback=%i\n",ndecomp, nback); */
/* translation of BRBDB from AUTO */

/* bvp_PrintVector(vp,"Before conrhs:"); */
  bvp_conrhs(mp, vp);
/*  bvp_PrintVector(vp,"After conrhs:"); */
  bvp_cpyrhs(mp, vp);
/*  bvp_PrintMatrix(mp,"After cpyrhs:",3); */
  bvp_redrhs(mp, vp);
/*  bvp_PrintVector(vp,"After redrhs:");  */
/*  bvp_PrintMatrix(mp,"After redrhs:",1); */
  bvp_dimrge(mp, vp, 0);
/*  bvp_PrintMatrix(mp,"After dimrge:",6); */
/*  bvp_PrintVector(vp,"After dimrge:"); */
  bvp_bcksub(mp, vp);
/*  bvp_PrintMatrix(mp,"After bcksub:",7); */
  bvp_infpar(mp, vp);
  
  return;
}

Entry(VAR) bvp_NormOfVector(Vector a) {
VAR s;
  s=sqrt(stdalg->ScalProd(a,a)/(NTST*NCOL));
  if (s>=1.0) s=1.0;
  return s;
}

Entry(void) bvp_setzero(bvpMatrixPtr mp, Vector vp) {
Int2 I,J;
  for (I=0; I<NA; I++) {
    for (J=0; J<NRA; J++) FA(I,J)=0.0;
  }
  for (I=0; I<NRC; I++) FC(I)=0.0;
  return;
}

Entry(void) bvp_SingSolve(bvpMatrixPtr mp, Vector vp) {
/* translation of BRBDN from AUTO */
  bvp_conpar(mp);  
/* bvp_PrintMatrix(mp,"After conpar in SingSolve:",0);  */
  bvp_copycp(mp);
/* bvp_PrintMatrix(mp,"Before reduce in SingSolve:",1); */
  bvp_reduce(mp);
/* bvp_PrintMatrix(mp,"After reduce in SingSolve:",1); */
  bvp_setzero(mp, vp);
/*  bvp_PrintVector(vp,"After setzero in SingSolve:"); */
  bvp_dimrge(mp, vp, 1);
/*  bvp_PrintMatrix(mp,"After dimrge in SingSolve:",6); */
/*  bvp_PrintVector(vp,"After dimrge in SingSolve:"); */
  bvp_bcksub(mp, vp);
/*  bvp_PrintMatrix(mp,"After bcksub in SingSolve:",7); */
  bvp_infpar(mp, vp);

  return;
}

#undef A
#undef B
#undef C
#undef D
#undef FA
#undef FC

#undef NDIM 
#undef NTST 
#undef NCOL 
#undef NFPR 

#undef IRF
#undef ICF
#undef A1
#undef A2
#undef BB
#undef CC
#undef S1
#undef S2
#undef ICF1
#undef ICF2

#undef FAA
#undef FCC
#undef XE
#undef E
#undef P0
#undef P1

#undef X

#undef NOV
#undef NA
#undef NRA
#undef NCA
#undef NCB
#undef NBC
#undef NRC
#undef VECTNUM