#define CS	1	/* check subscripts on */

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

#include "_linalg.h"	/* for access std vectors */

/* Vector and Matrix are type from _linalg.h */

#include <values.h>

#define Var		FloatHi
#define VarPtr		FloatHiPtr
typedef VarPtr VectPtr;

#ifdef WIN_MSWIN

static FILE *deb=NULL;

static void DebOpen(void) {
  if (!deb) deb=fopen("deb","wt");
}
#endif

#if CS
Local(Boolean) Check(Int2 sub, Int2 lim, Int2 line) {
  if (sub<0 || sub>=lim) {
    if (Message(MSG_OKC,"Index is out of range:\n"
                   " %i not in [0,%i)\n"
		   "_linalg.c, line %i",
		   (int)sub,(int)lim,(int)line) == ANS_CANCEL) {
      abort();
    } else return TRUE;
  }
  return FALSE;
}
#endif


/***************************************/

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

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

/* phase dim and number of mesh points */
Local(Int2) nphase,nmesh,nappend;

/* Decomp's work stuff */
Local(Matrix) T=NULL,Z1=NULL,Z2=NULL,A0=NULL,ZZ2=NULL,tt=NULL,ttfi=NULL;
Local(Vector) t=NULL,ttt=NULL,tttt=NULL,beta=NULL,beta1=NULL,wdet=NULL;

/* Package initialization function */
Entry(void) rd_Init(StdLinAlgDataPtr Stdalg, Int2 Nphase, Int2 Nmesh, Int2 Nappend) {
  /* save global data */
  stdalg=Stdalg;
  nphase=Nphase;
  nmesh=Nmesh;
  nappend=Nappend;
  /* Create Decomp's work matrices and vectors */
  if (!T) {
    T=stdalg->CreateMatrix(nphase,nphase);
    Z1=stdalg->CreateMatrix(nphase,nphase);
    Z2=stdalg->CreateMatrix(nphase,nphase);
    A0=stdalg->CreateMatrix(nphase,nphase);
    ZZ2=stdalg->CreateMatrix(nphase,nappend+1);
    tt=stdalg->CreateMatrix(nphase,nappend+1);
    ttfi=stdalg->CreateMatrix(nappend+1,nappend+1);
    t=stdalg->CreateVector(nphase);
    ttt=stdalg->CreateVector(nphase);
    tttt=stdalg->CreateVector(nphase);
    beta=stdalg->CreateVector(nappend+1);
    beta1=stdalg->CreateVector(nappend+1);
    wdet=stdalg->CreateVector(nphase*nmesh+nappend+1);
  }
}

/* Package clean up function */
Entry(void) rd_Term(void) {
  /* Decomp */
  T=stdalg->DestroyMatrix(T);
  Z1=stdalg->DestroyMatrix(Z1);
  Z2=stdalg->DestroyMatrix(Z2);
  A0=stdalg->DestroyMatrix(A0);
  ZZ2=stdalg->DestroyMatrix(ZZ2);
  tt=stdalg->DestroyMatrix(tt);
  ttfi=stdalg->DestroyMatrix(ttfi);
  t=stdalg->DestroyVector(t);
  ttt=stdalg->DestroyVector(ttt);
  tttt=stdalg->DestroyVector(tttt);
  beta=stdalg->DestroyVector(beta);
  beta1=stdalg->DestroyVector(beta1);
  wdet=stdalg->DestroyVector(wdet);
}
/* A matrix */
/*
   Logical structure of a matrix (n=nmesh):

   a[0] b[0] c[0]   0    0  ...    0      0      0      0      0   delta[0]
   a[1] b[1] c[1]   0    0  ...    0      0      0      0      0   delta[1]
     0  a[2] b[2] c[2]   0  ...    0      0      0      0      0   delta[2]
     0    0  a[3] b[3] c[3] ...    0      0      0      0      0   delta[3]
     .    .    .    .    .  ...    .      .      .      .      .
     0    0    0    0    0  ... a[n-4] b[n-4] c[n-4]    0      0   delta[n-4]
     0    0    0    0    0  ...    0   a[n-3] b[n-3] c[n-3]    0   delta[n-3]
     0    0    0    0    0  ...    0      0   a[n-2] b[n-2] c[n-2] delta[n-2]
     0    0    0    0    0  ...    0      0   a[n-1] b[n-1] c[n-1] delta[n-1]
   p[0] p[1] p[2] p[3] p[4] ... p[n-5] p[n-4] p[n-3] p[n-2] p[n-1] fi
   
   dimensions:
   a[i], b[i], c[i]   nphase    x nphase
   delta[i]           nphase    x 1+nappend
   p[i]               1+nappend x nphase
   fi                 1+nappend x 1+nappend
*/
typedef struct {
  Int2 n,m;			/* dimensions */
  Matrix R;			/* from Decomp to Solve */
  Matrix PNTR a;
  Matrix PNTR b;
  Matrix PNTR c;
  Matrix PNTR delta;
  Matrix PNTR p;
  Matrix PNTR fi;
} rdMatrix, PNTR rdMatrixPtr;

typedef rdMatrixPtr MatrPtr;

#define MATRSIZE(n,m)	(sizeof(rdMatrix)+		\
			(5*nmesh+1)*sizeof(Matrix)	\
			)


Entry(void) rd_PrintMatrix(MatrPtr mp, CharPtr title, Int2 Level) {
  int i;
  Char buf[20];
  if (mp==NULL) {
#ifdef WIN_MSWIN
    DebOpen();
    fprintf(deb,"\n%s=NULL\n",title);
#else
    printf("\n%s=NULL\n",title);
#endif
    return;
  }
#ifdef WIN_MSWIN
  DebOpen();
  fprintf(deb,"\n--------------------------------------\n%s\n",title);
#else
  printf("\n--------------------------------------\n%s\n",title);
#endif
  for (i=0; i<nmesh; i++) {
    sprintf(buf,"a[%i]=",i);
    stdalg->PrintMatrix(mp->a[i],buf,0);
    sprintf(buf,"b[%i]=",i);
    stdalg->PrintMatrix(mp->b[i],buf,0);
    sprintf(buf,"c[%i]=",i);
    stdalg->PrintMatrix(mp->c[i],buf,0);
    sprintf(buf,"delta[%i]=",i);
    stdalg->PrintMatrix(mp->delta[i],buf,0);
  }
  for (i=0; i<nmesh; i++) {
    sprintf(buf,"p[%i]=",i);
    stdalg->PrintMatrix(mp->p[i],buf,0);
  }
  stdalg->PrintMatrix(mp->fi[0],"fi=",0);
  stdalg->PrintMatrix(mp->R,"R=",0);
}

/* Creates a Matrix of dimensions n,m */
/* IMPORTANT NOTE: we always create a square matrix */
Entry(MatrPtr) rd_CreateMatrix(int n, int m DEBPAR) {
  rdMatrixPtr mp;
  Matrix PNTR ma;
  Int2 i;
  mp = MemNewDeb(MATRSIZE(d,d) DEBARG);
  mp->n=(Int2)n;
  mp->m=(Int2)m;
  mp->R=stdalg->CreateMatrix(nphase,nphase);
  ma=(Matrix PNTR)((CharPtr)mp+sizeof(rdMatrix));
  mp->a=ma;
  for (i=0; i<nmesh; i++) *ma++=stdalg->CreateMatrix(nphase,nphase);
  mp->b=ma;
  for (i=0; i<nmesh; i++) *ma++=stdalg->CreateMatrix(nphase,nphase);
  mp->c=ma;
  for (i=0; i<nmesh; i++) *ma++=stdalg->CreateMatrix(nphase,nphase);
  mp->delta=ma;
  for (i=0; i<nmesh; i++) *ma++=stdalg->CreateMatrix(nphase,nappend+1);
  mp->p=ma;
  for (i=0; i<nmesh; i++) *ma++=stdalg->CreateMatrix(nappend+1,nphase);
  mp->fi=ma;
  *ma++=stdalg->CreateMatrix(nappend+1,nappend+1);
  return mp;
}

