/*
   __zh
   Zero-Hopf Curve
   September 30, 1997
*/

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

#include "stagen.h"

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

#include "__zh.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_EIGV		4
#define CL_TEST		5
#define CL_UTEST	6

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

#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(__zh_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) E=NULL;		/* current set of eigenvalues */
Local(VAR) zero=0;

#include "_symdiff.c"

/* Working space for defining and test functions  */

Local(Matrix) A,B,C,C1,C2,AA,BB,BP,BT,BP1,BT1,BP2,BT2,D,DE,G;
Local(Vector) v1,vn,xp,xp1,V,W,F,Vp,L,VV,VB11,VB21,VB12,VB22,W0,Q,PP11,PP21,PP12,PP22;
Local(Vector) TV,K1,K2,K3,OK1,OK2,V11,W11,V12,V22,W12,W22,PIVOT,RL,OLD;
Local(Vector) Re=NULL,Im=NULL;
Local(VAR) test[2],dx,ddx;

#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=(__zh_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];
  /* 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+3));    
  if (staUtil->CheckParameters(staData->Active,npar,active,3+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 (!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;
  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 */
    E=MemNew(sizeof(FloatHi)*2*nphase);
    Re=staFunc->CreateVector(nphase);
    Im=staFunc->CreateVector(nphase);
    /* Allocate memory for working space used in defining and test functions */
    A=staFunc->CreateMatrix(nphase,nphase);
    B=staFunc->CreateMatrix(nphase,nphase);
    BP=staFunc->CreateMatrix(nphase+1,nphase+1);
    BT=staFunc->CreateMatrix(nphase+1,nphase+1);
    W0=staFunc->CreateVector(nphase);
    VB11=staFunc->CreateVector(nphase+1);
    VB21=staFunc->CreateVector(nphase+1);
    VB12=staFunc->CreateVector(m+2);
    VB22=staFunc->CreateVector(m+2);
    VV=staFunc->CreateVector(2*nphase);
    Q=staFunc->CreateVector(nphase);
    PP11=staFunc->CreateVector(nphase+1);
    PP21=staFunc->CreateVector(nphase+1);
    PP12=staFunc->CreateVector(m+2);
    PP22=staFunc->CreateVector(m+2);
    TV=staFunc->CreateVector(nphase+3);
        
    AA=staFunc->CreateMatrix(nphase,nphase);
    BB=staFunc->CreateMatrix(2*nphase,2*nphase);
    C=staFunc->CreateMatrix(nphase+nappend,nphase+3+nappend);
    C1=staFunc->CreateMatrix(nphase,nphase+3+nappend);
    C2=staFunc->CreateMatrix(2,nphase+3+nappend);
    x0=staFunc->CreateVector(nphase+3+nappend);
    x1=staFunc->CreateVector(nphase+3+nappend);
    v0=staFunc->CreateVector(nphase+3+nappend);
    v1=staFunc->CreateVector(2*nphase);
    vn=staFunc->CreateVector(nphase);
    xp=staFunc->CreateVector(nphase+3+nappend);
    xp1=staFunc->CreateVector(nphase+3+nappend);
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    F=staFunc->CreateVector(nphase);
    L=staFunc->CreateVector(nphase);
    Vp=staFunc->CreateVector(nphase+3+nappend);
    BifVector=MemNew(sizeof(FloatHi)*(4*nphase+2));
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);
    BP1=staFunc->CreateMatrix(m+1,m+1);
    BT1=staFunc->CreateMatrix(m+1,m+1);    
    BP2=staFunc->CreateMatrix(m+2,m+2);
    BT2=staFunc->CreateMatrix(m+2,m+2);
    G=staFunc->CreateMatrix(m,m);
    K1=staFunc->CreateVector(m+2);
    K2=staFunc->CreateVector(m+2);          
    K3=staFunc->CreateVector(m+2);
    OK1=staFunc->CreateVector(m+2);
    OK2=staFunc->CreateVector(m+2);
    V11=staFunc->CreateVector(nphase);
    W11=staFunc->CreateVector(nphase);
    V12=staFunc->CreateVector(m+2);
    V22=staFunc->CreateVector(m+2);
    W12=staFunc->CreateVector(m+2);
    W22=staFunc->CreateVector(m+2);
    PIVOT=staFunc->CreateVector(m);
    RL=staFunc->CreateVector(m+1);
    OLD=staFunc->CreateVector(m+1);    
  }
  Der2num=ind2_(nphase+npar-1,nphase+npar-1)+1;
  Der3num=ind3_(nphase+npar-1,nphase+npar-1,nphase+npar-1)+1;
  dx=_dx;
  ddx=pow(dx,3.0/4.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;
  m=nphase*(nphase-1)/2;
  stagenData->BifDataOutLen=0;
/* saving global vectors V11,W11,V12,V22,W12,W22 used in defining functions */
    if (OK&&P&&(stagenData->Reason==RF_CLEAN)) {
      stagenData->ContDataStaLen=(2*nphase+4*(m+2))*sizeof(FloatHi);
      stagenData->ContDataStaPtr=MemNew(stagenData->ContDataStaLen);
      MemCopy(stagenData->ContDataStaPtr,V11,(nphase)*sizeof(FloatHi));
      MemCopy(stagenData->ContDataStaPtr+(nphase),W11,(nphase)*sizeof(FloatHi));
      MemCopy(stagenData->ContDataStaPtr+2*nphase,V12,(m+2)*sizeof(FloatHi));
      MemCopy(stagenData->ContDataStaPtr+2*nphase+(m+2),V22,(m+2)*sizeof(FloatHi));
      MemCopy(stagenData->ContDataStaPtr+2*(nphase+(m+2)),W12,(m+2)*sizeof(FloatHi));
      MemCopy(stagenData->ContDataStaPtr+2*nphase+3*(m+2),W22,(m+2)*sizeof(FloatHi));
      stagenData->BifDataOutLen=0;
      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);
      E=MemFree(E);
      Re=staFunc->DestroyVector(Re);
      Im=staFunc->DestroyVector(Im);
      A=staFunc->DestroyMatrix(A);
      B=staFunc->DestroyMatrix(B);
      BP=staFunc->DestroyMatrix(BP);
      BT=staFunc->DestroyMatrix(BT);
      W0=staFunc->DestroyVector(W0);
      Q=staFunc->DestroyVector(Q);
      PP11=staFunc->DestroyVector(PP11);
      PP21=staFunc->DestroyVector(PP21);
      PP12=staFunc->DestroyVector(PP12);
      PP22=staFunc->DestroyVector(PP22);
      TV=staFunc->DestroyVector(TV);
      AA=staFunc->DestroyMatrix(AA);
      BB=staFunc->DestroyMatrix(BB);
      VV=staFunc->DestroyVector(VV);
      VB11=staFunc->DestroyVector(VB11);
      VB21=staFunc->DestroyVector(VB21);
      VB12=staFunc->DestroyVector(VB12);
      VB22=staFunc->DestroyVector(VB22);
      C=staFunc->DestroyMatrix(C);
      C1=staFunc->DestroyMatrix(C1);
      C2=staFunc->DestroyMatrix(C2);
      vn=staFunc->DestroyVector(vn);
      xp=staFunc->DestroyVector(xp);
      xp1=staFunc->DestroyVector(xp1);
      v1=staFunc->DestroyVector(v1);
      V=staFunc->DestroyVector(V);
      W=staFunc->DestroyVector(W);
      F=staFunc->DestroyVector(F);
      Vp=staFunc->DestroyVector(Vp);
      L=staFunc->DestroyVector(L);
      BifVector=MemFree(BifVector);
      D=staFunc->DestroyMatrix(D);
      DE=staFunc->DestroyMatrix(DE);
      BP1=staFunc->DestroyMatrix(BP1);
      BT1=staFunc->DestroyMatrix(BT1);
      BP2=staFunc->DestroyMatrix(BP2);
      BT2=staFunc->DestroyMatrix(BT2);
      G=staFunc->DestroyMatrix(G);
      K1=staFunc->DestroyVector(K1);
      K2=staFunc->DestroyVector(K2);     
      K3=staFunc->DestroyVector(K3);
      OK1=staFunc->DestroyVector(OK1);
      OK2=staFunc->DestroyVector(OK2);
      V11=staFunc->DestroyVector(V11);
      W11=staFunc->DestroyVector(W11);
      V12=staFunc->DestroyVector(V12);
      V22=staFunc->DestroyVector(V22);
      W12=staFunc->DestroyVector(W12);
      W22=staFunc->DestroyVector(W22);
      PIVOT=staFunc->DestroyVector(PIVOT);
      RL=staFunc->DestroyVector(RL);
      OLD=staFunc->DestroyVector(OLD);      
    }
} 

