/*
   m__pd1
   Flip (Period doubling) Curve for maps using bordering method
*/

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

#include "stagen.h"

#include "_util.h"
#include "_linalg.h"

#include "m__pd1.h"
#include "_conti.h"	/* continuer's (i.e. generator's) data */

/* Numbers of visible name classes;
   the order must be the same as in /names section */
#define CL_TIME		1
#define CL_PHASE	2
#define CL_PAR		3
#define CL_TEST         4
#define CL_UTEST	5
#define CL_STDFUN	6
#define CL_MULT		7

#define _eps     staData->eps
#define _dx      staData->dx
#define _branch  staData->branch   

#define PI	4.0*atan(1.0)
#define PID2	2.0*atan(1.0)   

Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(m__pd_DataPtr) staData=NULL;	/* global copy of our own data area */
Local(StdLinAlgDataPtr) staFunc=NULL;	/* global copy of pointer standard linear algebra struct */
Local(UtilitiesDataPtr) staUtil=NULL;	/* global copy of pointer to utility functions struct */
Local(Int2Ptr) active=NULL;        	/* index of active parameter(s) */
Local(Int2) nappend;			/* number of appended user functions */ 
Local(Int2) npar,nphase;		/* dims of par and phase spaces */
Local(Vector) x0=NULL;			/* approximate first point */
Local(Vector) x1=NULL;			/* current point */
Local(Vector) v0=NULL;			/* approximate tangent vector at first point */
Local(FloatHiPtr) X=NULL;		/* current set of phases used by Decoder */
Local(FloatHiPtr) P=NULL;		/* current set of params used by Decoder */
Local(FloatHiPtr) M=NULL;		/* current set of multipliers */
Local(VAR) zero=0;
Local(Vector) StdFun=NULL;
Local(int) itmap;
Local(FloatHiPtr) F,DF,DFF,D2,DD2,D3F,DD3F;		/* work space for supdiff */

#define _itmap       staData->itmap   	/* number of superpositions */

#include "_supdiff.c"

/* Working space for defining and test functions  */

Local(Matrix) A0,A,B,C,C1,C2,AA,AT,AH,BB,BP,BT,BIP,D,DE;
Local(Vector) vn,xp,xp1,V,W,F,Vp,L,VV,VB,V11,W11,Q,Qr,Qi,PP,Pr,Pi,aa,bb,cc,dd,rr;
Local(Vector) F0,F1,F2,a,x01,v01,x11;
Local(Vector) Re=NULL,Im=NULL,Mod=NULL,Arg=NULL;
Local(VAR) test[4],dx,ddx,dddx;

#define BifVector stagenData->BifDataOutPtr
#define ActiveUTestFuns     staData->ActiveUserTest
#define nUTest      	    stagenData->udFuncNum

/* Initialize function common for Starter and Decoder */
Local(Int2) Init(StagenDataPtr stagenptr) {
Int2 i,m;
  /* Set local pointers to data area and functions array */
  stagenData=stagenptr;
  staData=(m__pd_DataPtr)stagenptr->StaPtr;
  /* Use standard linear algebra */
  staFunc=(StdLinAlgDataPtr)stagenptr->StdLinAlg;
  staUtil=(UtilitiesDataPtr)stagenptr->Utilities;
  /* Get dimensions of Phase and Par spaces */
  npar=*stagenptr->ClassDim[CL_PAR];
  nphase=*stagenptr->ClassDim[CL_PHASE];
  itmap=_itmap;
  if (itmap<=0) {
    staUtil->ErrMessage("Superposition must be positive.");
    return 6;
  }
  /* Call udf SetInitPoint function, if any */
  if (!stagenData->ContDataGenLen && FuncParActive(staData->SetInitPoint)) {	/* there is UDF to set InitPoint */
    Int2 Index_X,Index_P;
    Index_X=stagenptr->ParIndex(staData,X);
    Index_P=stagenptr->ParIndex(staData,P);
    FuncParCall(staData->SetInitPoint)(staData);
    stagenptr->UpdatePar(Index_X);
    stagenptr->UpdatePar(Index_P);
  }
  /* Do all the checks related to active parameters and find them */
  for (i=0,nappend=0; i<nUTest; i++) if (ActiveUTestFuns[i]==UDF_APPEND) nappend++;
  active=MemNew(sizeof(Int2)*(nappend+2));    
  if (staUtil->CheckParameters(staData->Active,npar,active,2+nappend)) {
    active=MemFree(active);
    return 1;
  }
  /* Checks related to analytical/numerical differentiation of RHS */
  if (!stagenData->Der1 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Jacoby matrix is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 2;
  }
  if ((staData->ActiveSing[0] || staData->ActiveSing[1] || staData->ActiveSing[2]) && 
      !stagenData->Der2 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Hessian matrix is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 3;
  }
  /* Create a copy parameters */
  m=nphase*(nphase-1)/2;
  Der2num=ind2_(nphase+npar-1,nphase+npar-1)+1;
  Der3num=ind3_(nphase+npar-1,nphase+npar-1,nphase+npar-1)+1;
  if (!P) {
    P=MemNew(sizeof(FloatHi)*npar);
    MemCopy(P,staData->P,sizeof(FloatHi)*npar);
    /* Allocate memory for copy of phase vars used by Decoder */
    X=MemNew(sizeof(FloatHi)*nphase);
    /* Allocate memory for eigenvalues */
    M=MemNew(sizeof(FloatHi)*4*nphase);
    Re=staFunc->CreateVector(nphase);
    Im=staFunc->CreateVector(nphase);
    Mod=staFunc->CreateVector(nphase);
    Arg=staFunc->CreateVector(nphase);   
    /* Allocate memory for working space used in defining and test functions */
    A=staFunc->CreateMatrix(nphase,nphase);
    A0=staFunc->CreateMatrix(nphase,nphase+2+nappend);    
    B=staFunc->CreateMatrix(nphase,nphase);
    BP=staFunc->CreateMatrix(nphase+1,nphase+1);
    BT=staFunc->CreateMatrix(nphase+1,nphase+1);
    BIP=staFunc->CreateMatrix(m,m);
    V11=staFunc->CreateVector(nphase);
    W11=staFunc->CreateVector(nphase);
    VB=staFunc->CreateVector(nphase+1);
    VV=staFunc->CreateVector(2*nphase);
    Q=staFunc->CreateVector(nphase);
    PP=staFunc->CreateVector(nphase+1);

    AA=staFunc->CreateMatrix(nphase,nphase);
    BB=staFunc->CreateMatrix(2*nphase,2*nphase);
    C=staFunc->CreateMatrix(nphase+nappend,nphase+2+nappend);
    C1=staFunc->CreateMatrix(nphase,nphase+2+nappend);
    C2=staFunc->CreateMatrix(1,nphase+2+nappend);
    x0=staFunc->CreateVector(nphase+2+nappend);
    x1=staFunc->CreateVector(nphase+2+nappend);
    v0=staFunc->CreateVector(nphase+2+nappend);
    x01=staFunc->CreateVector(nphase+2+nappend);
    x11=staFunc->CreateVector(nphase+2+nappend);
    v01=staFunc->CreateVector(nphase+2+nappend);
    vn=staFunc->CreateVector(nphase);
    xp=staFunc->CreateVector(nphase+2+nappend);
    xp1=staFunc->CreateVector(nphase+2+nappend);
    StdFun=staFunc->CreateVector(1+2*nphase); /* L2_NORM, MIN&MAX */ 
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    L=staFunc->CreateVector(nphase);
    Vp=staFunc->CreateVector(nphase+2+nappend);
    BifVector=MemNew(sizeof(FloatHi)*(3*nphase+2));
    if (nphase == 1) staData->ActiveSing[1]=FALSE;
    if (nphase <= 2) staData->ActiveSing[2]=FALSE; 
    F=MemNew(sizeof(FloatHi)*nphase);  
    F0=staFunc->CreateVector(nphase);
    F1=staFunc->CreateVector(nphase);
    F2=staFunc->CreateVector(nphase);
    a=staFunc->CreateVector(nphase);
    DF=MemNew(sizeof(FloatHi)*nphase*(nphase+npar));  
    DFF=MemNew(sizeof(FloatHi)*nphase*(nphase+npar));  
    if (staData->ActiveSing[0]+staData->ActiveSing[1]+staData->ActiveSing[2]+staData->ActiveSing[3]) {
      D2=MemNew(sizeof(FloatHi)*nphase*Der2num);
      DD2=MemNew(sizeof(FloatHi)*nphase*Der2num);
      if (staData->ActiveSing[2]+staData->ActiveSing[3]) {
	D3F=MemNew(sizeof(FloatHi)*nphase*Der3num);
	DD3F=MemNew(sizeof(FloatHi)*nphase*Der3num);
      }
    }
  }
  dx=_dx;
  ddx=pow(dx,3.0/4.0);
  dddx=pow(dx,3.0/6.0);

  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  
  return 0;
}

