/*
   __h1.c
   Hopf(neutral saddle) point-->Hopf(neutral saddle) Curve 
   using double-bordering biproduct method by W. Govaerts and B. Sijnave
   September 29, 1997
*/

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

#include "stagen.h"

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

#include "__h1.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 _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(__h_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,AT,BB,BP,BT,D,DE,BP1,BT1;
Local(Vector) vn,xp,xp1,V,W,F,Vp,L,VV,VB,V2,W2,Q,Qr,Qi,PP1,PP2,Pr,Pi,aa,bb,cc,dd,rr;
Local(Vector) VB1,VB2,V1,W1,OLD,RL,K1,K2,K3,OK1,OK2,PIVOT;
Local(Vector) Re=NULL,Im=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=(__h_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+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;
  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(m+2,m+2);
    BT=staFunc->CreateMatrix(m+2,m+2);
    BP1=staFunc->CreateMatrix(m+1,m+1);
    BT1=staFunc->CreateMatrix(m+1,m+1);
    
    V1=staFunc->CreateVector(m+2);
    W1=staFunc->CreateVector(m+2);    
    V2=staFunc->CreateVector(m+2);
    W2=staFunc->CreateVector(m+2);
    VB=staFunc->CreateVector(m+2);
    VB1=staFunc->CreateVector(m+2);
    VB2=staFunc->CreateVector(m+2);    
    VV=staFunc->CreateVector(2*nphase);
    Q=staFunc->CreateVector(m);
    PP1=staFunc->CreateVector(m+2);
    PP2=staFunc->CreateVector(m+2);
    OLD=staFunc->CreateVector(m+1);
    RL=staFunc->CreateVector(m+1);
    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);
    PIVOT=staFunc->CreateVector(m);

    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);
    vn=staFunc->CreateVector(nphase);
    xp=staFunc->CreateVector(nphase+2+nappend);
    xp1=staFunc->CreateVector(nphase+2+nappend);
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    F=staFunc->CreateVector(nphase);
    L=staFunc->CreateVector(nphase);
    Vp=staFunc->CreateVector(nphase+2+nappend);
    BifVector=MemNew(sizeof(FloatHi)*(4*nphase+2));
    if (nphase==1) staData->ActiveSing[0]=FALSE;
    if ((nphase>1) && staData->ActiveSing[0]) {  /* Bautin detection ON*/
      AT=staFunc->CreateMatrix(nphase,nphase);
      Qr=staFunc->CreateVector(nphase);
      Qi=staFunc->CreateVector(nphase);
      Pr=staFunc->CreateVector(nphase);
      Pi=staFunc->CreateVector(nphase);
      aa=staFunc->CreateVector(nphase);
      bb=staFunc->CreateVector(nphase);
      cc=staFunc->CreateVector(nphase);
      dd=staFunc->CreateVector(nphase);
      rr=staFunc->CreateVector(nphase);
    }
    if (nphase<=3) staData->ActiveSing[3]=FALSE;
    if ((nphase>3) && staData->ActiveSing[3]) {  /*Double Hopf detection ON*/
      D=staFunc->CreateMatrix(nphase+1+nappend,nphase+2+nappend);
      DE=staFunc->CreateMatrix(nphase+2+nappend,nphase+2+nappend);
    }
  }
  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);
  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 max1,max2,d,a11,a12,a21,a22;

  stagenData->BifDataOutLen=0;
  m=nphase*(nphase-1)/2;
  if (OK&&P&&(stagenData->Reason==RF_CLEAN)) {
/* saving global vectors V1,V2,W1 and W2 used in defining functions */
    stagenData->ContDataStaLen=4*(m+2)*sizeof(FloatHi);
    stagenData->ContDataStaPtr=MemNew(stagenData->ContDataStaLen);
    MemCopy(stagenData->ContDataStaPtr,V1,(m+2)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+(m+2),W1,(m+2)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+2*(m+2),V2,(m+2)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+3*(m+2),W2,(m+2)*sizeof(FloatHi));    

/* reevaluate Jacobian */
    JacRHS (x1,C1);
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    
/* find pair in (g11,g12), (g21,g22) with max absolute value */
    max1=(fabs(VB1[m])>fabs(VB2[m]))?VB1[m]:VB2[m];
    max2=(fabs(VB1[m+1])>fabs(VB2[m+1]))?VB1[m+1]:VB2[m+1];
    if (fabs(max1)>fabs(max2))
     imax=m;
    else imax=m+1;    
    for(i=0;i<m;i++) {
      Q[i]=VB2[imax]*VB1[i]-VB1[imax]*VB2[i];
    }
       
/* find a component of Q with max(fabs) */
    imax=1;
    jmax=0;
    d=fabs(Q[0]);
    for (i=0; i<nphase; i++) {
      for (j=0; j<i; j++) {
        if (fabs(Q[i*(i-1)/2+j])>d) {
          imax=i;
	  jmax=j;
	  d=fabs(Q[i*(i-1)/2+j]);
        }
      }
    }
    d=Q[imax*(imax-1)/2+jmax];

/* compute eigenvectors */
    staFunc->ClearVector(V);
    staFunc->ClearVector(W);
    for (i=0; i<nphase; i++) {
      if (i<imax) V[i]=Q[imax*(imax-1)/2+i]/d;
      if (i>imax) V[i]=-Q[i*(i-1)/2+imax]/d;
      if (i<jmax) W[i]=-Q[jmax*(jmax-1)/2+i]/d;
      if (i>jmax) W[i]=Q[i*(i-1)/2+jmax]/d;
    }

/* compute -lambda*lambda */
    staFunc->MultiplyMatrixVector(A,V,F);  /* F=A*V  */
    a11=staFunc->ScalProd(V,F);
    a12=staFunc->ScalProd(W,F);
    staFunc->MultiplyMatrixVector(A,W,F);  /* F=A*W  */
    a21=staFunc->ScalProd(W,F);
    a22=staFunc->ScalProd(V,F);

    d=staFunc->ScalProd(V,V)*staFunc->ScalProd(W,W)
     -staFunc->ScalProd(V,W)*staFunc->ScalProd(V,W);
  
/* orthogonalize and normalize eigenvectors */
    staFunc->NormalizeVector(V);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),L);
    staFunc->NormalizeVector(L);
  
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1);
    MemCopy(BifVector,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase,L,sizeof(FloatHi)*nphase);
    BifVector[2*nphase]=(a11*a21-a22*a12)/d;
    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);
    BP1=staFunc->DestroyMatrix(BP1);
    BT1=staFunc->DestroyMatrix(BT1);
    V1=staFunc->DestroyVector(V1);
    W1=staFunc->DestroyVector(W1);    
    V2=staFunc->DestroyVector(V2);
    W2=staFunc->DestroyVector(W2);
    Q=staFunc->DestroyVector(Q);
    PP1=staFunc->DestroyVector(PP1);
    PP2=staFunc->DestroyVector(PP2);
    AA=staFunc->DestroyMatrix(AA);
    BB=staFunc->DestroyMatrix(BB);
    VV=staFunc->DestroyVector(VV);
    VB=staFunc->DestroyVector(VB);
    VB1=staFunc->DestroyVector(VB1);
    VB2=staFunc->DestroyVector(VB2);    
    OLD=staFunc->DestroyVector(OLD);
    RL=staFunc->DestroyVector(RL);    
    K1=staFunc->DestroyVector(K1);
    K2=staFunc->DestroyVector(K2);
    K3=staFunc->DestroyVector(K3);
    OK1=staFunc->DestroyVector(OK1);
    OK2=staFunc->DestroyVector(OK2); 
    PIVOT=staFunc->DestroyVector(PIVOT);
    C=staFunc->DestroyMatrix(C);
    C1=staFunc->DestroyMatrix(C1);
    C2=staFunc->DestroyMatrix(C2);
    vn=staFunc->DestroyVector(vn);
    xp=staFunc->DestroyVector(xp);
    xp1=staFunc->DestroyVector(xp1);
    V=staFunc->DestroyVector(V);
    W=staFunc->DestroyVector(W);
    F=staFunc->DestroyVector(F);
    Vp=staFunc->DestroyVector(Vp);
    L=staFunc->DestroyVector(L);
    BifVector=MemFree(BifVector);
    if ((nphase>1) && staData->ActiveSing[0]) {  /* Bautin detection ON*/
      AT=staFunc->DestroyMatrix(AT);
      Qr=staFunc->DestroyVector(Qr);
      Qi=staFunc->DestroyVector(Qi);
      Pr=staFunc->DestroyVector(Pr);
      Pi=staFunc->DestroyVector(Pi);
      aa=staFunc->DestroyVector(aa);
      bb=staFunc->DestroyVector(bb);
      cc=staFunc->DestroyVector(cc);
      dd=staFunc->DestroyVector(dd);
      rr=staFunc->DestroyVector(rr);
    }
    if ((nphase>3) && staData->ActiveSing[3]) {  /* Double Hopf detection ON*/
      D=staFunc->DestroyMatrix(D);
      DE=staFunc->DestroyMatrix(DE);
    }
  }
} 