/* Destroys a Matrix and free the memory it occupies */
Entry(MatrPtr) rd_DestroyMatrix(MatrPtr mp) {
  Matrix PNTR ma;
  Int2 i;
  if (mp) {
    stdalg->DestroyMatrix(mp->R);				/* R */
    ma=(Matrix PNTR)((CharPtr)mp+sizeof(rdMatrix));
    for (i=0; i<nmesh; i++) stdalg->DestroyMatrix(*ma++);	/* a */
    for (i=0; i<nmesh; i++) stdalg->DestroyMatrix(*ma++);	/* b */
    for (i=0; i<nmesh; i++) stdalg->DestroyMatrix(*ma++);	/* c */
    for (i=0; i<nmesh; i++) stdalg->DestroyMatrix(*ma++);	/* delta */
    for (i=0; i<nmesh; i++) stdalg->DestroyMatrix(*ma++);	/* p */
    stdalg->DestroyMatrix(*ma++);				/* fi */
    MemFree(mp);
  }
  return NULL;
}

/* Copies a Matrix */
Entry(MatrPtr) rd_CopyMatrix(MatrPtr from, MatrPtr to DEBPAR) {
  Matrix PNTR fma;
  Matrix PNTR tma;
  Int2 i;
  if (!to) to=rd_CreateMatrix(from->n,from->m DEBARG);
  stdalg->CopyMatrix(from->R,to->R);
  fma=(Matrix PNTR)((CharPtr)from+sizeof(rdMatrix));
  tma=(Matrix PNTR)((CharPtr)to+sizeof(rdMatrix));
  for (i=0; i<nmesh; i++) stdalg->CopyMatrix(*fma++,*tma++);	/* a */
  for (i=0; i<nmesh; i++) stdalg->CopyMatrix(*fma++,*tma++);	/* b */
  for (i=0; i<nmesh; i++) stdalg->CopyMatrix(*fma++,*tma++);	/* c */
  for (i=0; i<nmesh; i++) stdalg->CopyMatrix(*fma++,*tma++);	/* delta */
  for (i=0; i<nmesh; i++) stdalg->CopyMatrix(*fma++,*tma++);	/* p */
  stdalg->CopyMatrix(*fma++,*tma++);				/* fi */
  return to;
}

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

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

/* Multilpy Matrix by vector */
Entry(VectPtr) rd_MultiplyMatrixVector(MatrPtr mp, VectPtr vp, VectPtr vpres DEBPAR) {
  Var s;
  VarPtr vres;
  Int2 i,j,k,l,m,n;
  n=stdalg->GetVectorDim(vp);
  if (!vpres) vpres=stdalg->CreateVector(n);
  vres=vpres;
  /* a[0],b[0],c[0] 0 0 0 ... d[0] */
  for (i=0; i<nphase; i++) {
    s=0;
    for (j=0; j<nphase; j++)
      s+=mp->a[0][i][j]*vp[j]+
	 mp->b[0][i][j]*vp[nphase+j]+
	 mp->c[0][i][j]*vp[2*nphase+j];
    for (j=0; j<nappend+1; j++)
      s+=mp->delta[0][i][j]*vp[n-nappend-1+j];
    *vres++=s;
  }
  /* a[i], b[i], c[i] 0 0 0 ... d[i]   (i=1,2,...nmesh-2) */
  for (k=1,l=0; k<nmesh-1; k++,l+=nphase)
    for (i=0; i<nphase; i++) {
      s=0;
      for (j=0; j<nphase; j++)
	s+=mp->a[k][i][j]*vp[l+j]+
	   mp->b[k][i][j]*vp[l+nphase+j]+
	   mp->c[k][i][j]*vp[l+2*nphase+j];
      for (j=0; j<nappend+1; j++)
	s+=mp->delta[k][i][j]*vp[n-nappend-1+j];
      *vres++=s;
    }
  /* a[nmesh-1], b[nmesh-1], c[nmesh-1] 0 0 0 ... d[nmesh-1] */
  for (i=0,l=(nmesh-3)*nphase; i<nphase; i++) {
    s=0;
    for (j=0; j<nphase; j++)
      s+=mp->a[nmesh-1][i][j]*vp[l+j]+
	 mp->b[nmesh-1][i][j]*vp[l+nphase+j]+
	 mp->c[nmesh-1][i][j]*vp[l+2*nphase+j];
    for (j=0; j<nappend+1; j++)
      s+=mp->delta[nmesh-1][i][j]*vp[n-nappend-1+j];
    *vres++=s;
  }
  m=mp->n-nmesh*nphase;
  /* ASSERT: m==nappend+1 for square matrix, m==nappend otherwise */
  for (i=0; i<m; i++) {
    s=0;
    for (k=0,l=0; k<nmesh; k++,l+=nphase)
      for (j=0; j<nphase; j++)
	s+=mp->p[k][i][j]*vp[l+j];
    for (j=0; j<nappend+1; j++)
      s+=mp->fi[0][i][j]*vp[n-nappend-1+j];
    *vres++=s;
  }
  return vpres;
}

/* Appends a vector to a Matrix as the last row */
Entry(MatrPtr) rd_AppendMatrixVector(MatrPtr mp, VectPtr vp, MatrPtr mpres DEBPAR) {
  Int2 k,l,m=mp->m;
  if (!mpres) mpres=rd_CreateMatrix(m,m DEBARG);
  if (mpres->n<mpres->m) (mpres->n)++;	/* the check needed when mpres!=NULL on entry */
  rd_CopyMatrix(mp,mpres DEBARG);
  /* overwrite the last row */
  for (k=0,l=0; k<nmesh; k++,l+=nphase)
    MemMove(mpres->p[k][nappend],vp+l,sizeof(Var)*nphase);
  MemMove(mpres->fi[0][nappend],vp+l,sizeof(Var)*(nappend+1));
  return mpres;
}

Entry(Matrix) rd_GetMatrixA(MatrPtr mp, int n) {
  return mp->a[n];
}

Entry(Matrix) rd_GetMatrixB(MatrPtr mp, int n) {
  return mp->b[n];
}

Entry(Matrix) rd_GetMatrixC(MatrPtr mp, int n) {
  return mp->c[n];
}

Entry(Matrix) rd_GetMatrixDelta(MatrPtr mp, int n) {
  return mp->delta[n];
}

Entry(Matrix) rd_GetMatrixP(MatrPtr mp, int n) {
  return mp->p[n];
}

Entry(Matrix) rd_GetMatrixFi(MatrPtr mp) {
  return mp->fi[0];
}

/***************/
/* Determinant */

/* Copy all elements from im-th row of S->p[0] to wdet */
Local(void) CopyToWdet(MatrPtr S, Int2 im) {
  Int2 i,j,k;
  for (i=k=0; i<nmesh; i++)
    for (j=0; j<nphase; j++)
      wdet[k++]=S->p[i][im][j];
  for (j=0; j<nappend+1; j++)
    wdet[k++]=S->fi[0][im][j];
}