Local(Int2) JacRHS (Vector x, Matrix J);

Local(void) Term (Boolean OK) {
Int2 m;
Int2 i,j,imax,jmax;
VAR d,a11,a12,a21,a22;
  m=nphase*(nphase-1)/2;
  if (OK&&P&&(stagenData->Reason==RF_CLEAN)) {
/* saving global vectors V11 and W11 used in defining functions */
    stagenData->ContDataStaLen=2*nphase*sizeof(FloatHi);
    stagenData->ContDataStaPtr=MemNew(stagenData->ContDataStaLen);
    MemCopy(stagenData->ContDataStaPtr,V11,nphase*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+nphase,W11,nphase*sizeof(FloatHi));
    stagenData->BifDataOutLen=sizeof(FloatHi)*(nphase+1);
    BifVector[0]=itmap;
    MemCopy(BifVector+1,V11,sizeof(FloatHi)*nphase); 
    stagenData->ProcessPoint(stagenData->ProcFuncNum, x1, "Last point");
  } else {
    stagenData->ContDataStaLen=0;
    stagenData->BifDataOutLen=0;
  }
  x0=staFunc->DestroyVector(x0);
  x1=staFunc->DestroyVector(x1);
  v0=staFunc->DestroyVector(v0);
  active=MemFree(active);
  if (P) {
    P=MemFree(P);
    X=MemFree(X);
    M=MemFree(M);
    Re=staFunc->DestroyVector(Re);
    Im=staFunc->DestroyVector(Im);
    Mod=staFunc->DestroyVector(Mod);
    Arg=staFunc->DestroyVector(Arg); 
    StdFun=staFunc->DestroyVector(StdFun);    
    A=staFunc->DestroyMatrix(A);
    A0=staFunc->DestroyMatrix(A0);
    B=staFunc->DestroyMatrix(B);
    BP=staFunc->DestroyMatrix(BP);
    BT=staFunc->DestroyMatrix(BT);
    BIP=staFunc->DestroyMatrix(BIP);
    V11=staFunc->DestroyVector(V11);
    W11=staFunc->DestroyVector(W11);
    Q=staFunc->DestroyVector(Q);
    PP=staFunc->DestroyVector(PP);
    F0=staFunc->DestroyVector(F0);
    F1=staFunc->DestroyVector(F1);
    F2=staFunc->DestroyVector(F2);
    a=staFunc->DestroyVector(a);

    AA=staFunc->DestroyMatrix(AA);
    BB=staFunc->DestroyMatrix(BB);
    VV=staFunc->DestroyVector(VV);
    VB=staFunc->DestroyVector(VB);
    C=staFunc->DestroyMatrix(C);
    C1=staFunc->DestroyMatrix(C1);
    C2=staFunc->DestroyMatrix(C2);
    vn=staFunc->DestroyVector(vn);
    xp=staFunc->DestroyVector(xp);
    xp1=staFunc->DestroyVector(xp1);
    x01=staFunc->DestroyVector(x01);
    x11=staFunc->DestroyVector(x11);
    v01=staFunc->DestroyVector(v01);
    V=staFunc->DestroyVector(V);
    W=staFunc->DestroyVector(W);
    Vp=staFunc->DestroyVector(Vp);
    L=staFunc->DestroyVector(L);
    BifVector=MemFree(BifVector);    
    F=MemFree(F);
    DF=MemFree(DF);
    DFF=MemFree(DFF);
    if (staData->ActiveSing[0]+staData->ActiveSing[1]+staData->ActiveSing[2]+staData->ActiveSing[3]) {
      D2=MemFree(D2);
      DD2=MemFree(DD2);
      if (staData->ActiveSing[2]+staData->ActiveSing[3]) {
	D3F=MemFree(D3F);
	DD3F=MemFree(DD3F);
      }
    }
  }
} 

Local(Int2) JacDefFixed (Vector x, Matrix J);
Local(Int2) JacDefEig(Vector x, Matrix J);
Entry(Int2) mpd1_JacDefFlip (Vector x, Matrix J);
Entry(Int2) mpd1_DefFlip (Vector x, Vector J);


/***************************************/
/* Any point on LP curve -> pd Starter */