Local(Int2) JacDefEquilib (Vector x, Matrix J);
Local(Int2) JacDefEig(Vector x, Matrix J);
Entry(Int2) h1_JacDefHopf (Vector x, Matrix J);
Entry(Int2) h1_DefHopf (Vector x, Vector J);


/************/
/* Starters */
/* h -> h1 Starter */
Entry(Int2) h_h1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,Ret,m,p,q,r,s,imax;
  Vector v1,v2;
  Matrix D,DE,G;
  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(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(m+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(m+2),sizeof(FloatHi)*(m+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(m+2),sizeof(FloatHi)*(m+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(m+2),sizeof(FloatHi)*(m+2));    
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];    
    
    /* compute Q : needed for extend of Zero Hopf */
      if (JacRHS(x0,C1)) return 1; 	  /* A,C1,BP,V1,V2,W1,W2 - global */
      staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
        
      /* Fill in the bordered bi-product matrix BP */
      staFunc->ClearMatrix(BP);
      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) BP[i][j]=-A[p][s];       
             if (r != p && s == q) BP[i][j]=A[p][r];   
             if (r == p && s == q) BP[i][j]=A[p][p]+A[q][q]; 
             if (r == p && s != q) BP[i][j]=A[q][s];
             if (s == p) BP[i][j]=-A[q][r];  
           }
         }
       }
     }
     for (i=0;i<m;i++) {
       BP[m][i]=V1[i];
       BP[i][m]=W1[i];
     }
     for (i=0;i<m+2;i++) {
       BP[m+1][i]=V2[i];
       BP[i][m+1]=W2[i];
     }
     staFunc->CopyMatrix(BP,BT);
     Ret=staFunc->Decomp(BT);
     if (Ret != 0) {
       return 1;
     }
     staFunc->ClearVector(VB1);
     VB1[m]=1.0;
     staFunc->Solve(BT,VB1);
     staFunc->ClearVector(VB2);
     VB2[m+1]=1.0;
     staFunc->Solve(BT,VB2);                         
   /* find pair in (g11,g12), (g21,g22) with max absolute value */
     max1=(fabs(VB1[m])>fabs(VB2[m]))?VB1[m]:VB2[m];
     max2=(fabs(VB1[m+1])>fabs(VB2[m+1]))?VB1[m+1]:VB2[m+1];
     if (fabs(max1)>fabs(max2))
      imax=m;
     else imax=m+1;    
     for(i=0;i<m;i++) {
       Q[i]=VB2[imax]*VB1[i]-VB1[imax]*VB2[i];
     } 
  } 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);

    if (stagenData->BifDataInOk) {
/* Message(MSG_OK,"Eigenvectors exist"); 			*/
/* bifurcation data exist 					*/
/*  BifVector[0] to [nphase-1]:        Hopf continuation vector	*/
/*  BifVector[nphase] to [2*nphase-1]: global vector L 		*/
/*  BifVector[2*nphase]:	       -lambda*lambda   	*/

      MemCopy(V,stagenData->BifDataInPtr,sizeof(FloatHi)*nphase);
      MemCopy(L,stagenData->BifDataInPtr+nphase,sizeof(FloatHi)*nphase);
 
    } else {
     /* bifurcation data to be computed */
/*  Message(MSG_OK,"Eigenvectors have to be computed"); 	*/
      Ret=0;
{int i,j,imin,jmin;
 double r,s,d,phi,beta,gamma,lambda,omega;
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    /* find the pair of eigenvalues with zero sum */
    s=fabs(Re[0]+Re[1])+fabs(Im[0]+Im[1]);
    imin=0;
    jmin=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) {
          imin=i;
          jmin=j;
          s=d;
        }
      }
    }
    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) {
        staUtil->ErrMessage("Unable to find (+)eigenvector at the initial point");
        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) {
        staUtil->ErrMessage("Unable to find (-)eigenvector at the initial point");
        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);
    } else {
/* Hopf */
      lambda=Re[imin];
      omega=fabs(Im[imin]);
      /* 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) {
        staUtil->ErrMessage("Unable to find eigenvectors at the initial point");
        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=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]=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);
    }
polufinal:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
    }
}
  if (JacRHS(x0,C1)) return 1;  /* for safety */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  
/* Fill in the bordered bi-product matrix BP */
  staFunc->ClearMatrix(BP);
  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) BP[i][j]=-A[p][s];       
          if (r != p && s == q) BP[i][j]=A[p][r];   
          if (r == p && s == q) BP[i][j]=A[p][p]+A[q][q]; 
          if (r == p && s != q) BP[i][j]=A[q][s];
          if (s == p) BP[i][j]=-A[q][r];  
        }
      }
    }
  }