Entry(Int2) zh_JacDefZeroHopf (Vector x, Matrix J);
Entry(Int2) zh_DefZeroHopf (Vector x, Vector J);


/************/
/* Starters */
/* zh -> zh Starter */
Entry(Int2) zh_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,Ret,m,p,q,r,s;
  Vector v1,v2;
  Matrix D,DE;
  VAR max1,max2, alpha, beta, 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 (stagenptr->ContDataGenLen) {       /* "Extention" */
    MemCopy(V11,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase));
    MemCopy(W11,stagenptr->ContDataStaPtr+(nphase),sizeof(FloatHi)*(nphase));
    MemCopy(V12,stagenptr->ContDataStaPtr+(2*nphase),sizeof(FloatHi)*(m+2));
    MemCopy(V22,stagenptr->ContDataStaPtr+(2*nphase+(m+2)),sizeof(FloatHi)*(m+2));
    MemCopy(W12,stagenptr->ContDataStaPtr+2*(nphase+(m+2)),sizeof(FloatHi)*(m+2));
    MemCopy(W22,stagenptr->ContDataStaPtr+(2*nphase+3*(m+2)),sizeof(FloatHi)*(m+2));
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<3+nappend; i++) x0[nphase+i]=staData->P[active[i]];    
  } else {
    v1=staFunc->CreateVector(nphase+3+nappend);
    v2=staFunc->CreateVector(nphase+3+nappend);
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<3+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
 
    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];
/* use bifurcation data if available; otherwise use pivot information */
/* BifVector[0] to [nphase-1]:       null-vector V (A*V=0, <V,V>=1) */
    if (stagenData->BifDataInOk) {
      MemCopy(V11,stagenData->BifDataInPtr,sizeof(FloatHi)*nphase);
      staFunc->NormalizeVector(V11); 
    } else {
        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[nphase][i]=V11[i];
      BP[i][nphase]=W11[i];
    }
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) { 
       goto final;
    }    
    staFunc->ClearVector(VB11); 
    VB11[nphase]=1.0; 
    staFunc->Solve(BT,VB11);
    staFunc->CopyVectorElements(VB11,0,V11,0,nphase);   
    staFunc->NormalizeVector(V11);    
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret != 0) { 
      goto final;
    }
    staFunc->ClearVector(PP11);
    PP11[nphase]=1.0;
    staFunc->Solve(BT,PP11);
    staFunc->CopyVectorElements(PP11,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 bordering of the biproduct */  
   staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
   staFunc->ClearMatrix(BP2);
   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++;
       if (r==q) BP2[i][j]=-A[p][s];
       if (r!=p && s==q) BP2[i][j]=A[p][r];
       if (r==p && s==q) BP2[i][j]=A[p][p]+A[q][q];
       if (r==p && s!=q) BP2[i][j]=A[q][s];
       if (s==p) BP2[i][j]=-A[q][r];
      }
     }
    }
   } 
   BP2[m][m-2]=1.0;
   BP2[m+1][m-1]=1.0;
   staFunc->ClearVector(V12); 
   staFunc->ClearVector(V22);
   V12[m-2]=1.0;
   V22[m-1]=1.0;