Entry(Int2) mpd1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,Ret,m,p,q,r,s;
  Vector v1,v2;
  Matrix D,DE;
  VAR t;

  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  if (Init(stagenptr)) return 1;	/* initialize local variables */
  m=nphase*(nphase-1)/2;
  
  /* Create the first point for generator */
  if (stagenData->BifDataInOk) _itmap=stagenData->BifDataInPtr[0];  
  if (stagenptr->ContDataGenLen) {       /* "Extension" */
    MemCopy(V11,stagenptr->ContDataStaPtr,sizeof(FloatHi)*nphase);
    MemCopy(W11,stagenptr->ContDataStaPtr+nphase,sizeof(FloatHi)*nphase);
    stagenptr->ContDataStaPtr=MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
  } else {
    v1=staFunc->CreateVector(nphase+2+nappend);
    v2=staFunc->CreateVector(nphase+2+nappend);
    D=staFunc->CreateMatrix(nphase+1+nappend,nphase+2+nappend);
    DE=staFunc->CreateMatrix(nphase+2+nappend,nphase+2+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0;
    staFunc->ClearMatrix(BP);
    staFunc->CopyMatrixBlock(A,0,0,BP,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP[i][nphase]=i;
      BP[nphase][i]=i;
    }
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<=nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP[i][j]);
          }           
	}
      }
      if (k!=q) { /* switch colom k with colom  q */
        for (i=0;i<=nphase;i++) {
          t=BP[i][k];
          BP[i][k]=BP[i][q];
          BP[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */
        for (i=0;i<=nphase;i++) {
          t=BP[k][i];
          BP[k][i]=BP[p][i];
          BP[p][i]=t;
	}
      }    
      if (fabs(BP[k][k]) == 0.0) {
        staUtil->ErrMessage("Matrix has rank defect at least two");
        Term(FALSE);
        break;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP[i][k]/BP[k][k];
        for (j=kp1;j<nphase;j++) {
          BP[i][j]+=t*BP[k][j];
	}
      }
    } 
    p=(Int2)BP[nphase][nphase-1];
    q=(Int2)BP[nphase-1][nphase];
    staFunc->ClearVector(V11);
    V11[p]=1.0;
    staFunc->ClearVector(W11); 
    W11[q]=1.0;

/* perform one adaptation step for V11 and W11 */
    staFunc->ClearMatrix(BP);
    staFunc->CopyMatrixBlock(C1,0,0,BP,0,0,nphase,nphase); 
    for (i=0; i<nphase; i++) {
      BP[i][i]+=1.0;
      BP[nphase][i]=V11[i];
      BP[i][nphase]=W11[i];
    }
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) { 
      goto final;
    }    
    staFunc->ClearVector(VB); 
    VB[nphase]=1.0; 
    staFunc->Solve(BT,VB);
    staFunc->CopyVectorElements(VB,0,V11,0,nphase);   
    staFunc->NormalizeVector(V11);    
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) { 
      goto final;
    }
    staFunc->ClearVector(PP);
    PP[nphase]=1.0;
    staFunc->Solve(BT,PP);
    staFunc->CopyVectorElements(PP,0,W11,0,nphase);   
    staFunc->NormalizeVector(W11);
final:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
 
  /* Compute a tangent vector to the Flip curve */
    mpd1_DefFlip(x0,v0);
    Ret=mpd1_JacDefFlip(x0,D);

    for (i=0; i<nphase+2+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+1+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Flip curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }    
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  return 0;
}

/***************************************/
/* FF Point -> pd Starter */
Entry(Int2) mf_pd1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,Ret,p,q;
  Vector v1,v2;
  Matrix D,DE;
  VAR t;

  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  if (Init(stagenptr)) return 1;	/* initialize local variables */
  
  /* Create the first point for generator */

  if (stagenData->BifDataInOk) _itmap=stagenData->BifDataInPtr[0];  
  if (stagenptr->ContDataGenLen) {       /* "Extension" */
    MemCopy(V11,stagenptr->ContDataStaPtr,sizeof(FloatHi)*nphase);
    MemCopy(W11,stagenptr->ContDataStaPtr+nphase,sizeof(FloatHi)*nphase);
    stagenptr->ContDataStaPtr=MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
  } else {
    v1=staFunc->CreateVector(nphase+2+nappend);
    v2=staFunc->CreateVector(nphase+2+nappend);
    D=staFunc->CreateMatrix(nphase+1+nappend,nphase+2+nappend);
    DE=staFunc->CreateMatrix(nphase+2+nappend,nphase+2+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0; 
    staFunc->ClearMatrix(BP);
    staFunc->CopyMatrixBlock(A,0,0,BP,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP[i][nphase]=i;
      BP[nphase][i]=i;
    }
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<=nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP[i][j]);
          }           
	}
      }
      if (k!=q) { /* switch colom k with colom  q */
        for (i=0;i<=nphase;i++) {
          t=BP[i][k];
          BP[i][k]=BP[i][q];
          BP[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */
        for (i=0;i<=nphase;i++) {
          t=BP[k][i];
          BP[k][i]=BP[p][i];
          BP[p][i]=t;
	}
      }    
      if (fabs(BP[k][k]) == 0.0) {
        staUtil->ErrMessage("Matrix has rank defect at least two");
        Term(FALSE);
        break;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP[i][k]/BP[k][k];
        for (j=kp1;j<nphase;j++) {
          BP[i][j]+=t*BP[k][j];
	}
      }
    } 
    p=(Int2)BP[nphase][nphase-1];
    q=(Int2)BP[nphase-1][nphase];
    staFunc->ClearVector(V11); 
    V11[p]=1.0;
    staFunc->ClearVector(W11); 
    W11[q]=1.0;

/* perform one adaptation step for V11 and W11 */
    staFunc->ClearMatrix(BP);
    staFunc->CopyMatrixBlock(C1,0,0,BP,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP[i][i]+=1.0;
      BP[nphase][i]=V11[i];
      BP[i][nphase]=W11[i];
    }
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) { 
       goto final;
    }    
    staFunc->ClearVector(VB); 
    VB[nphase]=1.0; 
    staFunc->Solve(BT,VB);
    staFunc->CopyVectorElements(VB,0,V11,0,nphase);   
    staFunc->NormalizeVector(V11);    
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) { 
      goto final;
    }
    staFunc->ClearVector(PP);
    PP[nphase]=1.0;
    staFunc->Solve(BT,PP);
    staFunc->CopyVectorElements(PP,0,W11,0,nphase);   
    staFunc->NormalizeVector(W11);
final:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
 
  /* Compute a tangent vector to the Flip curve */
    mpd1_DefFlip(x0,v0);
    Ret=mpd1_JacDefFlip(x0,D);
    for (i=0; i<nphase+2+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+1+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Flip curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);

    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }    
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  return 0;
}