/* searching 2 rows with maximal pivots */  
  if (m >= 2) {
    G=staFunc->CreateMatrix(m+1,m+1);
    staFunc->CopyMatrixBlock(BP,0,0,G,0,0,m,m);
    for (i=0; i<m; i++) {
      G[i][m]=i;
      G[m][i]=i;
    } 
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<m-2;k++) {
      t=0.0;
      for (i=k;i<m;i++) {
	for (j=k;j<m;j++) {
    	  if (fabs(G[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(G[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=m;i++) {
          t=G[i][k];
          G[i][k]=G[i][q];
          G[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=m;i++) {
          t=G[k][i];         
          G[k][i]=G[p][i];
          G[p][i]=t;
	}
      }        
      if (fabs(G[k][k]) == 0.0) {
        staUtil->ErrMessage("Biproduct has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<m;i++) { /* perform elimination */
        t=-G[i][k]/G[k][k];
        for (j=kp1;j<m;j++) {
          G[i][j]+=t*G[k][j];
	}
      }
    } 
    p=(Int2)G[m-2][m];
    q=(Int2)G[m-1][m];   
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    W1[p]=1.0;
    W2[q]=1.0;
    BP[p][m]=1.0;
    BP[q][m+1]=1.0;
    p=(Int2)G[m][m-2];
    q=(Int2)G[m][m-1];   
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    V1[p]=1.0;
    V2[q]=1.0;
    BP[m][p]=1.0;
    BP[m+1][q]=1.0;    
    G=staFunc->DestroyMatrix(G);
      
/* solve first system to obtain new values for V1 */
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    staFunc->CopyMatrix(BP,BT);  
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(V1,K1);
    if (alpha < 0) {
      for (i=0;i<m+2;i++)
        K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,V1);
  
/* solve second system to obtain new values for W1 */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(W1,K1);
    if (alpha < 0) {
     for (i=0;i<m+2;i++)
       K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,W1);
/* adapt V2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }    
    staFunc->CopyVectorElements(W2,0,RL,0,m+1);
    staFunc->Decomp(BP1);
    staFunc->Solve(BP1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(V2,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(V2);
    staFunc->CopyVectorElements(RL,0,V2,0,m+1);
/* adapt W2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }      
    staFunc->TransposeMatrix(BP1,BT1);
    staFunc->Decomp(BT1);
    staFunc->ClearVector(RL);
    staFunc->CopyVectorElements(V2,0,RL,0,m+1);
    staFunc->Solve(BT1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(W2,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(W2);
    staFunc->CopyVectorElements(RL,0,W2,0,m+1);
  
    for (i=0;i<m;i++) {
      BP[m][i]=V1[i];
      BP[i][m]=W1[i];
    }
    for (i=0;i<m+2;i++) {
      BP[m+1][i]=V2[i];
      BP[i][m+1]=W2[i];
    }  
  } else {
     BP[0][m] = 1.0;
     BP[1][0] = 1.0;
     BP[m+1][m+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[m+1]=1.0;
    } 
  /* Compute a tangent vector to the Hopf curve */
    h1_DefHopf(x0,v0);
    Ret=h1_JacDefHopf(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 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;
}

/********************/
/* zh -> h1 Starter */
Entry(Int2) zh_h1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,Ret,m,p,q,r,s;
  Vector v1,v2;
  Matrix D,DE,G;
  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) {       /* "Continuation" */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(m+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(m+2),sizeof(FloatHi)*(m+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(m+2),sizeof(FloatHi)*(m+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(m+2),sizeof(FloatHi)*(m+2));    
    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);

    if (stagenData->BifDataInOk) {
    /* bifurcation data exist */
/*  BifVector[0] to [nphase-1]:          null-vector*/
/*  BifVector[nphase] to [2*nphase-1]:   Hopf continuation vector*/
/*  BifVector[2*nphase] to [3*nphase-1]: projection vector L */
/*  BifVector[3*nphase]:                 -lambda*lambda  */

      MemCopy(V,stagenData->BifDataInPtr+nphase,sizeof(FloatHi)*nphase);
      MemCopy(L,stagenData->BifDataInPtr+2*nphase,sizeof(FloatHi)*nphase);

    } else {
     /* bifurcation data to be computed */
      Ret=0;
{int i,j,imin,jmin;
 double r,s,d,phi,beta,gamma,lambda,omega;
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    /* find the pair of eigenvalues with zero sum */
    s=fabs(Re[0]+Re[1])+fabs(Im[0]+Im[1]);
    imin=0;
    jmin=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) {
          imin=i;
          jmin=j;
          s=d;
        }
      }
    }
    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) {
        staUtil->ErrMessage("Unable to find (+)eigenvector at the initial point");
        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) {
        staUtil->ErrMessage("Unable to find (-)eigenvector at the initial point");
        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);      
    } else {
/* Hopf */
      lambda=Re[imin];
      omega=fabs(Im[imin]);
      /* 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) {
        staUtil->ErrMessage("Unable to find eigenvectors at the initial point");
        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=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]=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);
    }
polufinal:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
    }
}

  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  
/* Fill in the bordered bi-product matrix BP */
  staFunc->ClearMatrix(BP);
  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) BP[i][j]=-A[p][s];       
          if (r != p && s == q) BP[i][j]=A[p][r];   
          if (r == p && s == q) BP[i][j]=A[p][p]+A[q][q]; 
          if (r == p && s != q) BP[i][j]=A[q][s];
          if (s == p) BP[i][j]=-A[q][r];  
        }
      }
    }
  } 
  if (m >= 2) {
    G=staFunc->CreateMatrix(m+1,m+1);
    staFunc->CopyMatrixBlock(BP,0,0,G,0,0,m,m);
    for (i=0; i<m; i++) {
      G[i][m]=i;
      G[m][i]=i;
    } 
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<m-2;k++) {
      t=0.0;
      for (i=k;i<m;i++) {
	for (j=k;j<m;j++) {
    	  if (fabs(G[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(G[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=m;i++) {
          t=G[i][k];
          G[i][k]=G[i][q];
          G[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=m;i++) {
          t=G[k][i];         
          G[k][i]=G[p][i];
          G[p][i]=t;
	}
      }        
      if (fabs(G[k][k]) == 0.0) {
        staUtil->ErrMessage("Biproduct has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<m;i++) { /* perform elimination */
        t=-G[i][k]/G[k][k];
        for (j=kp1;j<m;j++) {
          G[i][j]+=t*G[k][j];
	}
      }
    } 
    p=(Int2)G[m-2][m];
    q=(Int2)G[m-1][m];   
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    W1[p]=1.0;
    W2[q]=1.0;
    BP[p][m]=1.0;
    BP[q][m+1]=1.0;
    p=(Int2)G[m][m-2];
    q=(Int2)G[m][m-1];   
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    V1[p]=1.0;
    V2[q]=1.0;
    BP[m][p]=1.0;
    BP[m+1][q]=1.0;    
    G=staFunc->DestroyMatrix(G);
      
/* solve first system to obtain new values for V1 */
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    staFunc->CopyMatrix(BP,BT);  
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(V1,K1);
    if (alpha < 0) {
      for (i=0;i<m+2;i++)
        K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,V1);
  
/* solve second system to obtain new values for W1 */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(W1,K1);
    if (alpha < 0) {
     for (i=0;i<m+2;i++)
       K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,W1);
/* adapt V2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }    
    staFunc->CopyVectorElements(W2,0,RL,0,m+1);
    staFunc->Decomp(BP1);
    staFunc->Solve(BP1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(V2,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(V2);
    staFunc->CopyVectorElements(RL,0,V2,0,m+1);
/* adapt W2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }      
    staFunc->TransposeMatrix(BP1,BT1);
    staFunc->Decomp(BT1);
    staFunc->ClearVector(RL);
    staFunc->CopyVectorElements(V2,0,RL,0,m+1);
    staFunc->Solve(BT1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(W2,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(W2);
    staFunc->CopyVectorElements(RL,0,W2,0,m+1);
  
    for (i=0;i<m;i++) {
      BP[m][i]=V1[i];
      BP[i][m]=W1[i];
    }
    for (i=0;i<m+2;i++) {
      BP[m+1][i]=V2[i];
      BP[i][m+1]=W2[i];
    }   
  } else {
     BP[0][m] = 1.0;
     BP[1][0] = 1.0;
     BP[m+1][m+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[m+1]=1.0;  
    } 
  /* Compute a tangent vector to the Hopf curve */
    h1_DefHopf(x0,v0);
    Ret=h1_JacDefHopf(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 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;
}

/********************/
/* bt -> h1 Starter */
Entry(Int2) bt_h1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,Ret,m,p,q,r,s;
  Vector v1,v2;
  Matrix D,DE,G;
  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(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(m+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(m+2),sizeof(FloatHi)*(m+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(m+2),sizeof(FloatHi)*(m+2));
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(m+2),sizeof(FloatHi)*(m+2));
    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);

    if (stagenData->BifDataInOk) {
    /* bifurcation data exist */
/*  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(V,stagenData->BifDataInPtr+nphase,sizeof(FloatHi)*nphase);
      MemCopy(L,stagenData->BifDataInPtr,sizeof(FloatHi)*nphase);
      staFunc->NormalizeVector(V);
    } else {
     /* bifurcation data to be computed */
      Ret=0;
      staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(AA,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(AA,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) BB[i+nphase][i]=-1;
      i=staFunc->SngSolve(BB,_eps,VV);
      if (i != 2*nphase-2) {
        staUtil->ErrMessage("Unable to find (generalized)eigenvectors at the initial point");
        Ret=1;
        goto final;
      }
      staFunc->CopyVectorElements(VV,0,L,0,nphase);       /* null-vector */
      staFunc->NormalizeVector(L);
      staFunc->CopyVectorElements(VV,nphase,V,0,nphase);  /* generalized null-vector */
      staFunc->AddScaledVector(V,L,-staFunc->ScalProd(L,V),V);
      staFunc->NormalizeVector(V);
final:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
    }

  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
/* Fill in the bordered bi-product matrix BP */
  staFunc->ClearMatrix(BP);
  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) BP[i][j]=-A[p][s];       
          if (r != p && s == q) BP[i][j]=A[p][r];   
          if (r == p && s == q) BP[i][j]=A[p][p]+A[q][q]; 
          if (r == p && s != q) BP[i][j]=A[q][s];
          if (s == p) BP[i][j]=-A[q][r];  
        }
      }
    }
  }
/* new version for bordering */ 
  if (m >= 2) {
    G=staFunc->CreateMatrix(m+1,m+1);
    staFunc->CopyMatrixBlock(BP,0,0,G,0,0,m,m);
    for (i=0; i<m; i++) {
      G[i][m]=i;
      G[m][i]=i;
    } 
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<m-2;k++) {
      t=0.0;
      for (i=k;i<m;i++) {
	for (j=k;j<m;j++) {
    	  if (fabs(G[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(G[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=m;i++) {
          t=G[i][k];
          G[i][k]=G[i][q];
          G[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=m;i++) {
          t=G[k][i];         
          G[k][i]=G[p][i];
          G[p][i]=t;
	}
      }        
      if (fabs(G[k][k]) == 0.0) {
        staUtil->ErrMessage("Biproduct has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<m;i++) { /* perform elimination */
        t=-G[i][k]/G[k][k];
        for (j=kp1;j<m;j++) {
          G[i][j]+=t*G[k][j];
	}
      }
    } 
    p=(Int2)G[m-2][m];
    q=(Int2)G[m-1][m];   
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    W1[p]=1.0;
    W2[q]=1.0;
    BP[p][m]=1.0;
    BP[q][m+1]=1.0;
    p=(Int2)G[m][m-2];
    q=(Int2)G[m][m-1];   
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    V1[p]=1.0;
    V2[q]=1.0;
    BP[m][p]=1.0;
    BP[m+1][q]=1.0;    
    G=staFunc->DestroyMatrix(G);
      
/* solve first system to obtain new values for V1 */
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    staFunc->CopyMatrix(BP,BT);  
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(V1,K1);
    if (alpha < 0) {
      for (i=0;i<m+2;i++)
        K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,V1);
  
/* solve second system to obtain new values for W1 */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(W1,K1);
    if (alpha < 0) {
     for (i=0;i<m+2;i++)
       K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,W1);
/* adapt V2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }    
    staFunc->CopyVectorElements(W2,0,RL,0,m+1);
    staFunc->Decomp(BP1);
    staFunc->Solve(BP1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(V2,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(V2);
    staFunc->CopyVectorElements(RL,0,V2,0,m+1);
/* adapt W2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }      
    staFunc->TransposeMatrix(BP1,BT1);
    staFunc->Decomp(BT1);
    staFunc->ClearVector(RL);
    staFunc->CopyVectorElements(V2,0,RL,0,m+1);
    staFunc->Solve(BT1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(W2,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(W2);
    staFunc->CopyVectorElements(RL,0,W2,0,m+1);
  
    for (i=0;i<m;i++) {
      BP[m][i]=V1[i];
      BP[i][m]=W1[i];
    }
    for (i=0;i<m+2;i++) {
      BP[m+1][i]=V2[i];
      BP[i][m+1]=W2[i];
    } 
  } else {
     BP[0][m] = 1.0;
     BP[1][0] = 1.0;
     BP[m+1][m+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[m+1]=1.0;  
    } 

  /* Compute a tangent vector to the Hopf curve */
    h1_DefHopf(x0,v0);
    Ret=h1_JacDefHopf(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 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;
}

/*************************************************************/
/* 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 */
    CopyToX(x);
    CopyToP(x);
    H=stagenData->Der2(X,P,&zero);
    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) JacDefEquilib (Vector x, Matrix Jac);

/*************************************************/
/* Jacobian of defining functions for Hopf curve */  
Local(Int2) h_JacDefHopf (Vector x, Matrix C3, Matrix Jac) {
  int i;
  staFunc->ClearMatrix(Jac);
  if (JacDefEquilib(x,C)) return 1;
  staFunc->CopyMatrixBlock(C,0,0,Jac,0,0,nphase+nappend,nphase+2+nappend);
  if (JacDefEig(x,C3)) return 2;
  staFunc->CopyMatrixBlock(C3,0,0,Jac,nphase+nappend,0,nphase,2*nphase+3+nappend);
  for (i=0; i<nphase; i++) {
    Jac[2*nphase+nappend][i+nphase+2+nappend]=2.0*x[i+nphase+2+nappend];
    Jac[2*nphase+1+nappend][i+nphase+2+nappend]=L[i];
  }
  return 0;
}

/*******************/
/* dh -> h1 Starter */
Entry(Int2) dh_h1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,imin,jmin,imin0,jmin0,imin1,jmin1,Ret,m,p,q,r,s;
  VAR u,t,d,phi,beta,gamma,lambda,omega;
  VAR max1,max2, alpha;
  Vector xx0,v1,v2;
  Matrix C3,D,DE,G;

  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) {       /* "Continuation" */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(m+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(m+2),sizeof(FloatHi)*(m+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(m+2),sizeof(FloatHi)*(m+2));
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(m+2),sizeof(FloatHi)*(m+2));
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
  } else {
    xx0=staFunc->CreateVector(2*nphase+3+nappend);
    C3=staFunc->CreateMatrix(nphase,2*nphase+3+nappend);

    v1=staFunc->CreateVector(2*nphase+3+nappend);
    v2=staFunc->CreateVector(2*nphase+3+nappend);
    D=staFunc->CreateMatrix(2*nphase+2+nappend,2*nphase+3+nappend);
    DE=staFunc->CreateMatrix(2*nphase+3+nappend,2*nphase+3+nappend);

    MemCopy(xx0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) xx0[nphase+i]=staData->P[active[i]];
    JacRHS(xx0,C1);

    if (stagenData->BifDataInOk) {
    /* bifurcation data exist */
/* 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]  */

      if (_branch == 0) {
        /* primary branch */
        MemCopy(xx0+nphase+2+nappend,stagenData->BifDataInPtr,sizeof(FloatHi)*nphase);
        MemCopy(V,stagenData->BifDataInPtr,sizeof(FloatHi)*nphase);
        MemCopy(L,stagenData->BifDataInPtr+nphase,sizeof(FloatHi)*nphase);
        xx0[2*nphase+2+nappend]=stagenData->BifDataInPtr[2*nphase];
      } else {
        /* secondary branch */
        MemCopy(xx0+nphase+2+nappend,stagenData->BifDataInPtr+2*nphase+1,sizeof(FloatHi)*nphase);
        MemCopy(V,stagenData->BifDataInPtr+2*nphase+1,sizeof(FloatHi)*nphase);
        MemCopy(L,stagenData->BifDataInPtr+3*nphase+1,sizeof(FloatHi)*nphase);
        xx0[2*nphase+2+nappend]=stagenData->BifDataInPtr[4*nphase+1];
      }
    } else {
     /* bifurcation data to be computed */
      Ret=0;
      staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
      staFunc->EigenValues(AA,Re,Im);
      /* find the pair of eigenvalues with zero sum */
      u=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<u) {
            imin0=i;
            jmin0=j;
            u=d;
          }
        }
      }
      imin1=imin0+1;
      if (imin1==nphase) imin1=0; 
      jmin1=jmin0+1;
      if (jmin1==nphase) jmin1=0; 
      u=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<u) {
                imin1=i;
                jmin1=j;
                u=d;
              }
	    }
	  }
        }
      }
/* printf("imin0=%i, jmin0=%i, imin1=%i, jmin1=%i\n",imin0,jmin0,imin1,jmin1); */
      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) {
          staUtil->ErrMessage("Unable to find (+)eigenvector at the initial point");
          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) {
          staUtil->ErrMessage("Unable to find (-)eigenvector at the initial point");
          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);
        staFunc->CopyVectorElements(V,0,xx0,nphase+2,nphase);
        xx0[2*nphase+2]=-lambda*lambda;
      } else {
/* Hopf */
        lambda=Re[imin];
        omega=fabs(Im[imin]);
        /* 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) {
          staUtil->ErrMessage("Unable to find eigenvectors at the initial point");
          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);
        u=staFunc->ScalProd(W,W);
        t=staFunc->ScalProd(V,W);
        if (t!=0) {
          beta=sqrt(fabs(4.0*t*t+(u-d)*(u-d)));
          if (2.0*fabs(2.0*t) > beta) {
            phi=PID2-atan(fabs((u-d)/(2.0*t)));
          } else {
            phi=atan(fabs(2.0*r/(u-d)));
          }
          if (t < 0) phi=PI-phi;
          if ((u-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,xx0,nphase+2+nappend,nphase);
        staFunc->CopyVectorElements(VV,nphase,L,0,nphase);
        xx0[2*nphase+2+nappend]=omega*omega;
        staFunc->NormalizeVector(L);
      }
polufinal:
      if (Ret) {
        xx0=staFunc->DestroyVector(x0);
	C3=staFunc->DestroyMatrix(C3);
        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 Hopf curve by projecting the standard augmented
    tangent vector (see __h.c) */
    Ret=h_JacDefHopf(xx0,C3,D);
    for (i=0; i<2*nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[2*nphase+2+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Hopf curve");   
       xx0=staFunc->DestroyVector(x0);
       C3=staFunc->DestroyMatrix(C3);
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }
   /* project the point and the tangent vector */
    staFunc->CopyVectorElements(xx0,0,x0,0,nphase+2+nappend);    
    staFunc->CopyVectorElements(v2,0,v0,0,nphase+2+nappend);    
    staFunc->NormalizeVector(v0);

    xx0=staFunc->DestroyVector(xx0);
    C3=staFunc->DestroyMatrix(C3);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);

  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  
/* Fill in the bordered bi-product matrix BP */
  staFunc->ClearMatrix(BP);
  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) BP[i][j]=-A[p][s];       
          if (r != p && s == q) BP[i][j]=A[p][r];   
          if (r == p && s == q) BP[i][j]=A[p][p]+A[q][q]; 
          if (r == p && s != q) BP[i][j]=A[q][s];
          if (s == p) BP[i][j]=-A[q][r];  
        }
      }
    }
  } 