/* searching two rows with maximal pivots */
   staFunc->ClearMatrix(G);
   staFunc->CopyMatrixBlock(BP2,0,0,G,0,0,m,m);
   staFunc->ClearVector(PIVOT);
   PIVOT[m-1]=1.0;  
   for (k=0;k<=m-2;k++) {
    kp1=k+1;
    p=k;
    for (i=kp1;i<=m-1;i++) { 
     if (fabs(G[i][k]) > fabs(G[p][k]))
      p=i;
    }
    PIVOT[k]=p;
    if (p!=k) PIVOT[m-1]=-PIVOT[m-1];
    t=G[p][k];
    G[p][k]=G[k][k];
    G[k][k]=t; 
    if (t==0.0) continue;
    for (i=kp1;i<=m-1;i++) G[i][k]=-G[i][k]/t;
    for (j=kp1;j<=m-1;j++) {
     t=G[p][j];
     G[p][j]=G[k][j];
     G[k][j]=t;
     if (t != 0.0) for (i=kp1;i<=m-1;i++) G[i][j] += G[i][k]*t;
    }
   }
   staFunc->ClearVector(W12);
   for (i=0;i<m;i++) W12[i]=i;
   for (k=0;k<=m-2;k++) {
    i=(Int2)PIVOT[k];
    t=W12[i];
    W12[i]=W12[k];
    W12[k]=t;
   }
   p=(Int2)W12[m-2];
   q=(Int2)W12[m-1];
   BP2[p][m]=1.0;
   BP2[q][m+1]=1.0;
   staFunc->ClearVector(W12);
   staFunc->ClearVector(W22);
   W12[p]=1.0;
   W22[q]=1.0;
/* second bordering to prepare for g22 as test function H-BT */
   staFunc->ClearVector(K1);
   staFunc->ClearVector(K2);
   staFunc->CopyMatrix(BP2,BT2);
   K1[m]=1.0;
   K2[m+1]=1.0;
   staFunc->Decomp(BT2);
   staFunc->Solve(BT2,K1);
   staFunc->Solve(BT2,K2);
   staFunc->CopyVector(K1,OK1);
   staFunc->CopyVector(K1,K3);
   staFunc->CopyVector(K2,OK2);
   max1=(fabs(OK1[m])>fabs(OK2[m]))?fabs(OK1[m]):fabs(OK2[m]);
   max2=(fabs(OK1[m+1])>fabs(OK2[m+1]))?fabs(OK1[m+1]):fabs(OK2[m+1]);
   if (max1>max2) {
    alpha=OK1[m];
    beta=OK2[m];
   }
   else {
    alpha=OK1[m+1];
    beta=OK2[m+1];
   }
   staFunc->MultiplyVectorScalar(OK1,-beta,OK1);
   staFunc->AddScaledVector(OK1,OK2,alpha,K1);
   staFunc->NormalizeVector(K1);
   alpha=staFunc->ScalProd(V12,K1);
   if (alpha<0) {
    for (i=0;i<m+2;i++)
      K1[i]=-K1[i];
   }
   staFunc->CopyVector(K1,V12);

   staFunc->TransposeMatrix(BP2,BT2);
   staFunc->ClearVector(K1);
   staFunc->ClearVector(K2);
   K1[m]=1.0;
   K2[m+1]=1.0;
   staFunc->Decomp(BT2);
   staFunc->Solve(BT2,K1);
   staFunc->Solve(BT2,K2);
   staFunc->CopyVector(K1,OK1);
   staFunc->CopyVector(K1,K3); 
   staFunc->CopyVector(K2,OK2);
   max1=(fabs(OK1[m])>fabs(OK2[m]))?fabs(OK1[m]):fabs(OK2[m]);
   max2=(fabs(OK1[m+1])>fabs(OK2[m+1]))?fabs(OK1[m+1]):fabs(OK2[m+1]);
   if (max1>max2) {
    alpha=OK1[m];  
    beta=OK2[m];   
   }
   else {
    alpha=OK1[m+1];
    beta=OK2[m+1]; 
   }
   staFunc->MultiplyVectorScalar(OK1,-beta,OK1);
   staFunc->AddScaledVector(OK1,OK2,alpha,K1);  
   staFunc->NormalizeVector(K1);
   alpha=staFunc->ScalProd(W12,K1);
   if (alpha<0) {
    for (i=0;i<m+2;i++)
      K1[i]=-K1[i];
   }
   staFunc->CopyVector(K1,W12);
/* adapt V22 */
   staFunc->ClearMatrix(BP1);
   staFunc->CopyMatrixBlock(BP2,0,0,BP1,0,0,m,m);
   for (i=0;i<m;i++) {
    BP1[i][m]=W12[i]; 
    BP1[m][i]=V12[i]; 
   }
   staFunc->CopyVectorElements(W22,0,RL,0,m+1);
   Ret=staFunc->Decomp(BP1);
   if (Ret != 0) {
    return 1;
   }
   staFunc->Solve(BP1,RL);
   staFunc->NormalizeVector(RL);
   staFunc->CopyVectorElements(V22,0,OLD,0,m+1);
   alpha=staFunc->ScalProd(OLD,RL);
   if (alpha < 0) {
    for (i=0;i<m+1;i++)
      RL[i]=-RL[i];
   }
   staFunc->ClearVector(V22);
   staFunc->CopyVectorElements(RL,0,V22,0,m+1);