/*************************************************************/
/* Defining conditions for the eigenvector: (A*A+kappa*E)V=0 */
Entry(Int2) DefEig (Vector x, Vector y) {

  if (JacRHS(x,C1)) return 1;	    /* A,C1,V,W - global */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

  staFunc->CopyVectorElements(x, nphase+2+nappend, V, 0, nphase);
  staFunc->MultiplyMatrixVector(A,V,W);
  staFunc->MultiplyMatrixVector(A,W,V);
  staFunc->CopyVectorElements(x, nphase+2+nappend, W, 0, nphase);
  staFunc->AddScaledVector(V,W,x[2*nphase+2+nappend],V);
  staFunc->CopyVectorElements(V, 0, y, 0, nphase);

  return 0;
}

/****************************************************/
/* Jacobian of defining conditions for eigenvectors */
Local(Int2) JacDefEig (Vector x, Matrix Jac) {
int i,j,k,l;
double s;

  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2) {
    /* creation of DefEig Jacobian matrix Jac using symbolic second derivatives */
    H=SymHess(x);
    if (JacRHS(x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      for (l=0; l<nphase+2+nappend; l++) {
        for (s=0.0,j=0; j<nphase; j++) {
          for (k=0; k<nphase; k++)
            s+=(HessianElem(i,l,j)*A[j][k]+HessianElem(j,l,k)*A[i][j])*x[k+nphase+2+nappend];       
        }
        Jac[i][l]=s;
      }
    }
    for (i=0; i<nphase; i++) {
      for (l=0; l<nphase; l++) {
        for (s=0.0,j=0; j<nphase; j++) s+=A[i][j]*A[j][l];
        if (i==l) s+=x[2*nphase+2+nappend];
        Jac[i][l+nphase+2+nappend]=s;
      }
      Jac[i][2*nphase+2+nappend]=x[i+nphase+2+nappend];
    }
  } else {
    if (staFunc->Der(DefEig, nphase, x, _dx, Jac)) return 1;
  }
  return(0);
}

Local(Int2) JacDefFixed (Vector x, Matrix Jac);

/************************/
/* Some common routines */
Local(void) CopyToX(Vector vp) {
  MemCopy(X,vp,sizeof(FloatHi)*nphase);  
}

Local(void) CopyToP(Vector vp) {
Int2 i;
  for (i=0; i<2+nappend; i++) P[active[i]]=vp[nphase+i];
}

Local(void) CopyToM(void) {
  MemCopy(M,Mod,sizeof(VAR)*nphase);
  MemCopy(M+nphase,Arg,sizeof(VAR)*nphase);
  MemCopy(M+2*nphase,Re,sizeof(VAR)*nphase);
  MemCopy(M+3*nphase,Im,sizeof(VAR)*nphase);
}

/***********/
/* Decoder */

Entry(Int2) mpd1_Decoder(VoidPtr vpoint, FloatHiPtr PNTR indirect, Int2 oper) {
  switch (oper) {
    case DECODER_INIT:	/* point is a pointer to StagenData structure */
    case DECODER_TERM:
     break;
    case DECODER_DIM:	/* returns number of M-points in given G-point */
      return 1;
    default:		/* extract next M-point from G-point pointed by vpoint */
      /* vpoint is a Vector */
      indirect[CL_TIME]=&zero;
      indirect[CL_PHASE]=X;
      CopyToX(vpoint);
      indirect[CL_PAR]=P;
      CopyToP(vpoint);
      indirect[CL_MULT]=M;
      CopyToM();
      indirect[CL_TEST]=test;
      indirect[CL_STDFUN]=StdFun;      
  }
  return 0;
}

/***************************/
/* User-defined functions  */
Entry(Int2) mpd1_DefUser (Vector x, Vector y) {
Int2 i,j;
  for (i=j=0; i<nUTest; i++) 
    if (ActiveUTestFuns[i]==UDF_APPEND) {
      (stagenData->udFunc)(x,i,&y[j]);
      j++;
    }
  return 0;
}

/*******/
/* RHS */
Entry(Int2) mpd1_DefRHS (Vector x, Vector y) {
Int2 i;
  CopyToP(x);
  staFunc->CopyVectorElements(x,0,y,0,nphase);
  for (i=1; i<=_itmap; i++) {
    CopyToX(y);
    stagenData->Rhs(X,P,&zero,y);
  }
  return 0;
}

/**************************************/
/* Defining function for fixed points */
Local(Int2) DefFixed (Vector x, Vector y) {
Int2 i;
  mpd1_DefRHS(x,y);
  for (i=0; i<nphase; i++) y[i]-=x[i];
  if (nappend) mpd1_DefUser(x,y+nphase);
  return 0;
}

/********************/
/* Jacobian of RHS  */
Local(Int2) JacRHS (Vector x, Matrix Jac) {
  staFunc->CopyVectorElements(x,0,xp,0,nphase+2+nappend);     
  if (stagenData->Der1) {
    SymJac(xp,Jac);
  } else {
    if (staFunc->Der(mpd1_DefRHS, nphase, xp, _dx, Jac)) return 1;
  }
  return 0;
}

/****************************************************/
/* Jacobian of defining conditions for fixed points */
Local(Int2) JacDefFixed (Vector x, Matrix Jac) {
Int2 i;
  staFunc->ClearMatrix(Jac);
  JacRHS(x, Jac);
  for (i=0;i<nphase; i++) Jac[i][i]-=1.0;
  if (nappend) staFunc->Der(mpd1_DefUser, nappend, x, _dx, Jac+nphase);
  return 0;
}

/***************************/
/* Flip defining condition */
Entry(Int2) mpd1_DefBord (Vector x, Vector y) {
Int2 i,j,p,q,r,s;
Int2 Ret;
 
  if (JacRHS(x,C1)) return 1;	   
/* Fill in the bordered matrix BP */
  staFunc->ClearMatrix(BP);
  staFunc->CopyMatrixBlock(C1,0,0,BP,0,0,nphase,nphase);  
  for (i=0;i<nphase;i++) BP[i][i]+=1.0;

/* global V11 and W11 are assumed to be correct for bordering */
  for (i=0; i<nphase; i++) {
    BP[nphase][i]=V11[i];
  }
  for (i=0; i<nphase; i++) {
    BP[i][nphase]=W11[i];
  }

/* compute fold defining function */
  staFunc->CopyMatrix(BP,BT);
  Ret=staFunc->Decomp(BT);
  if (Ret != 0) {
    return 1;
  }
  staFunc->ClearVector(VB);
  VB[nphase]=1.0;
  staFunc->Solve(BT,VB);
  staFunc->CopyVectorElements(VB,0,Q,0,nphase);     /* global vector Q is computed here */
  staFunc->NormalizeVector(Q);
  y[0]=VB[nphase];
  return 0;
}