/* searching 2 rows with maximal pivots */  
  if (m >= 2) {
    G=staFunc->CreateMatrix(m+1,m+1);
    staFunc->CopyMatrixBlock(BP,0,0,G,0,0,m,m);
    for (i=0; i<m; i++) {
      G[i][m]=i;
      G[m][i]=i;
    } 
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<m-2;k++) {
      t=0.0;
      for (i=k;i<m;i++) {
	for (j=k;j<m;j++) {
    	  if (fabs(G[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(G[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=m;i++) {
          t=G[i][k];
          G[i][k]=G[i][q];
          G[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=m;i++) {
          t=G[k][i];         
          G[k][i]=G[p][i];
          G[p][i]=t;
	}
      }        
      if (fabs(G[k][k]) == 0.0) {
        staUtil->ErrMessage("Biproduct has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<m;i++) { /* perform elimination */
        t=-G[i][k]/G[k][k];
        for (j=kp1;j<m;j++) {
          G[i][j]+=t*G[k][j];
	}
      }
    } 
    p=(Int2)G[m-2][m];
    q=(Int2)G[m-1][m];   
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    W1[p]=1.0;
    W2[q]=1.0;
    BP[p][m]=1.0;
    BP[q][m+1]=1.0;
    p=(Int2)G[m][m-2];
    q=(Int2)G[m][m-1];   
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    V1[p]=1.0;
    V2[q]=1.0;
    BP[m][p]=1.0;
    BP[m+1][q]=1.0;    
    G=staFunc->DestroyMatrix(G);
      
/* solve first system to obtain new values for V1 */
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    staFunc->CopyMatrix(BP,BT);  
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(V1,K1);
    if (alpha < 0) {
      for (i=0;i<m+2;i++)
        K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,V1);
  
/* solve second system to obtain new values for W1 */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(W1,K1);
    if (alpha < 0) {
     for (i=0;i<m+2;i++)
       K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,W1);
/* adapt V2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }    
    staFunc->CopyVectorElements(W2,0,RL,0,m+1);
    staFunc->Decomp(BP1);
    staFunc->Solve(BP1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(V2,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(V2);
    staFunc->CopyVectorElements(RL,0,V2,0,m+1);
/* adapt W2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }      
    staFunc->TransposeMatrix(BP1,BT1);
    staFunc->Decomp(BT1);
    staFunc->ClearVector(RL);
    staFunc->CopyVectorElements(V2,0,RL,0,m+1);
    staFunc->Solve(BT1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(W2,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(W2);
    staFunc->CopyVectorElements(RL,0,W2,0,m+1);
  
    for (i=0;i<m;i++) {
      BP[m][i]=V1[i];
      BP[i][m]=W1[i];
    }
    for (i=0;i<m+2;i++) {
      BP[m+1][i]=V2[i];
      BP[i][m+1]=W2[i];
    }      
  } else {
     BP[0][m] = 1.0;
     BP[1][0] = 1.0;
     BP[m+1][m+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[m+1]=1.0;
    }
  } 
  /* 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 */
/*  printf("dh_h1_Starter left\n"); */
  return 0;
}

/************************/
/* 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) CopyToE(Vector Re, Vector Im) {
  MemCopy(E,Re,sizeof(FloatHi)*nphase);
  MemCopy(E+nphase,Im,sizeof(FloatHi)*nphase);
}


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

Entry(Int2) h1_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) h1_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) h1_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) {
  h1_DefRHS(x,y);
  if (nappend) h1_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(h1_DefRHS, nphase, xp, dx, Jac)) {
      return 1; }
  }
  return 0;
}
/***************************************************/
/* Jacobian of defining conditions for equilibria  */
 Local(Int2) JacDefEquilib (Vector x, Matrix Jac) {
  /* Message(MSG_OK,"JacDefEquilib entered"); */
  staFunc->ClearMatrix(Jac);
  JacRHS(x, Jac);
  if (nappend) staFunc->Der(h1_DefUser, nappend, x, dx, Jac+nphase);
  /* Message(MSG_OK,"JacDefEquilib left"); */
  return 0;
}
/***************************/
/* Hopf defining condition */
Entry(Int2) h1_DefBip (Vector x, Vector y) {
Int2 i,j,p,q,r,s,m,imax;
Int2 Ret;
VAR max1,max2;
  m=nphase*(nphase-1)/2;
   
  if (JacRHS(x,C1)) return 1;	  /* A,C1,BP,V1,V2,W1,W2 - global */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

/* Fill in the bordered bi-product matrix BP */
  staFunc->ClearMatrix(BP);
  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) BP[i][j]=-A[p][s];       
          if (r != p && s == q) BP[i][j]=A[p][r];   
          if (r == p && s == q) BP[i][j]=A[p][p]+A[q][q]; 
          if (r == p && s != q) BP[i][j]=A[q][s];
          if (s == p) BP[i][j]=-A[q][r];  
        }
      }
    }
  }
  for (i=0; i<m; i++) {
    BP[m][i]=V1[i];
    BP[i][m]=W1[i];
  }
  for (i=0; i<m+2; i++) {  
    BP[m+1][i]=V2[i];
    BP[i][m+1]=W2[i];
  }

/* compute Hopf defining function */
  staFunc->CopyMatrix(BP,BT);
  Ret=staFunc->Decomp(BT);
  if (Ret != 0) {
    return 1;
  }
  staFunc->ClearVector(VB1);
  VB1[m]=1.0;
  staFunc->Solve(BT,VB1);
  staFunc->ClearVector(VB2);
  VB2[m+1]=1.0;
  staFunc->Solve(BT,VB2);                         
  y[0]=(VB1[m]*VB2[m+1]-VB1[m+1]*VB2[m]); 
/* find pair in (g11,g12), (g21,g22) with max absolute value */
  max1=(fabs(VB1[m])>fabs(VB2[m]))?VB1[m]:VB2[m];
  max2=(fabs(VB1[m+1])>fabs(VB2[m+1]))?VB1[m+1]:VB2[m+1];
  if (fabs(max1)>fabs(max2))
   imax=m;
  else imax=m+1;    
  for(i=0;i<m;i++) {
    Q[i]=VB2[imax]*VB1[i]-VB1[imax]*VB2[i];
  }
  /* Message(MSG_OK,"DefBip left"); */
  return 0;
}

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

  m=nphase*(nphase-1)/2;
  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2) {
/* creation of DefBip Jacobian matrix Jac using symbolic second derivatives */
    CopyToX(x);
    CopyToP(x);
    H=stagenData->Der2(X,P,&zero);
/* take global vector VB computed at the computation of DefBip */
/* compute PP1 and PP2 using the transposed bordered bi-product */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->ClearVector(PP1);
    PP1[m]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[m+1]=1.0;
    staFunc->Solve(BT,PP2);

/* compute gradient of DefBip */
    for (l=0; l<nphase+2+nappend; l++) {
      staFunc->ClearMatrix(BT);
      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) BT[i][j]=-HessianElem(p,s,l);       
              if (r != p && s == q) BT[i][j]=HessianElem(p,r,l);   
              if (r == p && s == q) BT[i][j]=HessianElem(p,p,l)+HessianElem(q,q,l); 
              if (r == p && s != q) BT[i][j]=HessianElem(q,s,l);
              if (s == p) BT[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-=PP1[i]*VB1[j]*BT[i][j];
	  s12-=PP1[i]*VB2[j]*BT[i][j];
	  s21-=PP2[i]*VB1[j]*BT[i][j];
	  s22-=PP2[i]*VB2[j]*BT[i][j];
	}
      }
       Jac[0][l]=VB1[m]*s22+VB2[m+1]*s11-VB2[m]*s21-VB1[m+1]*s12; 
    }
  } else {
    if (staFunc->Der(h1_DefBip, 1, x, dx, Jac)) return 1;
  }
  return(0);
}