Local(Int2) Determinant (MatrPtr S,  VarPtr deteq,  VarPtr det) {
Int2 im,jm,i,j,k,l,isimb,imesh,nm,nm3,nm2,nm1,np,m,sign=1;
Var s,w;
Boolean zeroeq=FALSE;


np=nphase;
nm=nmesh;
  nm3=nmesh-3;
  nm2=nmesh-2;
  nm1=nmesh-1;
m=1+nappend;
*det=0;


#define d delta

  /* Pervie dve kletotchnih stroki. Pervii kletotchnii stolbets. */
  for(k=0;k<np;k++) {
    im=0;
    jm=0;
    /* Poisk glavnogo elementa */
    s=0.0;
    for(j=k; j<np; j++){      /* po matritse a[0] */
      for(i=k; i<np; i++) 
        if(fabs(S->a[0][i][j])>s) {s=fabs(S->a[0][i][j]);
	                           im=i;
	                           jm=j;
				   isimb=0;}
      /* po matritse a[1] */
      for(i=0; i<np; i++) 
        if(fabs(S->a[1][i][j])>s) {s=fabs(S->a[1][i][j]);
			           im=i;
				   jm=j;
				   isimb=1;}
    }
    if (s==0.0) {
      zeroeq=TRUE;
      for(j=k; j<np; j++) {      /* po matritse p[0] */
	for(i=0; i<m; i++) 
	  if(fabs(S->p[0][i][j])>s) {s=fabs(S->p[0][i][j]);
	                           im=i;
	                           jm=j;
				   isimb=-1;}
      }
      if (s==0.0) {
	*det=*deteq=0;
	return 0;
      }
    }
    /* Perestanovka stolbtsov */
    if(jm != k && k<np-1) {
      sign=-sign;
      for(i=k; i<np; i++){/* Perestanovka stolbtsov a[0] */
         w=S->a[0][i][jm];
         S->a[0][i][jm]=S->a[0][i][k];
         S->a[0][i][k]=w;}
      for(i=0; i<np; i++){/* Perestanovka stolbtsov a[1] */
         w=S->a[1][i][jm];
         S->a[1][i][jm]=S->a[1][i][k];
         S->a[1][i][k]=w;}
      for(i=0; i<m; i++){/* Perestanovka stolbtsov p[0] */
         w=S->p[0][i][j];
         S->p[0][i][jm]=S->p[0][i][k];
         S->p[0][i][k]=w;}
    }
    switch (isimb) {
    /* Perestanovka strok */	    
      case 0:		/* glavnaya stroka - v a[0] */		    
	if(im != k && k<np-1) {
	  sign=-sign;
	  for(j=k; j<np; j++){/* Perestanovka strok a[0] */
		      w=S->a[0][im][j];
		      S->a[0][im][j]=S->a[0][k][j];
 		      S->a[0][k][j]=w;}
	  for(j=0; j<np; j++){/* Perestanovka strok b[0] */
		      w=S->b[0][im][j];
		      S->b[0][im][j]=S->b[0][k][j];
		      S->b[0][k][j]=w;}
	  for(j=0; j<np; j++){/* Perestanovka strok c[0] */
		      w=S->c[0][im][j];
		      S->c[0][im][j]=S->c[0][k][j];
	              S->c[0][k][j]=w;}
	  for(j=0; j<m; j++){/* Perestanovka strok d[0] */
		      w=S->d[0][im][j];
		      S->d[0][im][j]=S->d[0][k][j];
		      S->d[0][k][j]=w;}
	}
	break;
      case 1:		/* glavnaya stroka - v a[1] */				  
	sign=-sign;
	for(j=k; j<np; j++){/* Obmen strok a[1] i a[0]*/
	     w=S->a[1][im][j];
             S->a[1][im][j]=S->a[0][k][j];
	     S->a[0][k][j]=w;}
	for(j=0; j<np; j++){/* Obmen strok b[1] i b[0] */
	     w=S->b[1][im][j];
	     S->b[1][im][j]=S->b[0][k][j];
	     S->b[0][k][j]=w;}
	for(j=0; j<np; j++){/* Obmen strok c[1] i c[0]*/
	     w=S->c[1][im][j];
             S->c[1][im][j]=S->c[0][k][j];
	     S->c[0][k][j]=w;}
	for(j=0; j<m; j++) {/* Obmen strok d[1] i d[0] */
	     w=S->d[1][im][j];
             S->d[1][im][j]=S->d[0][k][j];
	     S->d[0][k][j]=w;}
	break;
      case -1:		/* glavnaya stroka v p[0] */
	sign=-sign;
	CopyToWdet(S,im);	/* copy it to wdet */
	for(j=k; j<np; j++){/* Obmen strok p[0] i a[0]*/
	     w=S->p[0][im][j];
             S->p[0][im][j]=S->a[0][k][j];
	     S->a[0][k][j]=w;}
	for(j=0; j<np; j++){/* Obmen strok p[1] i b[0] */
	     w=S->p[1][im][j];
	     S->p[1][im][j]=S->b[0][k][j];
	     S->b[0][k][j]=w;}
	for(j=0; j<np; j++){/* Obmen strok p[2] i c[0]*/
	     w=S->p[2][im][j];
             S->p[2][im][j]=S->c[0][k][j];
	     S->c[0][k][j]=w;}
	for(j=0; j<m; j++) {/* Obmen strok fi[0] i d[0] */
	     w=S->fi[0][im][j];
             S->fi[0][im][j]=S->d[0][k][j];
	     S->d[0][k][j]=w;}
	for (i=3; i<nm; i++)	/* Obnulenie stroki im v p[3..nmesh] */
	  for(j=0; j<np; j++) S->p[i][im][j]=0;
	break;
    }		/* switch */
    if (k<np-1) {
      /* Otchistka stolbtsa a[0] */
      for(i=k+1;i<np;i++) {
        s=S->a[0][i][k]/S->a[0][k][k];
        for(j=k;j<np;j++) S->a[0][i][j]-=s*S->a[0][k][j];
        for(j=0;j<np;j++) S->b[0][i][j]-=s*S->b[0][k][j];
        for(j=0;j<np;j++) S->c[0][i][j]-=s*S->c[0][k][j];
        for(j=0;j<m;j++)  S->d[0][i][j]-=s*S->d[0][k][j];
      }
    }
    /* Otchistka stolbtsa a[1] */
    for(i=0;i<np;i++) {
      s=S->a[1][i][k]/S->a[0][k][k];
      for(j=k;j<np;j++) S->a[1][i][j]-=s*S->a[0][k][j];
      for(j=0;j<np;j++) S->b[1][i][j]-=s*S->b[0][k][j];
      for(j=0;j<np;j++) S->c[1][i][j]-=s*S->c[0][k][j];
      for(j=0;j<m;j++)  S->d[1][i][j]-=s*S->d[0][k][j];
    }
    /* Otchistka stolbtsa p[0] */
    for(i=0;i<m;i++) {
      s=S->p[0][i][k]/S->a[0][k][k];
      for(j=k;j<np;j++) S->p[0][i][j]-=s*S->a[0][k][j];
      for(j=0;j<np;j++) S->p[1][i][j]-=s*S->b[0][k][j]; 
      for(j=0;j<np;j++) S->p[2][i][j]-=s*S->c[0][k][j]; 
      for(j=0;j<m;j++)  S->fi[0][i][j]-=s*S->d[0][k][j];
      if (isimb==-1)
	for (l=3; l<nm; l++)
	  for(j=0;j<np;j++) S->p[l][i][j]-=s*wdet[l*np+j]; 
    }
    /* Determinant */
    if (S->a[0][k][k]<0) sign=-sign;
    *det+=log(fabs(S->a[0][k][k]));
  }	/* for (k=... */

  
  /* Kletotchnie stroki 1 <= imesh <= nmesh-4. Pervii kletotchnii stolbets. */
  for(imesh=1; imesh<nm3; imesh++) {
    for(k=0;k<np;k++){
      im=0;
      jm=0;
      /* Poisk glavnogo elementa */
      s=0.0;
      for(j=k; j<np; j++) {
        /* po matritse b[imesh] */
        for(i=k; i<np; i++) 
          if(fabs(S->b[imesh][i][j])>s) {s=fabs(S->b[imesh][i][j]);
				       im=i;
				       jm=j;
				       isimb=imesh;}
        /* po matritse a[imesh+1] */
	for(i=0; i<np; i++) 
          if(fabs(S->a[imesh+1][i][j])>s) {s=fabs(S->a[imesh+1][i][j]);
					 im=i;
				         jm=j;
				         isimb=imesh+1;}
      }	/* for (j=k; j<np; j++) */
      if (s==0.0) {
	zeroeq=TRUE;
	for(j=k; j<np; j++) {      /* po matritse p[imesh] */
	  for(i=0; i<m; i++) 
	    if(fabs(S->p[imesh][i][j])>s) {s=fabs(S->p[imesh][i][j]);
	                           im=i;
	                           jm=j;
				   isimb=-1;}
	}
	if (s==0.0) {
	  *det=*deteq=0;
	  return 0;
	}
      }
      /* Perestanovka stolbtsov */
      if(jm != k && k<np-1) {
	sign=-sign;
        for(i=k; i<np; i++){/* Perestanovka stolbtsov b[imesh] */
	  w=S->b[imesh][i][jm];
          S->b[imesh][i][jm]=S->b[imesh][i][k];
          S->b[imesh][i][k]=w;}
        for(i=0; i<np; i++){/* Perestanovka stolbtsov a[imesh+1] */
          w=S->a[imesh+1][i][jm];
          S->a[imesh+1][i][jm]=S->a[imesh+1][i][k];
          S->a[imesh+1][i][k]=w;}
        for(i=0; i<m; i++){/* Perestanovka stolbtsov p[imesh] */
          w=S->p[imesh][i][jm];
          S->p[imesh][i][jm]=S->p[imesh][i][k];
          S->p[imesh][i][k]=w;}
      }
      /* Perestanovka strok */	    
      if(isimb == imesh) {
        /* glavnaya stroka - v b[imesh] */		    
        if(im != k && k<np-1) {
	  sign=-sign;
	  for(j=k; j<np; j++){/* Perestanovka strok b[imesh] */
		      w=S->b[imesh][im][j];
		      S->b[imesh][im][j]=S->b[imesh][k][j];
		      S->b[imesh][k][j]=w;}
	  for(j=0; j<np; j++){/* Perestanovka strok c[imesh] */
		      w=S->c[imesh][im][j];
		      S->c[imesh][im][j]=S->c[imesh][k][j];
		      S->c[imesh][k][j]=w;}
	  for(j=0; j<m; j++){/* Perestanovka strok d[imesh] */
		      w=S->d[imesh][im][j];
		      S->d[imesh][im][j]=S->d[imesh][k][j];
		      S->d[imesh][k][j]=w;}
	}
      }
      if (isimb == imesh+1) {
        /* glavnaya stroka - v a[imesh+1] */
	sign=-sign;
	for(j=k; j<np; j++){/* Obmen strok a[imesh+1] i b[imesh] */
	      w=S->a[imesh+1][im][j];
              S->a[imesh+1][im][j]=S->b[imesh][k][j];
	      S->b[imesh][k][j]=w;}
	for(j=0; j<np; j++){/* Obmen strok b[imesh+1] i c[imesh] */
	      w=S->b[imesh+1][im][j];
	      S->b[imesh+1][im][j]=S->c[imesh][k][j];
	      S->c[imesh][k][j]=w;}
	for(j=0; j<np; j++){/* Zapominanie stroki c[imesh+1] v vectore t i ee obnulenie*/
	      t[j]=S->c[imesh+1][im][j];
              S->c[imesh+1][im][j]=0.0;}
	for(j=0; j<m; j++) {/* Obmen strok d[imesh+1] i d[imesh] */
	      w=S->d[imesh+1][im][j];
	      S->d[imesh+1][im][j]=S->d[imesh][k][j];
	      S->d[imesh][k][j]=w;}
      }
      if (isimb == -1) {
        /* glavnaya stroka - v p[imesh] */
	sign=-sign;
	CopyToWdet(S,im);	/* copy it to wdet */
	for(j=k; j<np; j++){/* Obmen strok p[imesh] i b[imesh]*/
	     w=S->p[imesh][im][j];
             S->p[imesh][im][j]=S->b[imesh][k][j];
	     S->b[imesh][k][j]=w;}
	for(j=0; j<np; j++){/* Obmen strok p[imesh+1] i c[imesh] */
	     w=S->p[imesh+1][im][j];
	     S->p[imesh+1][im][j]=S->c[imesh][k][j];
	     S->c[imesh][k][j]=w;}
	for(j=0; j<m; j++) {/* Obmen strok fi[0] i d[imesh] */
	     w=S->fi[0][im][j];
             S->fi[0][im][j]=S->d[imesh][k][j];
	     S->d[imesh][k][j]=w;}
	for (i=imesh+2; i<nm; i++)	/* Obnulenie stroki im v p[imesh+2..nmesh] */
	  for(j=0; j<np; j++) S->p[i][im][j]=0;
      }
      if (k<np-1) {
        /* Otchistka stolbtsa b[imesh] */
	for(i=k+1;i<np;i++) {
	  s=S->b[imesh][i][k]/S->b[imesh][k][k];
	  for(j=k;j<np;j++) S->b[imesh][i][j]-=s*S->b[imesh][k][j];
	  for(j=0;j<np;j++) S->c[imesh][i][j]-=s*S->c[imesh][k][j];
	  for(j=0;j<m;j++)  S->d[imesh][i][j]-=s*S->d[imesh][k][j];
        }
      }		/* if (k<np-1) */
      /* Otchistka stolbtsa a[imesh+1] */
      for(i=0;i<np;i++) {
	s=S->a[imesh+1][i][k]/S->b[imesh][k][k];
	for(j=k;j<np;j++) S->a[imesh+1][i][j]-=s*S->b[imesh][k][j];
	for(j=0;j<np;j++) S->b[imesh+1][i][j]-=s*S->c[imesh][k][j];
	if(isimb == imesh+1) 
	  for(j=0;j<np;j++) S->c[imesh+1][i][j]-=s*t[j];
	if (isimb==-1)
	  for(j=0;j<np;j++) S->c[imesh+1][i][j]-=s*wdet[(imesh+2)*np+j]; 
	for(j=0;j<m;j++)  S->d[imesh+1][i][j]-=s*S->d[imesh][k][j];
      }
      /* Otchistka stolbtsa p[imesh] */
      for(i=0;i<m;i++) {
	s=S->p[imesh][i][k]/S->b[imesh][k][k];
	for(j=k;j<np;j++) S->p[imesh][i][j]-=s*S->b[imesh][k][j];
	for(j=0;j<np;j++) S->p[imesh+1][i][j]-=s*S->c[imesh][k][j];
	if(isimb == imesh+1) 
	  for(j=0;j<np;j++) S->p[imesh+2][i][j]-=s*t[j]; 
	for(j=0;j<m;j++)  S->fi[0][i][j]-=s*S->d[imesh][k][j];
	if (isimb==-1)
	  for (l=imesh+2; l<nm; l++)
	    for(j=0;j<np;j++) S->p[l][i][j]-=s*wdet[l*np+j]; 
      }
      /* Determinant */
      if (S->b[imesh][k][k]<0) sign=-sign;
      *det+=log(fabs(S->b[imesh][k][k]));
    }		/* for(k=0;k<np;k++) */
  }		/*  for(imesh=1; imesh<nm3; imesh++) */

  /* Kletotchnaya stroka nm3. Pervii kletotchnii stolbets. */
  for(k=0;k<np;k++) {
    im=0;
    jm=0;
    /* Poisk glavnogo elementa */
    s=0.0;
    for(j=k; j<np; j++){
      /* po matritse b[nm3] */
      for(i=k; i<np; i++) 
        if(fabs(S->b[nm3][i][j])>s) {s=fabs(S->b[nm3][i][j]);
			             im=i;
				     jm=j;
				     isimb=nm3;}
      /* po matritse a[nm2] */
      for(i=0; i<np; i++) 
        if(fabs(S->a[nm2][i][j])>s) {s=fabs(S->a[nm2][i][j]);
			             im=i;
				     jm=j;
			             isimb=nm2;}
      /* po matritse a[nm1] */
      for(i=0; i<np; i++) 
        if(fabs(S->a[nm1][i][j])>s) {s=fabs(S->a[nm1][i][j]);
				     im=i;
				     jm=j;
	                             isimb=nm1;}
    }		/* for(j=k; j<np; j++) */
    if (s==0.0) {
      zeroeq=TRUE;
      for(j=k; j<np; j++) {      /* po matritse p[nm3] */
	 for(i=0; i<m; i++) 
	  if(fabs(S->p[nm3][i][j])>s) {s=fabs(S->p[nm3][i][j]);
	                           im=i;
	                           jm=j;
				   isimb=-1;}
      }
      if (s==0.0) {
	*det=*deteq=0;
	return 0;
      }
    }
    /* Perestanovka stolbtsov */
    if (jm != k && k<np-1) {
      sign=-sign;
      for(i=k; i<np; i++){/* Perestanovka stolbtsov b[nm3]  */
        w=S->b[nm3][i][jm];
        S->b[nm3][i][jm]=S->b[nm3][i][k];
        S->b[nm3][i][k]=w;}
     for(i=0; i<np; i++){/* Perestanovka stolbtsov a[nm2] */
        w=S->a[nm2][i][jm];
        S->a[nm2][i][jm]=S->a[nm2][i][k];
        S->a[nm2][i][k]=w;}
     for(i=0; i<np; i++){/* Perestanovka stolbtsov a[nm1] */
       w=S->a[nm1][i][jm];
       S->a[nm1][i][jm]=S->a[nm1][i][k];
       S->a[nm1][i][k]=w;}
     for(i=0; i<m; i++){/* Perestanovka stolbtsov p[nm3] */
       w=S->p[nm3][i][jm];
       S->p[nm3][i][jm]=S->p[nm3][i][k];
       S->p[nm3][i][k]=w;}
    }		/* if (jm != k && k<np-1) */	  
    /* Perestanovka strok */	    
    if(isimb == nm3) {
      /* glavnaya stroka - v b[nm3] */
      if(im != k && k<np-1) {
	sign=-sign;
	for(j=k; j<np; j++) {/* Perestanovka strok b[nm3] */
		      w=S->b[nm3][im][j];
		      S->b[nm3][im][j]=S->b[nm3][k][j];
		      S->b[nm3][k][j]=w;}
	for(j=0; j<np; j++) {/* Perestanovka strok c[nm3] */
		      w=S->c[nm3][im][j];
		      S->c[nm3][im][j]=S->c[nm3][k][j];
		      S->c[nm3][k][j]=w;}
	for(j=0; j<m; j++) {/* Perestanovka strok d[nm3] */
		      w=S->d[nm3][im][j];
		      S->d[nm3][im][j]=S->d[nm3][k][j];
		      S->d[nm3][k][j]=w;}
      }		/* if(im != k && k<np-1) */
    }		/* if(isimb == nm3) */
    if(isimb == nm2) {
      /* glavnaya stroka - v a[nm2] */
      sign=-sign;
      for(j=k; j<np; j++){/* Obmen strok a[nm2] i b[nm3] */
	    w=S->a[nm2][im][j];
            S->a[nm2][im][j]=S->b[nm3][k][j];
	    S->b[nm3][k][j]=w;}
      for(j=0; j<np; j++){/* Obmen strok b[nm2] i c[nm3] */
	    w=S->b[nm2][im][j];
            S->b[nm2][im][j]=S->c[nm3][k][j];
	    S->c[nm3][k][j]=w;}
      for(j=0; j<np; j++){/* Zapominanie stroki c[nm2] v vectore t i ee obnulenie */
	    t[j]=S->c[nm2][im][j];
            S->c[nm2][im][j]=0.0;}
      for(j=0; j<m; j++){/* Obmen strok d[nm2] i d[nm3] */
	    w=S->d[nm2][im][j];
            S->d[nm2][im][j]=S->d[nm3][k][j];
	    S->d[nm3][k][j]=w;}
    }		/* if(isimb == nm2) */
    if(isimb == nm1) {
      /* glavnaya stroka - v a[nm1] */				  
      sign=-sign;
      for(j=k; j<np; j++){/* obmen strok a[nm1] i b[nm3] */
	     w=S->a[nm1][im][j];
             S->a[nm1][im][j]=S->b[nm3][k][j];
	     S->b[nm3][k][j]=w;}
      for(j=0; j<np; j++){/* Obmen strok b[nm1] i c[nm3] */
	     w=S->b[nm1][im][j];
             S->b[nm1][im][j]=S->c[nm3][k][j];
	     S->c[nm3][k][j]=w;}
      for(j=0; j<np; j++){/* Zapominanie stroki c[nm1] v vectore t i ee obnulenie */
	     t[j]=S->c[nm1][im][j];
             S->c[nm1][im][j]=0.0;}
      for(j=0; j<m; j++){/* Obmen strok d[nm1] i d[nm3] */
	     w=S->d[nm1][im][j];
             S->d[nm1][im][j]=S->d[nm3][k][j];
	     S->d[nm3][k][j]=w;}
    }		/* if(isimb == nm1) */
    if (isimb==-1) {
      /* glavnaya stroka - v p[nm3] */				  
      sign=-sign;
      for(j=k; j<np; j++){/* obmen strok p[nm3] i b[nm3] */
	     w=S->p[nm3][im][j];
             S->p[nm3][im][j]=S->b[nm3][k][j];
	     S->b[nm3][k][j]=w;}
      for(j=0; j<np; j++){/* Obmen strok p[nm2] i c[nm3] */
	     w=S->p[nm2][im][j];
             S->p[nm2][im][j]=S->c[nm3][k][j];
	     S->c[nm3][k][j]=w;}
      for(j=0; j<np; j++){/* Zapominanie stroki p[nm1] v vectore t i ee obnulenie */
	     t[j]=S->p[nm1][im][j];
             S->p[nm1][im][j]=0.0;}
      for(j=0; j<m; j++){/* Obmen strok fi[0] i d[nm3] */
	     w=S->fi[0][im][j];
             S->fi[0][im][j]=S->d[nm3][k][j];
	     S->d[nm3][k][j]=w;}
    }
    if (k<np-1) {
      /* Otchistka stolbtsa b[nm3] */
      for(i=k+1;i<np;i++) {
	s=S->b[nm3][i][k]/S->b[nm3][k][k];
	for(j=k;j<np;j++) S->b[nm3][i][j]-=s*S->b[nm3][k][j];
	for(j=0;j<np;j++) S->c[nm3][i][j]-=s*S->c[nm3][k][j];
	for(j=0;j<m;j++)  S->d[nm3][i][j]-=s*S->d[nm3][k][j];
      }
    }		/* if (k<np-1) */
    /* Otchistka stolbtsa a[nm2] */
    for(i=0;i<np;i++) {
      s=S->a[nm2][i][k]/S->b[nm3][k][k];
      for(j=k;j<np;j++) S->a[nm2][i][j]-=s*S->b[nm3][k][j];
      for(j=0;j<np;j++) S->b[nm2][i][j]-=s*S->c[nm3][k][j];
      if(isimb != nm3) for(j=0;j<np;j++) S->c[nm2][i][j]-=s*t[j];
      for(j=0;j<m;j++)  S->d[nm2][i][j]-=s*S->d[nm3][k][j];
    }
    /* Otchistka stolbtsa a[nm1] */
    for(i=0;i<np;i++) {
      s=S->a[nm1][i][k]/S->b[nm3][k][k];
      for(j=k;j<np;j++) S->a[nm1][i][j]-=s*S->b[nm3][k][j];
      for(j=0;j<np;j++) S->b[nm1][i][j]-=s*S->c[nm3][k][j];
      if(isimb != nm3) for(j=0;j<np;j++) S->c[nm1][i][j]-=s*t[j];
      for(j=0;j<m;j++)  S->d[nm1][i][j]-=s*S->d[nm3][k][j];
    }
    /* Otchistka stolbtsa p[nm3] */
    for(i=0;i<m;i++) {
      s=S->p[nm3][i][k]/S->b[nm3][k][k];
      for(j=k;j<np;j++) S->p[nm3][i][j]-=s*S->b[nm3][k][j];
      for(j=0;j<np;j++) S->p[nm2][i][j]-=s*S->c[nm3][k][j];
      if(isimb != nm3) for(j=0;j<np;j++) S->p[nm1][i][j]-=s*t[j];
      for(j=0;j<m;j++)  S->fi[0][i][j]-=s*S->d[nm3][k][j];
    }
    /* Determinant */
    if (S->b[nm3][k][k]<0) sign=-sign;
    *det+=log(fabs(S->b[nm3][k][k]));
  }		/* for(k=0;k<np;k++) */
   
   
   
  /* Kletotchnaya stroka nm2. Pervii kletotchnii stolbets. */
  for(k=0;k<np;k++){
    im=0;
    jm=0;
    /* Poisk glavnogo elementa */
    s=0.0;
    for(j=k; j<np; j++){
      /* po matritse b[nm2] */
      for(i=k; i<np; i++) 
        if(fabs(S->b[nm2][i][j])>s) {s=fabs(S->b[nm2][i][j]);
		  	             im=i;
				     jm=j;
				     isimb=nm2;}
      /* po matritse b[nm1] */
      for(i=0; i<np; i++) 
	if(fabs(S->b[nm1][i][j])>s) {s=fabs(S->b[nm1][i][j]);
				     im=i;
				     jm=j;
				     isimb=nm1;}		    
    }
    if (s==0.0) {
      zeroeq=TRUE;
      for(j=k; j<np; j++) {      /* po matritse p[nm2] */
	 for(i=0; i<m; i++) 
	  if(fabs(S->p[nm2][i][j])>s) {s=fabs(S->p[nm2][i][j]);
	                           im=i;
	                           jm=j;
				   isimb=-1;}
      }
      if (s==0.0) {
	*det=*deteq=0;
	return 0;
      }
    }
    /* Perestanovka stolbtsov */
    if (jm != k && k<np-1) {
      sign=-sign;
      for(i=k; i<np; i++){/* Perestanovka stolbtsov b[nm2]  */
	w=S->b[nm2][i][jm];
	S->b[nm2][i][jm]=S->b[nm2][i][k];
	S->b[nm2][i][k]=w;}
      for(i=0; i<np; i++){/* Perestanovka stolbtsov b[nm2] */
	w=S->b[nm2][i][jm];
	S->b[nm2][i][jm]=S->b[nm2][i][k];
	S->b[nm2][i][k]=w;}
      for(i=0; i<m; i++){/* Perestanovka stolbtsov v p[nm2] */
	w=S->p[nm2][i][jm];
	S->p[nm2][i][jm]=S->p[nm2][i][k];
	S->p[nm2][i][k]=w;}
    }		/* if (jm != k && k<np-1) */	  
    /* Perestanovka strok */	    
    if(isimb == nm2) {
      /* glavnaya stroka - v b[nm2] */		    
      if(im != k && k<np-1) {
	sign=-sign;
	for(j=k; j<np; j++){/* Perestanovka strok b[nm2] */
		      w=S->b[nm2][im][j];
		      S->b[nm2][im][j]=S->b[nm2][k][j];
		      S->b[nm2][k][j]=w;}
	for(j=0; j<np; j++){/* Perestanovka strok c[nm2] */
		      w=S->c[nm2][im][j];
		      S->c[nm2][im][j]=S->c[nm2][k][j];
		      S->c[nm2][k][j]=w;}
	for(j=0; j<m; j++){/* Perestanovka strok d[nm2] */
		      w=S->d[nm2][im][j];
		      S->d[nm2][im][j]=S->d[nm2][k][j];
		      S->d[nm2][k][j]=w;}
      }		/* if(im != k && k<np-1) */
    }		/* if(isimb == nm3) */
    if (isimb==nm1) {/* glavnaya stroka - v b[nm1] */				  
      sign=-sign;
      for(j=k; j<np; j++){/* Obmen strok b[nm1] i b[nm2] */
	      w=S->b[nm1][im][j];
	      S->b[nm1][im][j]=S->b[nm2][k][j];
	      S->b[nm2][k][j]=w;}
      for(j=0; j<np; j++){/* Obmen strok c[nm1] i c[nm2] */
	      w=S->c[nm1][im][j];
	      S->c[nm1][im][j]=S->c[nm2][k][j];
	      S->c[nm2][k][j]=w;}
      for(j=0; j<m; j++){/* Obmen strok d[nm1] i d[nm2] */
	      w=S->d[nm1][im][j];
	      S->d[nm1][im][j]=S->d[nm2][k][j];
	      S->d[nm2][k][j]=w;}
    }		/* if(isimb == nm2) */
    if (isimb==-1) {
      /* glavnaya stroka - v p[nm2] */
      sign=-sign;
      for(j=k; j<np; j++){/* obmen strok p[nm2] i b[nm2] */
	     w=S->p[nm2][im][j];
             S->p[nm2][im][j]=S->b[nm2][k][j];
	     S->b[nm2][k][j]=w;}
      for(j=0; j<np; j++){/* Obmen strok p[nm1] i c[nm2] */
	     w=S->p[nm1][im][j];
             S->p[nm1][im][j]=S->c[nm2][k][j];
	     S->c[nm2][k][j]=w;}
      for(j=0; j<m; j++){/* Obmen strok fi[0] i d[nm2] */
	     w=S->fi[0][im][j];
             S->fi[0][im][j]=S->d[nm2][k][j];
	     S->d[nm2][k][j]=w;}
    }		/* if(isimb == -1) */
    if (k<np-1) {
      /* Otchistka stolbtsa b[nm2] */
      for(i=k+1;i<np;i++) {
	s=S->b[nm2][i][k]/S->b[nm2][k][k];
	for(j=k;j<np;j++) S->b[nm2][i][j]-=s*S->b[nm2][k][j];
	for(j=0;j<np;j++) S->c[nm2][i][j]-=s*S->c[nm2][k][j];
	for(j=0;j<m;j++)  S->d[nm2][i][j]-=s*S->d[nm2][k][j];
      }
    }		/* if (k<np-1) */
    /* Otchistka stolbtsa b[nm1] */
    for(i=0;i<np;i++) {
      s=S->b[nm1][i][k]/S->b[nm2][k][k];
      for(j=k;j<np;j++) S->b[nm1][i][j]-=s*S->b[nm2][k][j];
      for(j=0;j<np;j++) S->c[nm1][i][j]-=s*S->c[nm2][k][j];
      for(j=0;j<m;j++)  S->d[nm1][i][j]-=s*S->d[nm2][k][j];}              
    /* Otchistka stolbtsa v p[nm2] */
    for(i=0;i<m;i++) {
      s=S->p[nm2][i][k]/S->b[nm2][k][k];
      for(j=k;j<np;j++) S->p[nm2][i][j]-=s*S->b[nm2][k][j];
      for(j=0;j<np;j++) S->p[nm1][i][j]-=s*S->c[nm2][k][j]; 
      for(j=0;j<m;j++)  S->fi[0][i][j]-=s*S->d[nm2][k][j];}
    /* Determinant */
    if (S->b[nm2][k][k]<0) sign=-sign;
    *det+=log(fabs(S->b[nm2][k][k]));
  }		/* for(k=0;k<np;k++) */
   
      
  /* Kletotchnaya stroka nm1. Pervii kletotchnii stolbets. */
  for(k=0;k<np;k++){
    im=0;
    jm=0;
    /* Poisk glavnogo elementa */
    s=0.0;
    for(j=k; j<np; j++) {
      /* po matritse c[nm1] */
      for(i=k; i<np; i++) 
        if(fabs(S->c[nm1][i][j])>s) {s=fabs(S->c[nm1][i][j]);
				     im=i;
				     jm=j;
				     isimb=nm1;}
    }
    if (s==0.0) {
      zeroeq=TRUE;
      for(j=k; j<np; j++) {      /* po matritse p[nm1] */
	for(i=0; i<np; i++) 
	  if(fabs(S->p[nm1][i][j])>s) {s=fabs(S->p[nm1][i][j]);
	                           im=i;
	                           jm=j;
				   isimb=-1;}
      }
      if (s==0.0) {
	*det=*deteq=0;
	return 0;
      }
    }
    /* Perestanovka stolbtsov */
    if(jm != k && k<np-1) {
      sign=-sign;
      for(i=k; i<np; i++){/* Perestanovka stolbtsov c[nm1]  */
	w=S->c[nm1][i][jm];
	S->c[nm1][i][jm]=S->c[nm1][i][k];
	S->c[nm1][i][k]=w;}
      for(i=0; i<m; i++){/* Perestanovka stolbtsov p[nm1] */
	w=S->p[nm1][i][jm];
	S->p[nm1][i][jm]=S->p[nm1][i][k];
	S->p[nm1][i][k]=w;}
    }		/* if(jm != k && k<np-1) */
    /* Perestanovka strok */	    	    
    if (isimb==nm1)
      if(im != k && k<np-1) {
	sign=-sign;
	for(j=k; j<np; j++) {/* Perestanovka strok c[nm1] */
		   w=S->c[nm1][im][j];
		   S->c[nm1][im][j]=S->c[nm1][k][j];
		   S->c[nm1][k][j]=w;}
	for(j=0; j<m; j++){/* Perestanovka strok d[nm1] */
		   w=S->d[nm1][im][j];
		   S->d[nm1][im][j]=S->d[nm1][k][j];
		   S->d[nm1][k][j]=w;}
      }		/* if(im != k && k<np-1) */
    if (isimb==-1) {
      /* glavnaya stroka - v p[nm1] */
      sign=-sign;
      for(j=k; j<np; j++){/* obmen strok p[nm1] i c[nm1] */
	     w=S->p[nm1][im][j];
             S->p[nm1][im][j]=S->c[nm1][k][j];
	     S->c[nm1][k][j]=w;}
      for(j=0; j<m; j++){/* Obmen strok fi[0] i d[nm1] */
	     w=S->fi[0][im][j];
             S->fi[0][im][j]=S->d[nm1][k][j];
	     S->d[nm1][k][j]=w;}
    }
    if (k<np-1) {
      /* Otchistka stolbtsa c[nm1] */
      for(i=k+1;i<np;i++){
	s=S->c[nm1][i][k]/S->c[nm1][k][k];
	for(j=k;j<np;j++) S->c[nm1][i][j]-=s*S->c[nm1][k][j];
	for(j=0;j<m;j++)  S->d[nm1][i][j]-=s*S->d[nm1][k][j];}
    }		/* if (k<np-1) */
    /* Otchistka stolbtsa p[nm1] */
    for(i=0;i<m;i++) {
      s=S->p[nm1][i][k]/S->c[nm1][k][k];
      for(j=k;j<np;j++) S->p[nm1][i][j]-=s*S->c[nm1][k][j];
      for(j=0;j<m;j++) S->fi[0][i][j]-=s*S->d[nm1][k][j];}
    /* Determinant */
    if (S->c[nm1][k][k]<0) sign=-sign;
    *det+=log(fabs(S->c[nm1][k][k]));
  }		/* for(k=0;k<np;k++) */
  if (zeroeq) *deteq=0.0;
  else {
    if (*det<=0) *deteq=sign*exp(*det);
    else *deteq=sign*(1+*det);
  }

  /* Kletka nmesh, t.e "fi" */
  for(k=0;k<m-1;k++) {
    /* Poisk glavnogo elementa */
    s=0.0;
    im=0;
    jm=0;
    for(j=k; j<m; j++) {
      /* po matritse "fi" */
      for(i=k; i<m; i++) 
        if(fabs(S->fi[0][i][j])>s) {s=fabs(S->fi[0][i][j]);
			         im=i;
				 jm=j;}	    
    }
    if (s==0.0) return 6;
    /* Perestanovka stolbtsov */
    if(jm != k && k<np-1){
      sign=-sign;
      for(i=k; i<m; i++){/* Perestanovka stolbtsov "fi"  */
	w=S->fi[0][i][jm];
	S->fi[0][i][jm]=S->fi[0][i][k];
	S->fi[0][i][k]=w;}
    }			  
    /* Perestanovka strok */	    	    
    if(im != k) {
      sign=-sign;
      for(j=k; j<m; j++){/* Perestanovka strok "fi" */
		   w=S->fi[0][im][j];
		   S->fi[0][im][j]=S->fi[0][k][j];
		   S->fi[0][k][j]=w;}
    }
    /* Otchistka stolbtsa "fi" */
    for(i=k+1;i<m;i++) {
      s=S->fi[0][i][k]/S->fi[0][k][k];
      for(j=k;j<m;j++) S->fi[0][i][j]-=s*S->fi[0][k][j];}
    /* Determinant */
    if (S->fi[0][k][k]<0) sign=-sign;
    *det+=log(fabs(S->fi[0][k][k]));
  }		/* for(k=0;k<m-1;k++) */
  if (S->fi[0][m-1][m-1]<0) sign=-sign;
  *det+=log(fabs(S->fi[0][m-1][m-1]));
  if (*det<=0) *det=sign*exp(*det);
  else *det=sign*(1+*det);
#undef d
  return 0;
}