/***************************************/
/* Jacobian of fold defining condition */
Local(Int2) mpd1_JacDefBord (Vector x, Matrix Jac) {
int i,j,l;
VAR S;

  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2) {
/* creation of DefBord Jacobian matrix Jac using symbolic second derivatives */
    H=SymHess(x);
/* take global vector VB generated at the computation of mpd1_DefBord */
/* compute PP using the transposed bordered Jacobian matrix */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->Decomp(BT);
    staFunc->ClearVector(PP);
    PP[nphase]=1.0;
    staFunc->Solve(BT,PP);
/* compute gradient of mpd1_DefBord */
    for (l=0; l<nphase+2+nappend; l++) {
      for(i=0, S=0.0; i<nphase; i++) {
        for (j=0; j<nphase; j++) {
          S-=PP[i]*VB[j]*HessianElem(i,j,l);
        }
      }
      Jac[0][l]=S;
    }    
  } else {
    if (staFunc->Der(mpd1_DefBord, 1, x, _dx, Jac)) return 1;
  }
  return(0);
}

/*************************************/
/* Defining functions for Flip curve */

Entry(Int2) mpd1_DefFlip (Vector x, Vector y) {
  if (DefFixed(x,y)) return 1;
  if (mpd1_DefBord(x,y+nphase+nappend)) return 2;
  return 0;
}

/*************************************************/
/* Jacobian of defining functions for Flip curve */  
Entry(Int2) mpd1_JacDefFlip (Vector x, Matrix Jac) {
  int i;
  staFunc->ClearMatrix(Jac);
  if (JacDefFixed(x,C)) return 1;
  staFunc->CopyMatrixBlock(C,0,0,Jac,0,0,nphase+nappend,nphase+2+nappend);
  if (mpd1_JacDefBord(x,C2)) return 2;
  staFunc->CopyMatrixBlock(C2,0,0,Jac,nphase+nappend,0,1,nphase+2+nappend);
  return 0;
}

/*********************/
/* Default processor */
#define SCALE	45.0/atan(1.0)

#ifdef WIN32
static int _USERENTRY fp_SortMult(const void *f, const void *s) {
#else
Entry(int) fp_SortMult(const void *f, const void *s) {
#endif
VAR Mod1,Mod2;
  Mod1=((FloatHiPtr)f)[0]*((FloatHiPtr)f)[0]+((FloatHiPtr)f)[1]*((FloatHiPtr)f)[1];
  Mod2=((FloatHiPtr)s)[0]*((FloatHiPtr)s)[0]+((FloatHiPtr)s)[1]*((FloatHiPtr)s)[1];
  if (Mod1==Mod2) {
    return 0;
  }
  return (Mod1<Mod2) ? 1 : -1;
}

Entry(Int2) mpd1_ProcDefault(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {
VAR s;
Int2 i;

  if (Ind == -1) {  /* zero of user-defined function */
    BifVector[0]=(VAR)_itmap;
    stagenData->BifDataOutLen=sizeof(VAR);     
    return 0;
  }
  /* Assert: v1==NULL, msg==NULL */
  /* mpd1_DefFlip (x, xp1);  compute global Q for test functions */
  MemCopy(x1,x,sizeof(FloatHi)*(nphase+2+nappend));     /* for the last point */

  /* L2_NORM and Min&Max */
  for (i=0,s=0; i<nphase; i++) {
    s+=x[i]*x[i];
    StdFun[1+i]=x[i];
    StdFun[1+nphase+i]=x[i];
  }
  StdFun[0]=sqrt(s);
  /* Multipliers */
  if (staData->eigcomp) {
    if (JacRHS(x, A0)) return (1);     
    staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase); 
    staFunc->EigenValues(AA,Re,Im);
    for (i=0; i<nphase; i++) {
      M[2*i]=Re[i];
      M[2*i+1]=Im[i];
    }
    qsort(M,nphase,2*sizeof(VAR),fp_SortMult);
    for (i=0; i<nphase; i++) {
      Re[i]=M[2*i];
      Im[i]=M[2*i+1];
    }
    for (i=0; i<nphase; i++) {
      Mod[i]=sqrt(Re[i]*Re[i]+Im[i]*Im[i]);
      if (Mod[i]==0.0) {
        Arg[i]=0.0;
      } else {
        if (fabs(Re[i]/Mod[i]) > 0.5) {
          Arg[i]=atan(fabs(Im[i]/Re[i]));
        } else {
          if (Re[i]!=0.0) {
            Arg[i]=PID2-atan(fabs(Re[i]/Im[i]));
	  } else {
	    Arg[i]=PID2;
	  }
        }
      }
      if (Re[i] < 0) Arg[i]=PI-Arg[i];
      if (Im[i] < 0) Arg[i]=-Arg[i];
      Arg[i]*=SCALE;
    }
  }
  return 0;
}

/*********************************/
/* Test and processing functions */
typedef EntryPtr(Int2,TestFuncPtr,(Vector xx, Vector vv, FloatHiPtr TV));
typedef EntryPtr(Int2,ProcFuncPtr,(Int2 Ind, Vector xx, Vector vv, CharPtr *msg));

Local(Char) buf[120];

/* cubic coefficient */
Entry(Int2) mpd1_TestDegen(Vector x, Vector v, FloatHiPtr res) {
Int2 i,j,k,l,Ret=0;
VAR b,c;

/* it is assumed that global Q is correct */   

/*  extract coordinates, parameters, and the flip eigenvector V */
    staFunc->CopyVectorElements(x,0,xp,0,nphase+2+nappend);
    staFunc->CopyVectorElements(Q,0,V,0,nphase);
    staFunc->ClearVector(Vp);
    staFunc->CopyVectorElements(Q,0,Vp,0,nphase);
    
/*  compute adjoint eigenvector W */
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) {
      return 1;
    }
    staFunc->ClearVector(PP);
    PP[nphase]=1.0;
    staFunc->Solve(BT,PP);
    staFunc->CopyVectorElements(PP,0,W,0,nphase);

    staFunc->NormalizeVectorRelative(W,V);

/* compute the coefficient in the normal form x'=-x+cx^3 */
 
    staFunc->ClearVector(v0);
    staFunc->CopyVectorElements(V,0,v01,0,nphase);
    staFunc->CopyVector(x,x01);

/* compute vector a=B(V,V) */
    if (stagenData->Der2) {
        H=SymHess(x);
        for (i=0; i<nphase; i++) {
           for (a[i]=0.0,j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) a[i]+=HessianElem(i,j,k)*V[j]*V[k];    /* a=B(V,V) */
           }
        }
    } else {
        staFunc->AddScaledVector(x01,v01,ddx,x11);
        mpd1_DefRHS(x11,a);
        mpd1_DefRHS(x01,F0);
	staFunc->AddScaledVector(a,F0,-2.0,a);
        staFunc->AddScaledVector(x01,v01,-ddx,x11);
        mpd1_DefRHS(x11,F0);
	staFunc->AddScaledVector(a,F0,1.0,a);
	staFunc->MultiplyVectorScalar(a,1.0/ddx/ddx,a);
    };