/* adapt W22 */
   staFunc->ClearMatrix(BP1);
   staFunc->CopyMatrixBlock(BP2,0,0,BP1,0,0,m,m);
   for (i=0;i<m;i++) {
    BP1[i][m]=W12[i]; 
    BP1[m][i]=V12[i]; 
   }
   staFunc->TransposeMatrix(BP1,BT1);
   staFunc->ClearVector(RL);
   staFunc->CopyVectorElements(OLD,0,RL,0,m+1);
   Ret=staFunc->Decomp(BT1);
   if (Ret != 0) {
    return 1;
   }
   staFunc->Solve(BT1,RL);
   staFunc->NormalizeVector(RL);
   staFunc->CopyVectorElements(W22,0,OLD,0,m+1);
   alpha=staFunc->ScalProd(OLD,RL);
   if (alpha<0) {
    for (i=0;i<m+1;i++)
     RL[i]=-RL[i];
   }
   staFunc->ClearVector(W22);
   staFunc->CopyVectorElements(RL,0,W22,0,m+1);
/* border biproduct matrix with V12,V22,W12 and W22 */
   for (i=0;i<m;i++) {
    BP2[m][i]=V12[i]; 
    BP2[i][m]=W12[i]; 
   }
   for (i=0;i<m+2;i++) {
    BP2[m+1][i]=V22[i]; 
    BP2[i][m+1]=W22[i]; 
   }
  
  /* Compute a tangent vector to the BT curve */
    zh_DefZeroHopf(x0,v0);
    zh_JacDefZeroHopf(x0,D);
    for (i=0; i<nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+2+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 zero-Hopf 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;
}


typedef EntryPtr(Int2,TestFuncPtr,(Vector xx, Vector vv, FloatHiPtr TV));
typedef EntryPtr(Int2,ProcFuncPtr,(Int2 Ind, Vector xx, Vector vv, CharPtr *msg));


/************************/
/* 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<3+nappend; i++) P[active[i]]=vp[nphase+i];
}

Local(void) CopyToE(Vector Re, Vector Im) {
  MemCopy(E,Re,sizeof(FloatHi)*nphase);
  MemCopy(E+nphase,Im,sizeof(FloatHi)*nphase);
}


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

Entry(Int2) zh_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_EIGV]=E;
      CopyToE(Re,Im);
      indirect[CL_TEST]=test;
  }
  return 0;
}
/***************************/
/* User-defined functions  */
Entry(Int2) zh_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) zh_DefRHS (Vector x, Vector y) {
  CopyToX(x);
  CopyToP(x);
  stagenData->Rhs(X,P,&zero,y);	
  return 0;
}

/*************************************/
/* Defining function for equilibrium */
Local(Int2) DefEquilib (Vector x, Vector y) {
  zh_DefRHS(x,y);
  if (nappend) zh_DefUser(x,y+nphase);
  return 0;
}

/********************/
/* Jacobian of RHS  */
Local(Int2) JacRHS (Vector x, Matrix Jac) {
  staFunc->CopyVectorElements(x,0,xp,0,nphase+3+nappend);     
  if (stagenData->Der1) {
    SymJac(xp,Jac);
  } else {
    if (staFunc->Der(zh_DefRHS, nphase, xp, dx, Jac)) return 1;
  }
  return 0;
}
/***************************************************/
/* Jacobian of defining conditions for equilibria  */
Local(Int2) JacDefEquilib (Vector x, Matrix Jac) {
  staFunc->ClearMatrix(Jac);
  JacRHS(x, Jac);
  if (nappend) staFunc->Der(zh_DefUser, nappend, x, dx, Jac+nphase);
  return 0;
}
/***********************************/
/* Zero-Hopf defining conditions */
Entry(Int2) zh_DefBord (Vector x, Vector y) {
Int2 i,j,m,p,q,r,s,Ret;

/* Creation of biproduct matrix BP2 */
  m=nphase*(nphase-1)/2;
  if (JacRHS(x,C1)) return 1;	    /* A,C1,BP,V1,W1,Q - global */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->ClearMatrix(BP2);
  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++;
      if (r==q) BP2[i][j]=-A[p][s];
      if (r!=p && s==q) BP2[i][j]=A[p][r];
      if (r==p && s==q) BP2[i][j]=A[p][p]+A[q][q];
      if (r==p && s!=q) BP2[i][j]=A[q][s];
      if (s==p) BP2[i][j]=-A[q][r];
     }
    }
   }
  }
  for (i=0;i<m;i++) {
   BP2[m][i]=V12[i];
   BP2[i][m]=W12[i];
  }
  for (i=0;i<m+2;i++) {
   BP2[m+1][i]=V22[i];
   BP2[i][m+1]=W22[i];
  }  
  staFunc->CopyMatrix(BP2,BT2);
  Ret=staFunc->Decomp(BT2);
  if (Ret!=0) {
    return 1;
  }
  staFunc->ClearVector(VB12);
  VB12[m]=1.0;
  staFunc->Solve(BT2,VB12);
  staFunc->ClearVector(VB22);
  VB22[m+1]=1.0;
  staFunc->Solve(BT2,VB22);
  y[0]=VB12[m]*VB22[m+1]-VB12[m+1]*VB22[m];
  
/* Fill in the bordered Jacobian BP */
  staFunc->ClearMatrix(BP);
  staFunc->CopyMatrixBlock(C1,0,0,BP,0,0,nphase,nphase);