/* Returns det and cond numbers of a DECOMPOSED matrix */
Entry(Int2) rd_GetMatrixData(MatrPtr mp, FloatHiPtr deteq, FloatHiPtr det) {
  if (deteq && det) return Determinant(mp,deteq,det);
  return -1;
}

/* Transpose a matrix  */
Entry(MatrPtr) rd_TransposeMatrix(MatrPtr from, MatrPtr to DEBPAR) {
}

/* Clears a Matrix */
Entry(void) rd_ClearMatrix(MatrPtr mp) {
}

/* Copy block */
Entry(void) rd_CopyMatrixBlock(MatrPtr from, Int2 frow, Int2 fcol,
                            MatrPtr to, Int2 trow, Int2 tcol,
                            Int2 nrow, Int2 ncol) {
}

#undef MATRSIZE


/*****************************************************/
/* Numerical differentiation by finite differencies  */
typedef EntryPtr(Int2,rdDiffFun,(Vector xx, Vector yy));   /* any differentiable function */

Entry(Int2) rd_der (rdDiffFun fun,
         int n,              /* number of function components */
	 VectPtr x,
	 Var d,
	 MatrPtr A)
{
}

/***********************************/
/* Standard solver and decomposer  */

/*======================================*/
/* Decomposition for progonka procedure */
Entry(Int2) rd_Decomp(MatrPtr S DEBPAR) {
Int2 np,nm,j,i,k;

np=nphase;
nm=nmesh;


/* 1 */
stdalg->CopyMatrix(S->c[1],T);
if (stdalg->Decomp(T)) return 1;
stdalg->SolveMatrix(T,S->a[1],Z1,t);

/* 2 */		   
stdalg->SolveMatrix(T,S->b[1],Z2,t);

/* 3 */
stdalg->MultiplyMatrixMatrix(S->c[0],Z1,T);
stdalg->AddScaledMatrix(S->a[0],T,-1.0,T);
if (stdalg->Decomp(T)) return 3;
stdalg->MultiplyMatrixMatrix(S->c[0],Z2,Z1);
stdalg->AddScaledMatrix(Z1,S->b[0],-1.0,Z1);
stdalg->SolveMatrix(T,Z1,A0,t);

/* 4 */
stdalg->MultiplyMatrixMatrix(S->a[1],A0,T);
stdalg->AddScaledMatrix(T,S->b[1],1.0,T);
if (stdalg->Decomp(T)) return 4;
stdalg->SolveMatrix(T,S->a[1],Z1,t);
stdalg->CopyMatrix(T,S->c[0]);

/* 5 - mimo */

/* 6 */
stdalg->MultiplyMatrixMatrix(S->a[0],A0,T);
stdalg->AddScaledMatrix(T,S->b[0],1.0,S->R);
stdalg->MultiplyMatrixMatrix(S->R,Z1,T);
stdalg->AddScaledMatrix(S->a[0],T,-1.0,S->a[0]);
if (stdalg->Decomp(S->a[0])) return 6;
stdalg->CopyMatrix(A0,S->b[0]);

/* 7 */
stdalg->SolveMatrix(S->c[0],S->delta[1],ZZ2,t);

/* 8 */
stdalg->MultiplyMatrixMatrix(S->R,ZZ2,tt);
stdalg->AddScaledMatrix(tt,S->delta[0],-1.0,ZZ2);
stdalg->SolveMatrix(S->a[0],ZZ2,S->delta[0],t);

/* 9,10,11 */
for(i=1; i<nm-1; i++) {
                     /* 9 */
                     stdalg->MultiplyMatrixMatrix(S->a[i],S->b[i-1],T);
                     stdalg->AddScaledMatrix(T,S->b[i],1.0,T);
                     if(stdalg->Decomp(T)) return 9;
		     for(j=0; j<np; j++){for(k=0; k<np; k++) t[k]=-S->c[i][k][j];
 		                         stdalg->Solve(T,t);
		                         for(k=0; k<np; k++) S->b[i][k][j]=t[k];
                                        }
		     stdalg->CopyMatrix(T,S->c[i]);
		     
		     /* 10 */
                     stdalg->MultiplyMatrixMatrix(S->a[i],S->delta[i-1],tt);
                     stdalg->AddScaledMatrix(tt,S->delta[i],1.0,tt);
		     stdalg->AddScaledMatrix(NULL,tt,-1.0,ZZ2);
 		     stdalg->SolveMatrix(S->c[i],ZZ2,S->delta[i],t);
		     
		     /* 11 - mimo */
                    }		                         
/* 12 */
stdalg->MultiplyMatrixMatrix(S->a[nm-1],S->b[nm-3],T);
stdalg->AddScaledMatrix(T,S->b[nm-1],1.0,S->b[nm-1]);
stdalg->MultiplyMatrixMatrix(S->b[nm-1],S->delta[nm-2],tt);
stdalg->AddScaledMatrix(tt,S->delta[nm-1],1.0,S->delta[nm-1]);
stdalg->MultiplyMatrixMatrix(S->a[nm-1],S->delta[nm-3],tt);
stdalg->AddScaledMatrix(tt,S->delta[nm-1],1.0,S->delta[nm-1]);
stdalg->AddScaledMatrix(NULL,S->delta[nm-1],-1.0,S->delta[nm-1]);
stdalg->MultiplyMatrixMatrix(S->b[nm-1],S->b[nm-2],T);
stdalg->AddScaledMatrix(T,S->c[nm-1],1.0,S->c[nm-1]);
if (stdalg->Decomp(S->c[nm-1])) return 12;
stdalg->SolveMatrix(S->c[nm-1],S->delta[nm-1],tt,t);
stdalg->CopyMatrix(tt,S->delta[nm-1]);
    
/* 13 - mimo */

/* 14 - mimo */

/* 15 */
for(i=nm-2; i>-1; i--) {stdalg->MultiplyMatrixMatrix(S->b[i],S->delta[i+1],tt);
                        stdalg->AddScaledMatrix(tt,S->delta[i],1.0,S->delta[i]);
                       }	       	   
/* 16 */
for(i=0; i<nm; i++) {stdalg->MultiplyMatrixMatrix(S->p[i],S->delta[i],ttfi);
                     stdalg->AddScaledMatrix(S->fi[0],ttfi,1.0,S->fi[0]);
		    }
if (stdalg->Decomp(S->fi[0])) return 16;

/* 17 mimo */

return 0;
}