/* solve the system for a */
/*  compose Jacobian matrix */
    JacRHS(x,C1);
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->CopyMatrix(AA,B);
    for (i=0; i<nphase; i++) B[i][i]-=1.0;
    Ret=staFunc->Decomp(B);
    if (Ret) {
      *res=1.0;
      goto final;
    }
    staFunc->Solve(B,a);

/* compute F0=B(V,a) and  b=<W,F0> */
    if (stagenData->Der2) {
      for (i=0; i<nphase; i++) {
        for (F0[i]=0.0,j=0; j<nphase; j++) {
          for (k=0; k<nphase; k++) F0[i]+=HessianElem(i,j,k)*V[j]*a[k];    
        }
      }
    } else {
      staFunc->AddScaledVector(V,a,1.0,F1);
      staFunc->CopyVectorElements(F1,0,v01,0,nphase);
      staFunc->AddScaledVector(x01,v01,0.5*ddx,x11);
      mpd1_DefRHS(x11,F1);
      staFunc->AddScaledVector(x01,v01,-0.5*ddx,x11);
      mpd1_DefRHS(x11,F0);
      staFunc->AddScaledVector(F1,F0,1.0,F2);

      staFunc->AddScaledVector(V,a,-1.0,F1);
      staFunc->CopyVectorElements(F1,0,v01,0,nphase);
      staFunc->AddScaledVector(x01,v01,0.5*ddx,x11);
      mpd1_DefRHS(x11,F1);
      staFunc->AddScaledVector(x01,v01,-0.5*ddx,x11);
      mpd1_DefRHS(x11,F0);
      staFunc->AddScaledVector(F1,F0,1.0,F1);
      staFunc->AddScaledVector(F2,F1,-1.0,F0);
      staFunc->MultiplyVectorScalar(F0,1.0/ddx/ddx,F0);   
    }
    b=staFunc->ScalProd(W,F0); 
/* compute c */
    if (stagenData->Der3) {
      D3=SymD3(x);
      for (c=0.0,i=0; i<nphase; i++) {
         for (j=0; j<nphase; j++) {
            for (k=0; k<nphase; k++) {
	       for (l=0; l<nphase; l++) {
	         c+=W[i]*D3Elem(i,j,k,l)*V[j]*V[k]*V[l];
	       }
	    }
         }
      }
    } else {
      staFunc->CopyVectorElements(V,0,v01,0,nphase);
      staFunc->AddScaledVector(x01,v01,3*dddx/2,x11);
      mpd1_DefRHS(x11,F1);
      c=staFunc->ScalProd(W,F1);
      staFunc->AddScaledVector(x01,v01,dddx/2,x11);
      mpd1_DefRHS(x11,F1);
      c-=3.0*staFunc->ScalProd(W,F1);
      staFunc->AddScaledVector(x01,v01,-dddx/2,x11);
      mpd1_DefRHS(x11,F1);
      c+=3.0*staFunc->ScalProd(W,F1);
      staFunc->AddScaledVector(x01,v01,-3*dddx/2,x11);
      mpd1_DefRHS(x11,F1);
      c-=staFunc->ScalProd(W,F1);
      c/=(dddx*dddx*dddx);   
    }
    *res=-b/2+c/6;
final:    
  test[0]=*res;
  return Ret;
}


Entry(Int2) mpd1_ProcDegen(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
Int2 i;
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Degenerate flip (?)");
  else {
/*  provide bifurcation data for DF ->* starters */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,Q,sizeof(FloatHi)*nphase);
    stagenData->BifDataOutLen=sizeof(FloatHi)*(nphase+1);
    sprintf(buf,"Degenerate flip");
  }
  *msg=buf;
  return 0;
}


/* product <W,V> */
Entry(Int2) mpd1_TestRes12(Vector x, Vector v, FloatHiPtr res) {
Int2 i;

  int Ret=0;
  *res=1.0;
  if (nphase < 2) goto end;

/*  compose Jacobian matrix */
  staFunc->CopyVectorElements(x,0,xp,0,nphase+2+nappend);
  staFunc->CopyVectorElements(Q,0,V,0,nphase);
  if (JacRHS(x,C1)) {
    Ret=1;
    goto end;
  }
  staFunc->CopyMatrixBlock(C1,0,0,BT,0,0,nphase,nphase);
  for (i=0;i<nphase;i++) {
    BT[nphase][i]=V11[i]; 
    BT[i][nphase]=W11[i];
    BT[i][i]+=1.0;
  }
  
  Ret=staFunc->Decomp(BT);
  if (Ret) goto end;
  
  staFunc->ClearVector(VB);
  staFunc->CopyVectorElements(Q,0,VB,0,nphase);
  staFunc->Solve(BT,VB);
  *res=VB[nphase];
  
end:
  test[1]=*res;
  return Ret;
}

Entry(Int2) mpd1_ProcRes12(Int2 Ind, Vector x, Vector v, CharPtr *msg) {

  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"1:2 Resonance (?)");
  else {
Vector v1,W0,W1;
Matrix D;
Int2 i;
double d,u;
    v1=staFunc->CreateVector(2*nphase);
    W0=staFunc->CreateVector(nphase);
    W1=staFunc->CreateVector(nphase);
    D=staFunc->CreateMatrix(2*nphase,2*nphase);

/*  compute (generalized) null-eigenvectors  */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0;
    
    staFunc->CopyMatrixBlock(A,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(A,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"1:2 Resonance : eigenvectors(?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,V,0,nphase);       
    d=sqrt(staFunc->ScalProd(V,V));
    staFunc->NormalizeVector(V);
    staFunc->CopyVectorElements(Q,0,W,0,nphase);
    u=staFunc->ScalProd(V,W);
    staFunc->MultiplyVectorScalar(V,u,V);	/* null-vector as continued */
    staFunc->CopyVectorElements(v1,nphase,W,0,nphase); 
    staFunc->MultiplyVectorScalar(W,u/d,W);	/* generalized null-vector */
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W); 

/*  provide bifurcation data for R2 -> * starters */
/*  BifVector[0]:                        _itmap */
/*  BifVector[1] to [nphase]:            null-vector V (A*V=0, <V,V>=1)*/
/*  BifVector[nphase+1] to [2*nphase]:   generalized null-vector W (A*W=V, <W,V>=0) */
/*  BifVector[2*nphase+1]:               cos(theta)  */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,W,sizeof(FloatHi)*nphase);

    BifVector[2*nphase+1]=-1.0;
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+2);