/*************************************/
/* Defining functions for Hopf curve */

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

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

  return 0;
}

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

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

  /* Assert: v1==NULL, msg==NULL, Ind==0 */
  h1_DefHopf (x, xp1);  /* compute global Q for test functions */
  MemCopy(x1,x,sizeof(FloatHi)*(nphase+2+nappend));     /* for the last point */
  if (!staData->eigcomp) return 0;
  if (JacRHS(x,C1)) return (1); 

  staFunc->EigenValues(C1,Re,Im);
  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[80];

/* first Lyapunov coefficient */
Entry(Int2) h1_TestBautin(Vector x, Vector v, FloatHiPtr res) {
VAR beta,gamma,phi,omega,s,r,h;
VAR d,a11,a12,a21,a22;
Int2 i,j,k,l,Ret;  
Int2 imax,jmax;

/* reevaluate Jacobian */
   if (JacRHS (x,C1)) {
     return 1;
   }
   staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
/* it is assumed that Q is correct */
/* find a component of Q with maxfabs */
  imax=1;
  jmax=0;
  d=fabs(Q[0]);
  for (i=0; i<nphase; i++) {
    for (j=0; j<i; j++) {
      if (fabs(Q[i*(i-1)/2+j])>d) {
        imax=i;
	jmax=j;
	d=fabs(Q[i*(i-1)/2+j]);
      }
    }
  }
  d=Q[imax*(imax-1)/2+jmax];

/* compute eigenvectors */
  staFunc->ClearVector(V);
  staFunc->ClearVector(W);
  for (i=0; i<nphase; i++) {
    if (i<imax) V[i]=Q[imax*(imax-1)/2+i]/d;
    if (i>imax) V[i]=-Q[i*(i-1)/2+imax]/d;
    if (i<jmax) W[i]=-Q[jmax*(jmax-1)/2+i]/d;
    if (i>jmax) W[i]=Q[i*(i-1)/2+jmax]/d;
  }

/* compute -lambda*lambda */
  staFunc->MultiplyMatrixVector(A,V,F);  /* F=A*V  */
  a11=staFunc->ScalProd(V,F);
  a12=staFunc->ScalProd(W,F);
  staFunc->MultiplyMatrixVector(A,W,F);  /* F=A*W  */
  a21=staFunc->ScalProd(W,F);
  a22=staFunc->ScalProd(V,F);
  d=staFunc->ScalProd(V,V)*staFunc->ScalProd(W,W)
   -staFunc->ScalProd(V,W)*staFunc->ScalProd(V,W);
  d=(a11*a21-a22*a12)/d;

  if (d <= 0.0) { 	/* neutral saddle */
    *res=777;
  } else {       	/* Hopf */
    omega=sqrt(d);
    staFunc->NormalizeVector(V);
    staFunc->MultiplyMatrixVector(A,V,W);
    staFunc->MultiplyVectorScalar(W,-1.0/omega,W);

    /* normalize real and imaginary parts <V,V>=1, <V,W>=0 */
    d=staFunc->ScalProd(V,V);
    s=staFunc->ScalProd(W,W);
    r=staFunc->ScalProd(V,W);
    if (r==0) {
      staFunc->CopyVectorElements(V,0,VV,0,nphase);
      staFunc->CopyVectorElements(W,0,VV,nphase,nphase);
    } else {
      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]=V[i]*cos(phi)-W[i]*sin(phi);
        VV[i+nphase]=V[i]*sin(phi)+W[i]*cos(phi);
      }
    }
    staFunc->CopyVectorElements(VV,0,W,0,nphase);
    d=staFunc->ScalProd(W,W);
    gamma=1/sqrt(d);
    staFunc->MultiplyVectorScalar(VV,gamma,VV);
    staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);
    