/* 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 BT defining functions */
  staFunc->CopyMatrix(BP,BT);
  Ret=staFunc->Decomp(BT);
  if (Ret != 0) {
    return 1;
  }
  staFunc->ClearVector(VB11);
  VB11[nphase]=1.0;
  staFunc->Solve(BT,VB11);
  staFunc->CopyVectorElements(VB11,0,Q,0,nphase);   /* make global null-vector  AQ=0 */
  y[1]=VB11[nphase];
  return 0;
}

/***************************************/
/* Jacobian of ZH defining condition */
Local(Int2) JacDefBord (Vector x, Matrix Jac) {
int i,j,l,m,p,r,s,q;
VAR S,s11,s12,s21,s22;

  m=nphase*(nphase-1)/2;
  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2) {
/* creation of DefBord Jacobian matrix Jac using symbolic second derivatives */
    CopyToX(x);
    CopyToP(x);
    H=stagenData->Der2(X,P,&zero);
/* take global vectors VB12 and VB22 computed at the computation of DefBord */
/* compute PP12 and PP22 using the transposed bordered bi-product */
    staFunc->TransposeMatrix(BP2,BT2);
    staFunc->ClearVector(PP12);
    PP12[m]=1.0;
    staFunc->Decomp(BT2);
    staFunc->Solve(BT2,PP12);
    staFunc->ClearVector(PP22);
    PP22[m+1]=1.0;
    staFunc->Solve(BT2,PP22);
    
/* compute gradient of DefBord */
    for (l=0; l<nphase+3+nappend; l++) {
      staFunc->ClearMatrix(BT2);
      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++;
              if (r == q) BT2[i][j]=-HessianElem(p,s,l);
              if (r != p && s == q) BT2[i][j]=HessianElem(p,r,l);
              if (r == p && s == q) BT2[i][j]=HessianElem(p,p,l)+HessianElem(q,q,l);
              if (r == p && s != q) BT2[i][j]=HessianElem(q,s,l);
              if (s == p) BT2[i][j]=-HessianElem(q,r,l);
            }
          }
        } 
      }
      for (i=0,s11=0.0,s12=0.0,s21=0.0,s22=0.0; i<m; i++) {
        for (j=0; j<m; j++) {
          s11-=PP12[i]*VB12[j]*BT2[i][j];
          s12-=PP12[i]*VB22[j]*BT2[i][j];
          s21-=PP22[i]*VB12[j]*BT2[i][j];
          s22-=PP22[i]*VB22[j]*BT2[i][j];
        }
      }
      Jac[0][l]=VB12[m]*s22+VB22[m+1]*s11-VB22[m]*s21-VB12[m+1]*s12;
    }

/* take global vector Q generated at the computation of DefBord */
/* compute PP11 using the transposed bordered bi-product */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->Decomp(BT);

    staFunc->ClearVector(PP11);
    PP11[nphase]=1.0;
    staFunc->Solve(BT,PP11);
    staFunc->ClearVector(PP21);
    staFunc->CopyVectorElements(PP11,0,PP21,0,nphase);
    staFunc->Solve(BT,PP21);
    staFunc->CopyMatrix(BP,BT);
    staFunc->Decomp(BT);
      
/* compute gradient of DefBord */
    for (l=0; l<nphase+3+nappend; l++) {
      for (i=0,S=0.0; i<nphase; i++) {
        for (j=0; j<nphase; j++) {
	  S-=PP11[i]*Q[j]*HessianElem(i,j,l);
	}
      }
      Jac[1][l]=S;
    }
  } else {
    if (staFunc->Der(zh_DefBord, 2, x, dx, Jac)) return 1;
  }
  return(0);
}

/*************************************/
/* Defining functions for BT curve */

Entry(Int2) zh_DefZeroHopf (Vector x, Vector y) {
  if (DefEquilib(x,y)) return 1;
  if (zh_DefBord(x,y+nphase+nappend)) return 2;
  return 0;
}

/*************************************************/
/* Jacobian of defining functions for BT curve */  
Entry(Int2) zh_JacDefZeroHopf (Vector x, Matrix Jac) {
  staFunc->ClearMatrix(Jac);
  if (JacDefEquilib(x,C)) return 1;
  staFunc->CopyMatrixBlock(C,0,0,Jac,0,0,nphase+nappend,nphase+3+nappend);
  if (JacDefBord(x,C2)) return 2;
  staFunc->CopyMatrixBlock(C2,0,0,Jac,nphase+nappend,0,2,nphase+3+nappend);
  return 0;
}

/*********************/
/* Default processor */