/* compute the adjoint eigenvectors */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0;
    staFunc->TransposeMatrix(A,B);
    staFunc->ClearMatrix(D);
    staFunc->CopyMatrixBlock(B,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(B,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"1:2 Resonance : adjoint eigenvectors (?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,W1,0,nphase);       /* adjoint null-vector W1 */
    staFunc->CopyVectorElements(v1,nphase,W0,0,nphase);  /* adjoint generalized null-vector W0 */
    d=staFunc->ScalProd(V,W0);
    staFunc->MultiplyVectorScalar(W1,1.0/d,W1);
    staFunc->AddScaledVector(W0,W1,-staFunc->ScalProd(W0,W),W0);
    staFunc->MultiplyVectorScalar(W0,1.0/d,W0);

    sprintf(buf,"1:2 Resonance");

final:
    D=staFunc->DestroyMatrix(D);
    v1=staFunc->DestroyVector(v1);
    W0=staFunc->DestroyVector(W0);
    W1=staFunc->DestroyVector(W1);
  }
  *msg=buf;
  return 0;
}

/* unit-product test function */
Entry(Int2) mpd1_TestFlipNS(Vector x, Vector v, FloatHiPtr res) {

  int i,j,p,q,r,s,m;

  *res=1.0;
  if (nphase < 2) goto end;

  if (JacRHS(x,C1)) return 1;
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

  m=nphase*(nphase-1)/2;
  staFunc->ClearMatrix(BIP);
  i=-1;
  for (p=1; p<nphase; p++) {
    for (q=0; q<p; q++) {
      i++;
      j=-1;
      for (r=1; r<nphase; r++) {
        for (s=0; s<r; s++) {
          j++;
          BIP[i][j]=A[p][r]*A[q][s]-A[q][r]*A[p][s];
        }
      }
    }
  }
  
  for (i=0;i<m;i++) BIP[i][i]-=1.0;
  
  if (staFunc->Decomp(BIP)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(BIP,res,NULL);
  }
end:
  test[2]=*res;
  return 0;
}