/* compute real and imaginary part of the adjoint eigenvector */
    staFunc->CopyMatrixBlock(A,0,0,AA,0,0,nphase,nphase);
    staFunc->TransposeMatrix(AA,AT);
    staFunc->ClearMatrix(BB);
    staFunc->CopyMatrixBlock(AT,0,0,BB,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(AT,0,0,BB,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) {
       BB[i][i+nphase]=-omega;
       BB[i][i]-=0.0;
       BB[i+nphase][i]=omega;
       BB[i+nphase][i+nphase]-=0.0;
    }
    if (staFunc->SngSolve(BB,_eps,VV) != 2*nphase-2) {
       *res=1.0;
       test[0]=*res;
       return 1;
    }
    staFunc->CopyVectorElements(VV,0,Pr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Pi,0,nphase);

    /* normalize the adjoint vector <Pr,Qr>+<Pi,Qi>=1, <Pr,Qi>-<Pi,Qr>=0 */
    d=staFunc->ScalProd(Pr,Qr);
    s=staFunc->ScalProd(Pi,Qi);
    r=staFunc->ScalProd(Pr,Qi);
    h=staFunc->ScalProd(Pi,Qr);

/*  phi=atan((r-h)/(d+s));  */
    beta=sqrt(fabs((r-h)*(r-h)+(d+s)*(d+s)));
    if (2.0*fabs(d+s) > beta) {
        phi=atan(fabs((r-h)/(d+s)));
      } else {
        phi=PID2-atan(fabs((d+s)/(r-h)));
      }
    if ((r-h) < 0) phi=PI-phi;
    if ((d+s) < 0) phi=-phi;

    for (i=0; i<nphase; i++) {
       VV[i]=Pr[i]*cos(phi)-Pi[i]*sin(phi);
       VV[i+nphase]=Pr[i]*sin(phi)+Pi[i]*cos(phi);
    }
    staFunc->CopyVectorElements(VV,0,Pr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Pi,0,nphase);
    r=staFunc->ScalProd(Pr,Qr);
    h=staFunc->ScalProd(Pi,Qi);
    beta=1.0/(r+h);
    staFunc->MultiplyVectorScalar(Pr,beta,Pr);
    staFunc->MultiplyVectorScalar(Pi,beta,Pi);

/* compute the first Lyapunov coefficient L1 */
/* quadratic terms */
      if (stagenData->Der2)  {
        CopyToX(x);
        CopyToP(x);
        H=stagenData->Der2(X,P,&zero);	/* The 1st elem must be time! */
	for (i=0; i<nphase; i++) {
           for (s=0.0,r=0.0,h=0.0,j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) {
                 s+=HessianElem(i,j,k)*Qr[j]*Qr[k];
                 r+=HessianElem(i,j,k)*Qi[j]*Qi[k];
                 h+=HessianElem(i,j,k)*Qr[j]*Qi[k];
	      }
           }
	   aa[i]=s;	
	   bb[i]=r;	
	   cc[i]=-2.0*h;	
	}
      } else {
        staFunc->ClearVector(v0);
        staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
        h1_DefRHS(x,aa);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        h1_DefRHS(x1,bb);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        h1_DefRHS(x1,cc);
	staFunc->AddScaledVector(bb,aa,-2.0,rr);
	staFunc->AddScaledVector(rr,cc,1.0,aa);
	staFunc->MultiplyVectorScalar(aa,1.0/ddx/ddx,aa);

        staFunc->CopyVectorElements(Qi,0,v0,0,nphase);
        h1_DefRHS(x,rr);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        h1_DefRHS(x1,bb);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        h1_DefRHS(x1,cc);
	staFunc->AddScaledVector(bb,rr,-2.0,rr);
	staFunc->AddScaledVector(rr,cc,1.0,bb);
	staFunc->MultiplyVectorScalar(bb,1.0/ddx/ddx,bb);

        staFunc->AddScaledVector(Qr,Qi,1.0,rr);
        staFunc->CopyVectorElements(rr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,rr);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,cc);
	staFunc->AddScaledVector(rr,cc,1.0,dd);

        staFunc->AddScaledVector(Qr,Qi,-1.0,rr);
        staFunc->CopyVectorElements(rr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,rr);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,cc);
	staFunc->AddScaledVector(rr,cc,1.0,rr);
        staFunc->AddScaledVector(dd,rr,-1.0,cc);
	staFunc->MultiplyVectorScalar(cc,-2.0/ddx/ddx,cc);   /* cc=-2*c */
      }

      staFunc->CopyMatrixBlock(A,0,0,AA,0,0,nphase,nphase);
      staFunc->AddScaledVector(aa,bb,1.0,rr);
      Ret=staFunc->Decomp(AA);
      if (Ret != 0) {    /* Zero-Hopf case */
        *res=777.0;
	test[0]=*res;
	return 0;
      }
      staFunc->Solve(AA,rr);
      staFunc->MultiplyVectorScalar(rr,-2.0,rr);		/* rr=-2*r */
      
      staFunc->AddScaledVector(bb,aa,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,VV,0,nphase);
      staFunc->CopyVectorElements(cc,0,VV,nphase,nphase);
      staFunc->ClearMatrix(BB);
      staFunc->CopyMatrixBlock(A,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(A,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) {
        BB[i][i+nphase]=2.0*omega;
        BB[i][i]-=0.0;
        BB[i+nphase][i]=-2.0*omega;
        BB[i+nphase][i+nphase]-=0.0;
      }
      Ret=staFunc->Decomp(BB);
      if (Ret != 0) {
        *res=1.0;
	test[0]=*res;
	return 3;      
      }
      staFunc->Solve(BB,VV);
      staFunc->CopyVectorElements(VV,0,aa,0,nphase);          	/* aa=Sr */
      staFunc->CopyVectorElements(VV,nphase,bb,0,nphase);	/* bb=Si */

      if (stagenData->Der2)  {
        for (s=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) 
	         s+=HessianElem(i,j,k)*(
		   Pr[i]*(Qr[j]*rr[k]+Qr[j]*aa[k]+Qi[j]*bb[k])+
		   Pi[i]*(Qi[j]*rr[k]+Qr[j]*bb[k]-Qi[j]*aa[k]));
           }
        }
      } else {
        staFunc->ClearVector(v0);

/*  -2*<Pr,B(Qr,r)>  */
        staFunc->AddScaledVector(Qr,rr,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qr,rr,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
/*  -2*<Pi,B(Qi,r)>  */
        staFunc->AddScaledVector(Qi,rr,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qi,rr,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
/*  <Pr,B(Qr,Sr)>  */
        staFunc->AddScaledVector(Qr,aa,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qr,aa,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
/*  <Pr,B(Qi,Si)>  */
        staFunc->AddScaledVector(Qi,bb,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qi,bb,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
/*  <Pi,B(Qr,Si)>  */
        staFunc->AddScaledVector(Qr,bb,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qr,bb,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
/*  <Pi,B(Qr,Sr)>  */
        staFunc->AddScaledVector(Qi,aa,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qi,aa,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
	s/=(ddx*ddx); 
      }            
/* cubic terms */
     if (stagenData->Der3) {
        CopyToX(x);
        CopyToP(x);
        D3=stagenData->Der3(X,P,&zero);	/* The 1st elem must be time! */
        for (r=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++) {
	           r+=D3Elem(i,j,k,l)*(Pr[i]*Qr[j]*(Qr[k]*Qr[l]+Qi[k]*Qi[l])+
		                       Pi[i]*Qi[l]*(Qi[j]*Qi[k]+Qr[j]*Qr[k]));
	         }
	      }
           }
        }	
     } else {
        staFunc->ClearVector(v0);
	
	staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
	
	staFunc->CopyVectorElements(Qi,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
	
        staFunc->AddScaledVector(Qr,Qi,1.0,aa);
        staFunc->CopyVectorElements(aa,0,v0,0,nphase);
        staFunc->AddScaledVector(Pr,Pi,1.0,bb);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/6;

        staFunc->AddScaledVector(Qr,Qi,-1.0,aa);
        staFunc->CopyVectorElements(aa,0,v0,0,nphase);
        staFunc->AddScaledVector(Pr,Pi,-1.0,bb);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        h1_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        h1_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/6;
        r/=(dddx*dddx*dddx);                                         
     }
    *res=(s+r)/(2*omega);
  }
  test[0]=*res;
  return 0;
}
Entry(Int2) h1_ProcBautin(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
Int2 i,j,imax,jmax;
VAR d,a11,a12,a21,a22;
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Generalized Hopf (?)");  
  else {
/* reevaluate Jacobian */
   if (JacRHS (x,C1)) return 1;
   staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

/* it is assumed that Q is correct */
/* find a component of Q with maxfabs */
  imax=1;
  jmax=0;
  d=fabs(Q[0]);
  for (i=0; i<nphase; i++) {
    for (j=0; j<i; j++) {
      if (fabs(Q[i*(i-1)/2+j])>d) {
        imax=i;
	jmax=j;
	d=fabs(Q[i*(i-1)/2+j]);
      }
    }
  }
  d=Q[imax*(imax-1)/2+jmax];

/* compute eigenvectors */
  staFunc->ClearVector(V);
  staFunc->ClearVector(W);
  for (i=0; i<nphase; i++) {
    if (i<imax) V[i]=Q[imax*(imax-1)/2+i]/d;
    if (i>imax) V[i]=-Q[i*(i-1)/2+imax]/d;
    if (i<jmax) W[i]=-Q[jmax*(jmax-1)/2+i]/d;
    if (i>jmax) W[i]=Q[i*(i-1)/2+jmax]/d;
  }

/* compute -lambda*lambda */
  staFunc->MultiplyMatrixVector(A,V,F);  /* F=A*V  */
  a11=staFunc->ScalProd(V,F);
  a12=staFunc->ScalProd(W,F);
  staFunc->MultiplyMatrixVector(A,W,F);  /* F=A*W  */
  a21=staFunc->ScalProd(W,F);
  a22=staFunc->ScalProd(V,F);

  d=staFunc->ScalProd(V,V)*staFunc->ScalProd(W,W)
   -staFunc->ScalProd(V,W)*staFunc->ScalProd(V,W);
  
/* orthogonalize and normalize eigenvectors */
    staFunc->NormalizeVector(V);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),L);
    staFunc->NormalizeVector(L);
  
/* provide bifurcation data for gh -> h starter */
/*  BifVector[0] to [nphase-1]:  	 Hopf continuation vector*/
/*  BifVector[nphase] to [2*nphase-1]: 	 global vector L */
/*  BifVector[2*nphase]:                 -lambda*lambda  */

    MemCopy(BifVector,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase,L,sizeof(FloatHi)*nphase);
    BifVector[2*nphase]=(a11*a21-a22*a12)/d;
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1);
    sprintf(buf,"Generalized Hopf: omega=%g",sqrt(BifVector[2*nphase]));    
  }  
  *msg=buf;
  return 0;
}