Entry(Int2) zh_ProcDefault(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {

  zh_DefBord (x, xp);					/* update of global Q and VB1 */
  /* Assert: v1==NULL, msg==NULL, Ind==0 */
  MemCopy(x1,x,sizeof(FloatHi)*(nphase+3+nappend));     /* for the last point */
  if (!staData->eigcomp) return 0;
  if (JacRHS(x,C1)) return (1); 
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->EigenValues(A,Re,Im);
  return 0;
}

/*********************************/
/* Test and processing functions */
Local(Char) buf[80];

/* Test functions */
/* Test function to rank deficiency 2 of BP2*/
Entry(Int2) zh_TestR2(Vector x, Vector v, FloatHiPtr res) {
int m;

  m=nphase*(nphase-1)/2;

  *res=VB22[m+1];
  test[0]=*res;
  
  return 0;
}

Entry(Int2) zh_ProcDH(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Double Hopf (?)"); 
  else {  
     {int i,j,imin,jmin,imin0,jmin0,imin1,jmin1,branch;
 VAR Ret,r,s,d,phi,beta,gamma,lambda,omega,kappa;
    Ret=0;
    JacRHS(x,C1);
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    /* find two pairs of eigenvalues with zero sum */
    s=fabs(Re[0]+Re[1])+fabs(Im[0]+Im[1]);
    imin0=0;
    jmin0=1;
    for (i=0; i<nphase; i++) {
      for (j=i+1; j<nphase; j++) {
        d=fabs(Re[i]+Re[j])+fabs(Im[i]+Im[j]);
        if (d<s) {
          imin0=i;
          jmin0=j;
          s=d;
        }
      }
    }
    imin1=imin0+1;
    if (imin1==nphase) imin1=0; 
    jmin1=jmin0+1;
    if (jmin1==nphase) jmin1=0; 
    s=fabs(Re[imin1]+Re[jmin1])+fabs(Im[imin1]+Im[jmin1]);
    for (i=0; i<nphase; i++) {
      if (i!=imin0) {
        for (j=i+1; j<nphase; j++) {
	  if (j!=jmin0) {
            d=fabs(Re[i]+Re[j])+fabs(Im[i]+Im[j]);
            if (d<s) {
              imin1=i;
              jmin1=j;
              s=d;
            }
	  }
	}
      }
    }
/* printf("imin0=%i, jmin0=%i, imin1=%i, jmin1=%i\n",imin0,jmin0,imin1,jmin1); */ 

  for (branch=0; branch<2; branch++) {
    if (branch == 0) {
    /* primary branch */
      imin=imin0;
      jmin=jmin0;
    } else {
    /* secondary branch */
      imin=imin1;
      jmin=jmin1;
    }
/* printf("imin=%i, jmin=%i\n",imin,jmin); */
    if (fabs(Im[imin])+fabs(Im[jmin])==0) {
/* Neutral saddle */
      lambda=fabs(Re[imin]);
      /* compute real eigenvectors associated with lambda and -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,V) != nphase-1) {
        Ret=1;
        goto polufinal;
      } 
      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,W) != nphase-1) {
        printf("no real eigenvectors(+)\n");
        Ret=1;
        goto polufinal;
      } 
      staFunc->NormalizeVector(V);   
      staFunc->NormalizeVector(W);
      staFunc->AddScaledVector(V,W,-1.0,L);             
      staFunc->NormalizeVector(L);
      staFunc->AddScaledVector(V,W,1.0,V);             
      staFunc->NormalizeVector(V);
      kappa=-lambda*lambda;
    } else {
/* Hopf */
      lambda=Re[imin];
      omega=fabs(Im[imin]);
      /* compute real and imaginary part of the eigenvector */
      staFunc->ClearMatrix(BB);
      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) {
        Ret=1;
        goto polufinal;
      }
      staFunc->CopyVectorElements(VV,0,V,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,W,0,nphase);
      /* normalize real and imaginary parts */
        d=staFunc->ScalProd(V,V);
        s=staFunc->ScalProd(W,W);
        r=staFunc->ScalProd(V,W);
        if (r!=0) {
          beta=sqrt(fabs(4.0*r*r+(s-d)*(s-d)));
          if (2.0*fabs(2.0*r) > beta) {
            phi=atan(fabs((s-d)/(2.0*r)));
          } else {
            phi=PID2-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]=V[i]*cos(phi)-W[i]*sin(phi);
            VV[i+nphase]=V[i]*sin(phi)+W[i]*cos(phi);
          }
        }
        staFunc->CopyVectorElements(VV,0,V,0,nphase);
        d=staFunc->ScalProd(V,V);
        gamma=1/sqrt(d);
        staFunc->MultiplyVectorScalar(VV,gamma,VV);

        staFunc->CopyVectorElements(VV,0,V,0,nphase);
        staFunc->CopyVectorElements(VV,nphase,L,0,nphase);
        staFunc->NormalizeVector(L);
        kappa=omega*omega;
    }
polufinal:
    if (Ret) {
        sprintf(buf,"Double Hopf: eigenvectors (?)"); 
        *msg=buf;
        return 0;
    }
/* provide bifurcation data for dh->h starter 
  Primary branch: 
  BifVector[0] to [nphase-1]:          1-Hopf continuation vector 
  BifVector[nphase] to [2*nphase-1]:   1-global vector L 
  BifVector[2*nphase]:                 -lambda[1]*lambda[1]  
  Secondary branch: 
  BifVector[2*nphase+1] to [3*nphase]: 2-Hopf continuation vector
  BifVector[3*nphase+1] to [4*nphase]: 2-global vector L 
  BifVector[4*nphase+1]:               -lambda[2]*lambda[2]  
*/
    MemCopy(BifVector+(2*nphase+1)*branch,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+(2*nphase+1)*branch,L,sizeof(FloatHi)*nphase);
    BifVector[2*nphase+(2*nphase+1)*branch]=kappa;    
  } /* for branch */
  stagenData->BifDataOutLen=sizeof(FloatHi)*(4*nphase+2);
}
  if (BifVector[2*nphase]<=0.0 && BifVector[4*nphase+1]<=0.0)  sprintf(buf,
    "Double neutral saddle: lambda[1]=%g, lambda[2]=%g",sqrt(fabs(BifVector[2*nphase])),
                                                        sqrt(fabs(BifVector[4*nphase+1])));
  if (BifVector[2*nphase]>0.0 && BifVector[4*nphase+1]>0.0)  sprintf(buf,
    "Double Hopf: omega[1]=%g, omega[2]=%g",sqrt(BifVector[2*nphase]),
                                            sqrt(BifVector[4*nphase+1]));
  if (BifVector[2*nphase]<=0.0 && BifVector[4*nphase+1]>0.0)  sprintf(buf,
    "Neutral saddle + Hopf: lambda[1]=%g, omega[2]=%g",sqrt(fabs(BifVector[2*nphase])),
                                                       sqrt(BifVector[4*nphase+1]));
  if (BifVector[2*nphase]>0.0 && BifVector[4*nphase+1]<=0.0)  sprintf(buf,
    "Hopf + neutral saddle: omega[1]=%g ,lambda_2=%g",sqrt(BifVector[2*nphase]),
                                                      sqrt(fabs(BifVector[4*nphase+1])));
  } 
  *msg=buf;
  return 0; 
}