Entry(Int2) mpd1_ProcFlipNS(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;

 if (Ind) {
    sprintf(buf,"Flip + Neimark-Sacker or neutral saddle (?)");
  } else {
/*  provide bifurcation data for FN->* starters */
/*  BifVector[0]: 			 _itmap */
/*  BifVector[1] to [nphase]:          flip vector  */
/*  BifVector[nphase+1] to [3*nphase]: NS vectors */
/*  BifVector[3*nphase+1]:             kappa */
    BifVector[0]=(VAR)_itmap;
{
Int2 i,j,k,l,imin,jmin,Ret;
VAR d,s,s1,r,r1,h,beta,gamma,phi,lambda,lambda1,omega,lambda2,omega2,theta;
Matrix AT,BB;
Vector V,Qr,Qi,Pr,Pi;
   
    AT=staFunc->CreateMatrix(nphase,nphase);
    BB=staFunc->CreateMatrix(2*nphase,2*nphase);
    V=staFunc->CreateVector(2*nphase);
    Qr=staFunc->CreateVector(nphase);
    Qi=staFunc->CreateVector(nphase);
    Pr=staFunc->CreateVector(nphase);
    Pi=staFunc->CreateVector(nphase);

/*  compose Jacobian matrix */
    JacRHS(x,C1);
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    for (i=0; i<nphase; i++) {
      M[2*i]=Re[i];
      M[2*i+1]=Im[i]; 
    }
    qsort(M,nphase,2*sizeof(VAR),fp_SortMult);
    for (i=0; i<nphase; i++) {
      Re[i]=M[2*i];
      Im[i]=M[2*i+1];
    }

    /* find the pair of multipliers on the unit circle */
    s=fabs(Re[0]*Re[1]-Im[0]*Im[1]-1.0)+fabs(Re[0]*Im[1]+Re[1]*Im[0]);
    imin=0;
    jmin=1;
    for (i=0; i<nphase; i++) {
      for (j=i+1; j<nphase; j++) {
        d=fabs(Re[i]*Re[j]-Im[i]*Im[j]-1.0)+fabs(Re[i]*Im[j]+Re[j]*Im[i]);
        if (d<s) {
          imin=i;
          jmin=j;
          s=d;
        }
      }
    }
    if (fabs(Im[imin])+fabs(Im[jmin])==0) {
/* Flip + Neutral saddle */
      lambda=Re[imin];
      lambda1=Re[jmin];
      /* compute real eigenvectors associated with lambda and lambda1=1/lambda */
      staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=lambda;
      if (staFunc->SngSolve(AA,_eps,Qr) != nphase-1) {
        sprintf(buf,"Flip + neutral saddle (?) (cannot find ns-eigenvector)");
        goto final;
      } 
      staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=lambda1;
      if (staFunc->SngSolve(AA,_eps,Qi) != nphase-1) {
        sprintf(buf,"Flip + neutral saddle (?) (cannot find ns-eigenvector)");
        goto final;
      } 
      staFunc->NormalizeVector(Qr);   
      staFunc->NormalizeVectorRelative(Qi,Qr);    /* ensures smallest possible angle */
      staFunc->NormalizeVector(Qi);

      /* provide bifurcation data for FN->* starter */
      staFunc->AddScaledVector(Qr,Qi,1.0,Pr);
      staFunc->NormalizeVector(Pr);
      MemCopy(BifVector+nphase+1,Pr,sizeof(FloatHi)*nphase);
      staFunc->AddScaledVector(Qr,Qi,-1.0,Pr);
      staFunc->NormalizeVector(Pr);
      MemCopy(BifVector+2*nphase+1,Pr,sizeof(FloatHi)*nphase);
      BifVector[3*nphase+1]=(lambda+lambda1)/2.0;  /* kappa */
      stagenData->BifDataOutLen=sizeof(FloatHi)*(3*nphase+1+1);      
      sprintf(buf,"Flip + neutral saddle: lambda=%g",lambda);
   } else {
/* Flip + Neimark-Sacker */
      if (fabs(Re[imin]) > 0.5) {
        theta=atan(fabs(Im[imin]/Re[imin]));
      } else {
        theta=PID2-atan(fabs(Re[imin]/Im[imin]));
      }
      if (Re[imin] < 0) theta=PI-theta;
      if (Im[imin] < 0) theta=-theta;
      lambda=cos(theta);
      omega=sin(theta);
      lambda2=cos(2*theta);
      omega2=sin(2*theta); 

      /* compute real and imaginary part of the eigenvector */
      staFunc->CopyMatrixBlock(C1,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(C1,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) {
        BB[i][i+nphase]=omega;
        BB[i][i]-=lambda;
        BB[i+nphase][i]=-omega;
        BB[i+nphase][i+nphase]-=lambda;
      }
      if (staFunc->SngSolve(BB,_eps,VV) != 2*nphase-2) {
        sprintf(buf,"Flip + Neimark-Sacker (?) (singular BB)");
        goto final;
      }
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);
      /* normalize real and imaginary parts <Qr,Qr>=1, <Qr,Qi>=0 */ 
      d=staFunc->ScalProd(Qr,Qr);
      s=staFunc->ScalProd(Qi,Qi);
      r=staFunc->ScalProd(Qr,Qi);
      if (r!=0) {
           beta=sqrt(fabs(4.0*r*r+(s-d)*(s-d)));
           if (2.0*fabs(2.0*r) > beta) {
             phi=PID2-atan(fabs((s-d)/(2.0*r)));
           } else {
             phi=atan(fabs(2.0*r/(s-d)));
           }
           if (r < 0) phi=PI-phi;
           if ((s-d) < 0) phi=-phi;
           phi/=2.0;
           for (i=0; i<nphase; i++) {
              VV[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
              VV[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
           }
      }
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      d=staFunc->ScalProd(Qr,Qr);
      gamma=1/sqrt(d);
      staFunc->MultiplyVectorScalar(VV,gamma,VV);
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);
      /* provision of bif data */
      MemCopy(BifVector+nphase+1,VV,sizeof(FloatHi)*2*nphase);
      BifVector[3*nphase+1]=cos(theta);     /* kappa */
      stagenData->BifDataOutLen=sizeof(FloatHi)*(3*nphase+2); 
      sprintf(buf,"Flip + Neimark-Sacker: theta=%g",theta);
    }
final:;
  AT=staFunc->DestroyMatrix(AT);
  BB=staFunc->DestroyMatrix(BB);
  V=staFunc->DestroyVector(V);
  Qr=staFunc->DestroyVector(Qr);
  Qi=staFunc->DestroyVector(Qi);
  Pr=staFunc->DestroyVector(Pr);
  Pi=staFunc->DestroyVector(Pi);
}
  } 
  *msg=buf;
  return 0;
}

Entry(Int2) mpd1_TestFlipFold(Vector x, Vector v, FloatHiPtr res) {
Int2 i;

  int Ret=0;
  *res=1;
  if (nphase < 2) goto end;

/* compose Jacobian matrix */
  staFunc->CopyVectorElements(x,0,xp,0,nphase+2+nappend);
  staFunc->CopyVectorElements(Q,0,V,0,nphase);
  if (JacRHS(x,C1)) {
    Ret=1;
    goto end;
  }
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

/* det(A-I) */
  for (i=0;i<nphase;i++) A[i][i]-=1.0;
  if (staFunc->Decomp(A)) {
    *res=0.0;
    Ret=1.0;
  } else {
      staFunc->GetMatrixData(A,res,NULL);
  }      
end:
  test[3]=*res;
  return Ret;
}

Entry(Int2) mpd1_ProcFlipFold(Int2 Ind, Vector x, Vector v, CharPtr *msg) {

  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Fold + flip (?)");
  else {
Vector v1,W0,W1;
Matrix D;
Int2 i;
double d,u;
    v1=staFunc->CreateVector(2*nphase);
    W0=staFunc->CreateVector(nphase);
    W1=staFunc->CreateVector(nphase);
    D=staFunc->CreateMatrix(2*nphase,2*nphase);

/*  compute (generalized) null-eigenvectors  */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    
    staFunc->CopyMatrixBlock(A,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(A,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"Fold + flip: eigenvectors(?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,V,0,nphase);       
    d=sqrt(staFunc->ScalProd(V,V));
    staFunc->NormalizeVector(V);
    staFunc->CopyVectorElements(Q,0,W,0,nphase);
    u=staFunc->ScalProd(V,W);
    staFunc->MultiplyVectorScalar(V,u,V);	/* null-vector as continued */
    staFunc->CopyVectorElements(v1,nphase,W,0,nphase); 
    staFunc->MultiplyVectorScalar(W,u/d,W);	/* generalized null-vector */
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W); 

/*  provide bifurcation data for bt->* starters */
/*  BifVector[0] to [nphase-1]:          null-vector V (A*V=0, <V,V>=1)*/
/*  BifVector[nphase] to [2*nphase-1]:   generalized null-vector W (A*W=V, <W,V>=0) */
/*  BifVector[2*nphase]:                 -lambda*lambda  */
    MemCopy(BifVector,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase,W,sizeof(FloatHi)*nphase);

    BifVector[2*nphase]=0.0;
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1);

/* compute the adjoint eigenvectors */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    staFunc->TransposeMatrix(A,B);
    staFunc->ClearMatrix(D);
    staFunc->CopyMatrixBlock(B,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(B,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"Fold + flip: adjoint eigenvectors (?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,W1,0,nphase);       /* adjoint null-vector W1 */
    staFunc->CopyVectorElements(v1,nphase,W0,0,nphase);  /* adjoint generalized null-vector W0 */
    d=staFunc->ScalProd(V,W0);
    staFunc->MultiplyVectorScalar(W1,1.0/d,W1);
    staFunc->AddScaledVector(W0,W1,-staFunc->ScalProd(W0,W),W0);
    staFunc->MultiplyVectorScalar(W0,1.0/d,W0);

    sprintf(buf,"Fold - flip");

final:
    D=staFunc->DestroyMatrix(D);
    v1=staFunc->DestroyVector(v1);
    W0=staFunc->DestroyVector(W0);
    W1=staFunc->DestroyVector(W1);
  }
  *msg=buf;
  return 0;
}

/****************************************************/
/*  Adaptation of the bordering vectors V11 and W11 */
Entry(Int2) mpd1_Adapter(Vector x, Vector v) {
int Ret;  

/* global normalized Q assumed to be correct */
  staFunc->CopyVector(Q,V11);

  staFunc->TransposeMatrix(BP,BT);
  Ret=staFunc->Decomp(BT);
  if (Ret != 0) {
    return 1;
  }
  staFunc->ClearVector(PP);
  PP[nphase]=1.0;
  staFunc->Solve(BT,PP);
  staFunc->CopyVectorElements(PP,0,W11,0,nphase);
  staFunc->NormalizeVector(W11);
  return 0;
}