/* -lambda^2  */
Entry(Int2) h1_TestBogdanovTakens(Vector x, Vector v, FloatHiPtr res) {
Int2 i,j,imax,jmax;
VAR d,a11,a12,a21,a22;
  *res=1;
  if (nphase == 1) goto end;

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

/* it is assumed that Q is correct */
/* find a component of Q with max(fabs) */
  imax=1;
  jmax=0;
  d=fabs(Q[0]);
  for (i=0; i<nphase; i++) {
    for (j=0; j<i; j++) {
      if (fabs(Q[i*(i-1)/2+j])>d) {
        imax=i;
	jmax=j;
	d=fabs(Q[i*(i-1)/2+j]);
      }
    }
  }
  d=Q[imax*(imax-1)/2+jmax];

/* compute eigenvectors */
  staFunc->ClearVector(V);
  staFunc->ClearVector(W);
  for (i=0; i<nphase; i++) {
    if (i<imax) V[i]=Q[imax*(imax-1)/2+i]/d;
    if (i>imax) V[i]=-Q[i*(i-1)/2+imax]/d;
    if (i<jmax) W[i]=-Q[jmax*(jmax-1)/2+i]/d;
    if (i>jmax) W[i]=Q[i*(i-1)/2+jmax]/d;
  }

/* compute -lambda*lambda */
  staFunc->MultiplyMatrixVector(A,V,F);  /* F=A*V  */
  a11=staFunc->ScalProd(V,F);
  a12=staFunc->ScalProd(W,F);
  staFunc->MultiplyMatrixVector(A,W,F);  /* F=A*W  */
  a21=staFunc->ScalProd(W,F);
  a22=staFunc->ScalProd(V,F);

  d=staFunc->ScalProd(V,V)*staFunc->ScalProd(W,W)
   -staFunc->ScalProd(V,W)*staFunc->ScalProd(V,W);
  
  *res=(a11*a21-a22*a12)/d;
end:
  test[1]=*res;  
  return 0;
}