/*======================================*/
/* Solving for progonka procedure */
Entry(void) rd_Solve(MatrPtr S, VectPtr F) {
Int2 i,np,nm;

np=nphase;
nm=nmesh;

/* 1 - mimo */

/* 2 - mimo */
		   
/* 3 - mimo */

/* 4 - mimo */

/* 5 */
stdalg->CopyVectorElements(F,1*np,ttt,0,np);
stdalg->Solve(S->c[0],ttt);

/* 6 */
stdalg->MultiplyMatrixVector(S->R,ttt,t);
stdalg->CopyVectorElements(F,0*np,ttt,0,np);
stdalg->AddScaledVector(ttt,t,-1.0,ttt);
stdalg->Solve(S->a[0],ttt);
stdalg->CopyVectorElements(ttt,0,F,0*np,np);

/* 7 - mimo */

/* 8 - mimo */

/* 9 - mimo */

/* 10 - mimo */

/* 11 */
for(i=1; i<nm-1; i++) {
                     stdalg->MultiplyMatrixArray(S->a[i],F+(i-1)*np,t);
		     stdalg->CopyVectorElements(F,i*np,ttt,0,np);
                     stdalg->AddScaledVector(ttt,t,-1.0,ttt);
                     stdalg->Solve(S->c[i],ttt);
		     stdalg->CopyVectorElements(ttt,0,F,i*np,np);
		    }         
/* 12 - mimo */

/* 13 */
stdalg->CopyVectorElements(F,(nm-2)*np,ttt,0,np);
stdalg->MultiplyMatrixVector(S->b[nm-1],ttt,t);
stdalg->CopyVectorElements(F,(nm-1)*np,ttt,0,np);
stdalg->AddScaledVector(ttt,t,-1.0,t);
stdalg->CopyVectorElements(F,(nm-3)*np,ttt,0,np);
stdalg->MultiplyMatrixVector(S->a[nm-1],ttt,tttt);
stdalg->AddScaledVector(t,tttt,-1.0,ttt);
stdalg->Solve(S->c[nm-1],ttt);
stdalg->CopyVectorElements(ttt,0,F,(nm-1)*np,np);

/* 14 */
for(i=nm-2; i>-1; i--) {
			stdalg->MultiplyMatrixArray(S->b[i],F+(i+1)*np,t);
                        stdalg->CopyVectorElements(F,i*np,ttt,0,np);
                        stdalg->AddScaledVector(t,ttt,1.0,ttt);
			stdalg->CopyVectorElements(ttt,0,F,i*np,np);
                       }	       	   
/* 15 - mimo */

/* 16 */
MemMove(beta,F+nm*np,sizeof(VAR)*(nappend+1));
for(i=0; i<nm; i++) {
                       stdalg->MultiplyMatrixArray(S->p[i],F+i*np,beta1);
		       stdalg->AddScaledVector(beta,beta1,-1.0,beta);
                      }		      
stdalg->Solve(S->fi[0],beta);
stdalg->CopyVectorElements(beta,0,F,nm*np,nappend+1);

/* 17 */
for(i=0; i<nm; i++) {
                       stdalg->MultiplyMatrixVector(S->delta[i],beta,t);
                       stdalg->CopyVectorElements(F,i*np,ttt,0,np);
                       stdalg->AddScaledVector(ttt,t,1.0,ttt);
                       stdalg->CopyVectorElements(ttt,0,F,i*np,np);
		      }         
}
/*======================================*/


Entry(Int2) rd_SngSolve (MatrPtr AA, FloatHi EPS, VectPtr XX) {
}


Entry(Int2) rd_CharPolyCoeffs(MatrPtr AA, VectPtr PP) {
}

Entry(Int2) rd_dbalan(MatrPtr AA, Int2 *LOW, Int2 *IGH) {
}


Entry(Int2) rd_delmhe(MatrPtr AA, Int2 LOW, Int2 IGH) {
}

Entry(Int2) rd_dhqr(MatrPtr HH, Int2 LOW, Int2 IGH, VectPtr WWR, VectPtr WWI) {
}

Entry(Int2) rd_EigenValues(MatrPtr A, VectPtr Re, VectPtr Im) {
}

Entry(Var) rd_NormOfVector(VectPtr a) {
  return sqrt(stdalg->ScalProd(a,a)/nmesh);
}