/* Test function for Bogdanov-Takens */
Entry(Int2) zh_TestBT(Vector x, Vector v, FloatHiPtr res) {

    staFunc->CopyMatrix(BP,BT);
    staFunc->Decomp(BT);
    staFunc->ClearVector(VB11);
    staFunc->CopyVectorElements(Q,0,VB11,0,nphase);
    staFunc->Solve(BT,VB11);
        
    *res=VB11[nphase];
    test[1]=*res;
    
    return 0;
}

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

  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Hopf-BT (?)");
  else {
    Vector v1,W0,W1;
    Matrix D;
    Int2 i,j,k;
    VAR a20,b20,b11,f0,f1,f2,f3,f4,f5,d;
     v1=staFunc->CreateVector(2*nphase);
     W0=staFunc->CreateVector(nphase);
     W1=staFunc->CreateVector(nphase);
     D=staFunc->CreateMatrix(2*nphase,2*nphase);

/*   compute (generalized)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,"Bogdanov-Takens: eigenvectors (?)");
       goto final;
     }
     staFunc->CopyVectorElements(v1,0,V,0,nphase);       /* null-vector */
     d=sqrt(staFunc->ScalProd(V,V));
     staFunc->NormalizeVector(V);
     staFunc->CopyVectorElements(v1,nphase,W,0,nphase);  /* generalized null-vector */
     staFunc->MultiplyVectorScalar(W,1.0/d,W);
     staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W);

/*  provide bifurcation data for zh->* 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,"Bogdanov-Takens: 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);

/* compute the coefficients of the normal form : x'=y; y'=ax^2+bxy */
/*  a20=<W0,B(V,V)>, b20=<W1,B(V,V)>,  b11=<W1,B(V,W)> */
      if (stagenData->Der2)  {
         CopyToX(x);
         CopyToP(x);
         H=stagenData->Der2(X,P,&zero);	/* The 1st elem must be time! */
         for (a20=0,b20=0,b11=0,i=0; i<nphase; i++) {
            for (j=0; j<nphase; j++) {
               for (k=0; k<nphase; k++) { 
	          d=HessianElem(i,j,k);
                  a20+=W0[i]*d*V[j]*V[k];
                  b20+=W1[i]*d*V[j]*V[k];
                  b11+=W1[i]*d*V[j]*W[k];
 	       }
            }
         }

      } else {
 	 staFunc->ClearVector(v0);
/*  a20=<W0,B(V,V)>, b20=<W1,B(V,V)> */
 	 staFunc->CopyVectorElements(V,0,v0,0,nphase);
         zh_DefRHS(x,F);
         f0=staFunc->ScalProd(W0,F);
         f3=staFunc->ScalProd(W1,F);
         staFunc->AddScaledVector(x,v0,ddx,x1);
         zh_DefRHS(x1,F);
         f1=staFunc->ScalProd(W0,F);
         f4=staFunc->ScalProd(W1,F);
         staFunc->AddScaledVector(x,v0,-ddx,x1);
         zh_DefRHS(x1,F);
         f2=staFunc->ScalProd(W0,F);
         f5=staFunc->ScalProd(W1,F);
	 a20=(f1-2*f0+f2)/ddx/ddx;
	 b20=(f4-2*f3+f5)/ddx/ddx;

/*  b11=<W1,B(V,W)>  */
         staFunc->AddScaledVector(V,W,1.0,F);
         staFunc->CopyVectorElements(F,0,v0,0,nphase);
         staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
         zh_DefRHS(x1,F);
	 b11=staFunc->ScalProd(W1,F);
         staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
         zh_DefRHS(x1,F);
	 b11+=staFunc->ScalProd(W1,F);
         staFunc->AddScaledVector(V,W,-1.0,F);
         staFunc->CopyVectorElements(F,0,v0,0,nphase);
         staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
         zh_DefRHS(x1,F);
	 b11-=staFunc->ScalProd(W1,F);
         staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
         zh_DefRHS(x1,F);
	 b11-=staFunc->ScalProd(W1,F);
	 b11/=(ddx*ddx); 
      }
     if (b20>=0)
       sprintf(buf,"Bogdanov-Takens: a=%g, b=%g",b20/2,a20+b11);
     else
       sprintf(buf,"Bogdanov-Takens: a=%g, b=%g",-b20/2,-(a20+b11));
final:
     staFunc->DestroyMatrix(D);
     staFunc->DestroyVector(v1);
     W0=staFunc->DestroyVector(W0);
     W1=staFunc->DestroyVector(W1);
  } 
  *msg=buf;
  return 0;
}


Entry(Int2) zh_ProcTZ(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
  if (Ind==0) {
    sprintf(buf,"Triple zero eigenvalue");
  } else {
    sprintf(buf,"Triple zero eigenvalue (?)");
  }
  *msg=buf;
  return 0;
}