Entry(Int2) h1_ProcBogdanovTakens(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Bogdanov-Takens (?)");
  else {
Vector v1,W0,W1;
Matrix D;
Int2 i,j,k;
double 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 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,"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);
        h1_DefRHS(x,F);
        f0=staFunc->ScalProd(W0,F);
        f3=staFunc->ScalProd(W1,F);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        h1_DefRHS(x1,F);
        f1=staFunc->ScalProd(W0,F);
        f4=staFunc->ScalProd(W1,F);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        h1_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);
        h1_DefRHS(x1,F);
	b11=staFunc->ScalProd(W1,F);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_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);
        h1_DefRHS(x1,F);
	b11-=staFunc->ScalProd(W1,F);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        h1_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;
}

/* det of the Jacobian matrix */
Entry(Int2) h1_TestZeroHopf(Vector x, Vector v, FloatHiPtr res) {
 /* Message(MSG_OK,"TestZeroHopf entered"); */
  *res=1.0;
  if (nphase < 2) goto end;
  if (JacRHS (x,C1)) return 1;
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  if (staFunc->Decomp(A)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(A,res,NULL);
  } 
end:
  test[2]=*res; 
  /* Message(MSG_OK,"TestZeroHopf left"); */
  return 0;
}

Entry(Int2) h1_ProcZeroHopf(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
Int2 i,j,imax,jmax;
VAR d,a11,a12,a21,a22;
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Zero-Hopf (?)"); 
  else {
/* reevaluate Jacobian */
   if (JacRHS (x,C1)) return 1;
   staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

/* it is assumed that Q is correct */
/* find a component of Q with maxfabs */
  imax=1;
  jmax=0;
  d=fabs(Q[0]);
  for (i=0; i<nphase; i++) {
    for (j=0; j<i; j++) {
      if (fabs(Q[i*(i-1)/2+j])>d) {
        imax=i;
	jmax=j;
	d=fabs(Q[i*(i-1)/2+j]);
      }
    }
  }
  d=Q[imax*(imax-1)/2+jmax];

/* compute eigenvectors */
  staFunc->ClearVector(V);
  staFunc->ClearVector(W);
  for (i=0; i<nphase; i++) {
    if (i<imax) V[i]=Q[imax*(imax-1)/2+i]/d;
    if (i>imax) V[i]=-Q[i*(i-1)/2+imax]/d;
    if (i<jmax) W[i]=-Q[jmax*(jmax-1)/2+i]/d;
    if (i>jmax) W[i]=Q[i*(i-1)/2+jmax]/d;
  }

/* compute -lambda*lambda */
  staFunc->MultiplyMatrixVector(A,V,F);  /* F=A*V  */
  a11=staFunc->ScalProd(V,F);
  a12=staFunc->ScalProd(W,F);
  staFunc->MultiplyMatrixVector(A,W,F);  /* F=A*W  */
  a21=staFunc->ScalProd(W,F);
  a22=staFunc->ScalProd(V,F);

  d=staFunc->ScalProd(V,V)*staFunc->ScalProd(W,W)
   -staFunc->ScalProd(V,W)*staFunc->ScalProd(V,W);
  
/* orthogonalize and normalize eigenvectors */
    staFunc->NormalizeVector(V);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),L);
    staFunc->NormalizeVector(L);
  
/* provide bifurcation data for zh -> h starter */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    if (staFunc->SngSolve(A,_eps,W) != nphase-1) {
      sprintf(buf,"Zero-Hopf: null-vector (?)"); 
      *msg=buf;
      return 0;
    }   
    staFunc->NormalizeVector(W);
/*  BifVector[0] to [nphase-1]:          null-vector */
/*  BifVector[nphase] to [2*nphase-1]:   Hopf continuation vector */
/*  BifVector[2*nphase] to [3*nphase-1]: global vector L */
/*  BifVector[3*nphase]:                 -lambda*lambda  */
    MemCopy(BifVector,W,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+2*nphase,L,sizeof(FloatHi)*nphase);
    BifVector[3*nphase]=(a11*a21-a22*a12)/d;
    stagenData->BifDataOutLen=sizeof(FloatHi)*(3*nphase+1);
    
    if (BifVector[3*nphase]>0) {
      sprintf(buf,"Zero-Hopf: omega=%g",sqrt(BifVector[3*nphase]));
    } else {
      sprintf(buf,"Zero eigenvalue + Neutral saddle: lambda=%g",sqrt(fabs(BifVector[3*nphase])));
    }      
  }  
  *msg=buf;
  return 0;
}

/*  DH test function */
Entry(Int2) h1_TestDoubleHopf(Vector x, Vector v, FloatHiPtr res) {
Int2 m;
  m=nphase*(nphase-1)/2;
  *res=VB2[m+1]; /* g22 = 0 */
  test[3]=*res;
  return 0;
}

Entry(Int2) h1_ProcDoubleHopf(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;
 double 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=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]=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;
}

/**************************************************/
/*  Adaptation of the bordering vectors V1,V2,W1,W2 */
Entry(Int2) h1_Adapter(Vector x, Vector v) {
int i,m;
VAR alpha,beta,max1,max2;
  m=nphase*(nphase-1)/2;
/* adapt only in case m >= 2 ; for m = 1 no adaptation is needed */
  if (m >= 2) {
/* solve first system to obtain new values for V1 */
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    staFunc->CopyMatrix(BP,BT);  
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(V1,K1);
    if (alpha < 0) {
      for (i=0;i<m+2;i++)
        K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,V1);
  
/* solve second system to obtain new values for W1 */
    staFunc->TransposeMatrix(BP,BT);
    staFunc->ClearVector(K1);
    staFunc->ClearVector(K2);
    K1[m]=1.0;
    K2[m+1]=1.0;
    staFunc->Decomp(BT);
    staFunc->Solve(BT,K1);
    staFunc->Solve(BT,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(W1,K1);
    if (alpha < 0) {
     for (i=0;i<m+2;i++)
       K1[i]=-K1[i];
    }
    staFunc->CopyVector(K1,W1);
/* adapt V2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }    
    staFunc->CopyVectorElements(W2,0,RL,0,m+1);
    staFunc->Decomp(BP1);
    staFunc->Solve(BP1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(V2,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(V2);
    staFunc->CopyVectorElements(RL,0,V2,0,m+1);
/* adapt W2 */
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,m,m);
    for (i=0;i<m;i++) {
     BP1[i][m] = W1[i];
     BP1[m][i] = V1[i];
    }
    staFunc->TransposeMatrix(BP1,BT1);
    staFunc->Decomp(BT1);
    staFunc->ClearVector(RL);
    staFunc->CopyVectorElements(V2,0,RL,0,m+1);
    staFunc->Solve(BT1,RL);
    staFunc->NormalizeVector(RL);
    staFunc->CopyVectorElements(W2,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(W2);
    staFunc->CopyVectorElements(RL,0,W2,0,m+1);
  } /* if m >=2 */
/* adapt test function for DH if set */
  if (staData->ActiveSing[3]) { 
   for (i=0;i<m;i++) {
      BP[m][i]=V1[i];
      BP[i][m]=W1[i];
    }
    for (i=0;i<m+2;i++) {
      BP[m+1][i]=V2[i];
      BP[i][m+1]=W2[i];
    }
    staFunc->CopyMatrix(BP,BT);
    staFunc->Decomp(BT);
    staFunc->ClearVector(VB2);
    VB2[m+1]=1.0;
    staFunc->Solve(BT,VB2);   
    test[3]=VB2[m+1];
  }   
  return 0;
}