/**************************************************/
/*  Adaptation of the bordering vectors V1 and W1 */
Entry(Int2) zh_Adapter(Vector x, Vector v) {
int i,m;
Int2 Ret;
VAR alpha,beta,max1,max2;
 
  m=nphase*(nphase-1)/2;   
/* global Q is assumed to be correct */
  staFunc->CopyVector(Q,V11);
  staFunc->NormalizeVector(V11);

/* global BP is assumed to be correct */
  staFunc->TransposeMatrix(BP,BT);
  staFunc->Decomp(BT);

  staFunc->ClearVector(VB21);
  VB21[nphase]=1.0;
  staFunc->Solve(BT,VB21);
  staFunc->CopyVectorElements(VB21,0,W11,0,nphase);
  staFunc->NormalizeVector(W11);

  staFunc->ClearVector(K1);
  staFunc->ClearVector(K2);
  staFunc->CopyMatrix(BP2,BT2);
  K1[m]=1.0;
  K2[m+1]=1.0;
  staFunc->Decomp(BT2);
  staFunc->Solve(BT2,K1);
  staFunc->Solve(BT2,K2);
  staFunc->CopyVector(K1,OK1);
  staFunc->CopyVector(K1,K3);
  staFunc->CopyVector(K2,OK2);
  max1=(fabs(OK1[m])>fabs(OK2[m]))?fabs(OK1[m]):fabs(OK2[m]);
  max2=(fabs(OK1[m+1])>fabs(OK2[m+1]))?fabs(OK1[m+1]):fabs(OK2[m+1]);
  if (max1>max2) {
   alpha=OK1[m];
   beta=OK2[m];
  }
  else {
   alpha=OK1[m+1];
   beta=OK2[m+1];
  }
  staFunc->MultiplyVectorScalar(OK1,-beta,OK1);
  staFunc->AddScaledVector(OK1,OK2,alpha,K1);
  staFunc->NormalizeVector(K1);
  alpha=staFunc->ScalProd(V12,K1);
  if (alpha<0) {
   for (i=0;i<m+2;i++)
     K1[i]=-K1[i];
  }
  staFunc->CopyVector(K1,V12);
  
  staFunc->TransposeMatrix(BP2,BT2);
  staFunc->ClearVector(K1);
  staFunc->ClearVector(K2);
  K1[m]=1.0;
  K2[m+1]=1.0;
  staFunc->Decomp(BT2);
  staFunc->Solve(BT2,K1);
  staFunc->Solve(BT2,K2);
  staFunc->CopyVector(K1,OK1);
  staFunc->CopyVector(K1,K3);
  staFunc->CopyVector(K2,OK2);
  max1=(fabs(OK1[m])>fabs(OK2[m]))?fabs(OK1[m]):fabs(OK2[m]);
  max2=(fabs(OK1[m+1])>fabs(OK2[m+1]))?fabs(OK1[m+1]):fabs(OK2[m+1]);
  if (max1>max2) {
   alpha=OK1[m];
   beta=OK2[m];
  }
  else {
   alpha=OK1[m+1];
   beta=OK2[m+1];
  }
  staFunc->MultiplyVectorScalar(OK1,-beta,OK1);
  staFunc->AddScaledVector(OK1,OK2,alpha,K1);
  staFunc->NormalizeVector(K1);
  alpha=staFunc->ScalProd(W12,K1);
  if (alpha<0) {
   for (i=0;i<m+2;i++)
     K1[i]=-K1[i];
  }
  staFunc->CopyVector(K1,W12);
/* adapt V22 */   
  staFunc->ClearMatrix(BP1);
  staFunc->CopyMatrixBlock(BP2,0,0,BP1,0,0,m,m);
  for (i=0;i<m;i++) {
   BP1[i][m]=W12[i];
   BP1[m][i]=V12[i];
  }
  staFunc->CopyVectorElements(W22,0,RL,0,m+1);
  Ret=staFunc->Decomp(BP1);
  if (Ret != 0) {
   return 1;
  }
  staFunc->Solve(BP1,RL);
  staFunc->NormalizeVector(RL);
  staFunc->CopyVectorElements(V22,0,OLD,0,m+1);
  alpha=staFunc->ScalProd(OLD,RL);
  if (alpha < 0) {
   for (i=0;i<m+1;i++)
     RL[i]=-RL[i];
  }
  staFunc->ClearVector(V22);
  staFunc->CopyVectorElements(RL,0,V22,0,m+1);
/* adapt W22 */   
  staFunc->ClearMatrix(BP1);
  staFunc->CopyMatrixBlock(BP2,0,0,BP1,0,0,m,m);
  for (i=0;i<m;i++) {
   BP1[i][m]=W12[i];
   BP1[m][i]=V12[i];
  }
  staFunc->TransposeMatrix(BP1,BT1);
  staFunc->ClearVector(RL);
  staFunc->CopyVectorElements(V22,0,RL,0,m+1);
  Ret=staFunc->Decomp(BT1);
  if (Ret != 0) {
   return 1;
  }
  staFunc->Solve(BT1,RL);
  staFunc->NormalizeVector(RL);
  staFunc->CopyVectorElements(W22,0,OLD,0,m+1);
  alpha=staFunc->ScalProd(OLD,RL);
  if (alpha<0) {
   for (i=0;i<m+1;i++)
    RL[i]=-RL[i];
  }    
  staFunc->ClearVector(W22);
  staFunc->CopyVectorElements(RL,0,W22,0,m+1);

  for (i=0;i<m;i++) { 
   BP2[m][i]=V12[i];
   BP2[i][m]=W12[i];
  }
  for (i=0;i<m+2;i++) {  
   BP2[m+1][i]=V22[i];
   BP2[i][m+1]=W22[i];
  }
  staFunc->CopyMatrix(BP2,BT2);
  staFunc->Decomp(BT2);
  staFunc->ClearVector(VB22);
  VB22[m+1]=1.0;
  staFunc->Solve(BT2,VB22);
  test[0]=VB22[m+1];  

  return 0;
}