source: git/Singular/LIB/control.lib @ 2301b7

version="$Id: control.lib,v 1.37 2007-01-23 15:03:29 Singular Exp $";
category="System and Control Theory";
info="
LIBRARY:  control.lib     Algebraic analysis tools for System and Control Theory

AUTHORS:  Oleksandr Iena       yena@mathematik.uni-kl.de
@*        Markus Becker        mbecker@mathematik.uni-kl.de
@*        Viktor Levandovskyy  levandov@mathematik.uni-kl.de

SUPPORT: Forschungsschwerpunkt 'Mathematik und Praxis' (Project of Dr. E. Zerz
and V. Levandovskyy), Uni Kaiserslautern

MAIN PROCEDURES:
  control(R);     analysis of controllability-related properties of R (using Ext modules)
  controlDim(R);  analysis of controllability-related properties of R (using dimension)
  autonom(R);     analysis of autonomy-related properties of R (using Ext modules)
  autonomDim(R);  analysis of autonomy-related properties of R (using dimension)

COMPONENT PROCEDURES:
  leftKernel(R);       a left kernel of R
  rightKernel(R);      a right kernel of R
  leftInverse(R);      a left inverse of R
  rightInverse(R);     a right inverse of R
  smith(M);            a Smith form of a module M
  colrank(M);          a column rank of M as of matrix
  genericity(M);       analysis of the genericity of parameters
  canonize(L);         Groebnerification for modules in the output of control or autonomy procs
  iostruct(R);         computes an IO-structure of behavior given by a module R
  findTorsion(R, I);   generators of the submodule of a module R, annihilated by the ideal I

AUXILIARY PROCEDURES:
  controlExample(s);   set up an example from the mini database inside of the library
  view();              well-formatted output of lists, modules and matrices
";

LIB "homolog.lib";
LIB "poly.lib";
LIB "primdec.lib";
LIB "matrix.lib";

//---------------------------------------------------------------
static proc Opt_Our()
"USAGE: Opt_Our();
RETURN: intvec, where previous options are stored
PURPOSE: save previous options and set customized options
"
{
  intvec v;
  v=option(get);
  option(redSB);
  option(redTail);
  return (v);
}

//-------------------------------------------------------------------------

static proc space(int n)
"USAGE:spase(n); n is an integer (number of needed spaces)
RETURN:  string consisting of n spaces
NOTE:  the procedure is used in the procedure 'view' to have a better formatted output
"{
  int i;
  string s="";
  for(i=1;i<=n;i++)
  {
    s=s+" ";
  };
  return(s);
};
//-----------------------------------------------------------------------------
proc view(M)
"USAGE:  view(M);   M is of any type
RETURN:  no return value
PURPOSE:  procedure for (well-) formatted output of modules, matrices, lists of modules, matrices; shows everything even if entries are long
NOTE:  in case of other types( not 'module', 'matrix', 'list') works just as standard 'print' procedure
EXAMPLE:  example view; shows an example
"{
  // to be replaced with something more feasible
  if ( (typeof(M)=="module")||(typeof(M)=="matrix") )
  {
  int @R=nrows(M);
  int @C=ncols(M);
  int i;
  int j;
  list MaxLength=list();
  int Size=0;
  int max;
  string s;

  for(i=1;i<=@C;i++)
  {
    max=0;

    for(j=1;j<=@R;j++)
    {
      Size=size( string( M[j,i] ) );
      if( Size>max )
      {
        max=Size;
      };
    };
    MaxLength[i] = max;
  };

  for(i=1;i<=@R;i++)
  {
    s="";
    for(j=1;j<@C;j++)
    {
      s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) ) +",";
    };

    s=s+string(M[i,j])+space( MaxLength[j]-size( string( M[i,j] ) ) );

    if (i!=@R)
    {
      s=s+",";
    };
    print(s);
  };

  return();
  };

  if(typeof(M)=="list")
  {
    int sz=size(M);
    int i;
    for(i=1;i<=sz;i++)
    {
      view(M[i]);
      print("");
    };

    return();
  };
  print(M);
  return();
}
example
{"EXAMPLE:";echo = 2;
  ring r;
  list L;
  matrix M[1][3] = x2+x,y3-y,z5-4z+7;
  L[1] = "a matrix:";
  L[2] = M;
  L[3] = "an ideal:";
  L[4] = ideal(M);
  view(L);
};
//--------------------------------------------------------------------------
proc rightKernel(matrix M)
"USAGE:  rightKernel(M);   M a matrix
RETURN:  module
PURPOSE: computes the right kernel of matrix M (a module of all elements v such that Mv=0)
EXAMPLE: example rightKernel; shows an example
"{
  return(modulo(M,std(0)));
}
example
{"EXAMPLE:";echo = 2;
  ring r = 0,(x,y,z),dp;
  matrix M[1][3] = x,y,z;
  print(M);
  matrix R = rightKernel(M);
  print(R);
  // check:
  print(M*R);
};
//-------------------------------------------------------------------------
proc leftKernel(matrix M)
"USAGE:  leftKernel(M);   M a matrix
RETURN:  module
PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0)
EXAMPLE:  example leftKernel; shows an example
"
{
  return( transpose( modulo( transpose(M),std(0) ) ) );
}
example
{"EXAMPLE:";echo = 2;
  ring r= 0,(x,y,z),dp;
  matrix M[3][1] = x,y,z;
  print(M);
  matrix L = leftKernel(M);
  print(L);
  // check:
  print(L*M);
};
//------------------------------------------------------------------------
proc leftInverse(module M)
"USAGE:  leftInverse(M);  M a module
RETURN:  module
PURPOSE: computes such a matrix L, that LM = Id;
EXAMPLE: example leftInverse; shows an example
NOTE: exists only in the case when M is free submodule
"
{
  // it works also for the NC case;
  int NCols = ncols(M);
  module Id = freemodule(NCols);
  module N  = transpose(M);
  intvec old_opt=Opt_Our();
  Id = std(Id);
  matrix T;
  // check the correctness (Id \subseteq M)
  // via dimension: dim (M) = -1!
  int d = dim_Our(N);
  if (d != -1)
  {
    // No left inverse exists
    return(matrix(0));
  }
  matrix T2 = lift(N, Id);
  T2 =  transpose(T2);
  option(set,old_opt);  // set the options back
  return(T2);
}
example
{
  "EXAMPLE:";echo =2;
  // a trivial example:
  ring r = 0,(x,z),dp;
  matrix M[2][1] = 1,x2z;
  print(M);
  print( leftInverse(M) );
  kill r;
  // derived from the example TwoPendula:
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  matrix U[3][1];
  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
  U[2,1]=(-L1)*Dt^4+(g)*Dt^2;
  U[3,1]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
  module M = module(U);
  module L = leftInverse(M);
  print(L);
  // check
  print(L*M);
};
//-----------------------------------------------------------------------
proc rightInverse(module R)
"USAGE:  rightInverse(M);  M a module
RETURN:  module
PURPOSE: computes such a matrix L, that ML = Id
EXAMPLE: example rightInverse; shows an example
NOTE: exists only in the case when M is free submodule
"
{
  return(transpose(leftInverse(transpose(R))));
}
example
{ "EXAMPLE:";echo =2;
  // a trivial example:
  ring r = 0,(x,z),dp;
  matrix M[1][2] = 1,x2+z;
  print(M);
  print( rightInverse(M) );
  kill r;
  // derived from the TwoPendula example:
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  matrix U[1][3];
  U[1,1]=(-L2)*Dt^4+(g)*Dt^2;
  U[1,2]=(-L1)*Dt^4+(g)*Dt^2;
  U[1,3]=(L1*L2)*Dt^4+(-g*L1-g*L2)*Dt^2+(g^2);
  module M = module(U);
  module L = rightInverse(M);
  print(L);
  // check
  print(M*L);
};
//-----------------------------------------------------------------------
static proc dim_Our(module R)
{
  int d;
  if (attrib(R,"isSB")<>1)
  {
    R = std(R);
  }
  d = dim(R);
  return(d);
}
//-----------------------------------------------------------------------
static proc Ann_Our(module R)
{
  return(Ann(R));
}
//-----------------------------------------------------------------------
static proc Ext_Our(int i, module R, list #)
{
  // mimicking 'Ext_R' from homolog.lib
  int ExtraArg = ( size(#)>0 );
  if (ExtraArg)
  {
    return( Ext_R(i,R,#[1]) );
  }
  else
  {
    return( Ext_R(i,R) );
  }
}
//------------------------------------------------------------------------
static proc is_zero_Our
{
  //just a copy of 'is_zero' from "poly.lib" patched with GKdim
  int d = dim_Our(std(#[1]));
  int a = ( d==-1 );
  if( size(#) >1 ) { list L=a,d; return(L); }
  return(a);
  //  return( is_zero(R) ) ;
};
//------------------------------------------------------------------------
static proc control_output(int i, int NVars, module R, module Ext_1, list Gen)
"USAGE:  control_output(i, NVars, R, Ext_1),
PURPOSE: where
@*         i is integer (number of first nonzero Ext or a number of variables in a basering + 1 in case that all the Exts are zero),
@*           NVars:  integer, number of variables in a base ring,
@*           R:  module R (cokernel representation),
@*           Ext_1:  module, the first Ext(its cokernel representation)
RETURN:  list with all the contollability properties of the system which is to be returned in 'control' procedure
NOTE:  this procedure is used in 'control' procedure
"{
  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
  // Is not clear what to do with gl.dim of qrings
  string DofS = "dimension of the system:";
  string Fn   = "number of first nonzero Ext:";
  string Gen_mes = "Parameter constellations which might lead to a non-controllable system:";

  module RK = rightKernel(R);
  int d=dim_Our(std(transpose(R)));

  if (i==1)
  {
    return(
            list ( Fn,
                   i,
                  "not controllable , image representation for controllable part:",
                    RK,
                  "kernel representation for controllable part:",
                   leftKernel( RK ),
                  "obstruction to controllability",
                   Ext_1,
                  "annihilator of torsion module (of obstruction to controllability)",
                   Ann_Our(Ext_1),
                   DofS,
                   d
                 )
          );
  };

  if(i>NVars)
  {
    return( list(  Fn,
                   -1,
                  "strongly controllable(flat), image representation:",
                   RK,
                  "left inverse to image representation:",
                   leftInverse(RK),
                   DofS,
                   d,
                   Gen_mes,
                   Gen)
          );
  };

  //
  //now i<=NVars
  //

  if( (i==2) )
  {
    return( list( Fn,
                  i,
                 "controllable, not reflexive, image representation:",
                  RK,
                  DofS,
                  d,
                  Gen_mes,
                  Gen)
          );
  };

  if( (i>=3) )
  {
    return( list ( Fn,
                   i,
                  "reflexive, not strongly controllable, image representation:",
                   RK,
                   DofS,
                   d,
                      Gen_mes,
                   Gen)
          );
  };
};
//-------------------------------------------------------------------------

proc control(module R)
"USAGE:  control(R);  R a module (R is the matrix of the system of equations to be investigated)
RETURN:  list
PURPOSE: compute the list of all the properties concerning controllability of the system (behavior), represented by the matrix R
NOTE: the properties and corresponding data like controllability, flatness, dimension of the system, degree of controllability, kernel and image representations, genericity of parameters, obstructions to controllability, annihilator of torsion submodule and left inverse are investigated
EXAMPLE:  example control; shows an example
"
{
  int i;
  int NVars=nvars(basering);
  // TODO: NVars to be replaced with the global hom. dimension of basering!!!
  int ExtIsZero;
  intvec v=Opt_Our();
  module R_std=std(R);
  module Ext_1 = std(Ext_Our(1,R_std));

  ExtIsZero=is_zero_Our(Ext_1);
  i=1;
  while( (ExtIsZero) && (i<=NVars) )
  {
    i++;
    ExtIsZero = is_zero_Our( Ext_Our(i,R_std) );
  };
  matrix T=lift(R,R_std);
  list l=genericity(T);
  option(set,v);

  return( control_output( i, NVars, R, Ext_1, l ) );
}
example
{"EXAMPLE:";echo = 2;
  // a WindTunnel example
  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
  module R;
  R = [D1+a, -k*a*delta, 0, 0],
      [0, D1, -1, 0],
      [0, omega^2, D1+2*zeta*omega, -omega^2];
  R=transpose(R);
  view(R);
  view(control(R));
};
//--------------------------------------------------------------------------
proc controlDim(module R)
"USAGE:  controlDim(R); R a module (R is the matrix of the system of equations to be investigated)
RETURN: list
PURPOSE: computes list of all the properties concerning controllability of the system (behavior), represented by the  matrix R
EXAMPLE:  example controlDim; shows an example
NOTE: the properties and corresponding data like controllability, flatness, dimension of the system, degree of controllability, kernel and image representations, genericity of parameters, obstructions to controllability, annihilator of torsion submodule and left inverse are investigated.
@* This procedure is analogous to 'control' but uses dimension calculations.
@* The implemented approach works for full row rank matrices only (the check is done automatically).
"
{
  if( nrows(R) != colrank(transpose(R)) )
  {
    return ("controlDim cannot be applied, since R does not have full row rank");
  }
  intvec     v = Opt_Our();
  module R_std = std(R);
  int        d = dim_Our(R_std);
  int    NVars = nvars(basering);
  int        i = NVars-d;
  module Ext_1 = std(Ext_Our(1,R_std));
  matrix     T = lift(R,R_std);
  list       l = genericity(T);
  option(set, v);
  return( control_output( i, NVars, R, Ext_1, l));
}
example
{"EXAMPLE:";echo = 2;
  //a WindTunnel example
  ring A = (0,a, omega, zeta, k),(D1, delta),dp;
  module R;
  R = [D1+a, -k*a*delta, 0, 0],
      [0, D1, -1, 0],
      [0, omega^2, D1+2*zeta*omega, -omega^2];
  R=transpose(R);
  view(R);
  view(controlDim(R));
};
//------------------------------------------------------------------------
proc colrank(module M)
"USAGE: colrank(M); M a matrix/module
RETURN: int
PURPOSE: compute the column rank of M as of matrix
NOTE: this procedure uses Bareiss algorithm
"{
  // NOte continued:
  // which might not terminate in some cases
  module M_red  = bareiss(M)[1];
  int NCols_red = ncols(M_red);
  return (NCols_red);
}
example
{"EXAMPLE: ";echo = 2;
  // de Rham complex
  ring r=0,(D(1..3)),dp;
  module R;
  R=[0,-D(3),D(2)],
    [D(3),0,-D(1)],
    [-D(2),D(1),0];
  R=transpose(R);
  colrank(R);
};

//------------------------------------------------------------------------
static proc autonom_output( int i, int NVars, module RC, int R_rank )
"USAGE:  proc autonom_output(i, NVars, RC, R_rank)
           i:  integer, number of first nonzero Ext or
               just number of variables in a base ring + 1 in case that all the Exts are zero
           NVars:  integer, number of variables in a base ring
           RC: module, kernel-representation of controllable part of the system
           R_rank: integer, column rank of the representation matrix
PURPOSE: compute all the autonomy properties of the system which is to be returned in 'autonom' procedure
RETURN:  list
NOTE:  this procedure is used in 'autonom' procedure
"
{
  int d=NVars-i;//that is the dimension of the system
  string DofS="dimension of the system:";
  string Fn = "number of first nonzero Ext:";
  if(i==0)
  {
    return( list(  Fn,
                   i,
                  "not autonomous",
                   "kernel representation for controllable part",
                   RC,
                   "column rank of the matrix",
                   R_rank,
                   DofS,
                   d )
          );
  };

  if( i>NVars )
  {
    return( list( Fn,
                  -1,
                  "trivial",
                  DofS,
                  d )
          );
  };

  //
  //now i<=NVars
  //


  if( i==1 )
  // in case that NVars==1 there is no sense to consider the notion
  // of strongly autonomous behavior, because it does not imply
  // that system is overdetermined in this case
  {
    return( list ( Fn,
                   i,
                  "autonomous, not overdetermined",
                   DofS,
                   d )
          );
  };

  if( i==NVars )
  {
    return( list(  Fn,
                   i,
                  "strongly autonomous(fin. dimensional),in particular overdetermined",
                   DofS,
                   d)
          );
  };

  if( i<NVars )
  {
    return( list ( Fn,
                   i,
                  "overdetermined, not strongly autonomous",
                   DofS,
                   d)
          );
  };
};
//--------------------------------------------------------------------------
proc autonomDim(module R)
"USAGE:  autonomDim(R);   R a module (R is a matrix of the system of equations which is to be investigated)
RETURN: list
PURPOSE: computes the list of all the properties concerning autonomy of the system (behavior), represented by the matrix R
NOTE: the properties and corresponding data like autonomy resp. strong autonomy, dimension of the system, autonomy degree, kernel representation and (over)determinacy are investigated.
@* This procedure is analogous to 'autonom' but uses dimension calculations
EXAMPLE:  example autonomDim; shows an example
"
{
  int d;
  int NVars  = nvars(basering);
  module RT  = transpose(R);
  module RC;  //for computation of controllable part if if exists
  int R_rank = ncols(R);
  d          = dim_Our( std(RT) );  //this is the dimension of the system
  int      i = NVars-d;  //First non-zero Ext
  if( d==0 )
    {
      RC = leftKernel(rightKernel(R));
      R_rank=colrank(R);
    }
  return( autonom_output(i,NVars,RC,R_rank) );
}
example
{"EXAMPLE:"; echo = 2;
  // Cauchy1 example
  ring r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
            [s3,-s4],
            [s4, s3];
  R=transpose(R);
  view( R );
  view( autonomDim(R) );
};
//----------------------------------------------------------
proc autonom(module R)
"USAGE:  autonom(R);   R a module (R is a matrix of the system of equations which is to be investigated)
RETURN:  list
PURPOSE: find all the properties concerning autonomy of the system (behavior) represented by the  matrix R
NOTE: the properties and corresponding data like autonomy resp. strong autonomy, dimension of the system, autonomy degree, kernel representation and (over)determinacy are investigated
EXAMPLE: example autonom; shows an example
"
{
  int NVars=nvars(basering);
  int ExtIsZero;
  module RT=transpose(R);
  module RC;
  int R_rank=ncols(R);
  ExtIsZero=is_zero_Our(Ext_Our(0,RT));
  int i=0;
  while( (ExtIsZero)&&(i<=NVars) )
  {
    i++;
    ExtIsZero = is_zero_Our(Ext_Our(i,RT));
  };
  if (i==0)
  {
    RC = leftKernel(rightKernel(R));
    R_rank=colrank(R);
  }
  return(autonom_output(i,NVars,RC,R_rank));
}
example
{"EXAMPLE:"; echo = 2;
  // Cauchy
  ring r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
            [s3,-s4],
            [s4, s3];
  R=transpose(R);
  view( R );
  view( autonom(R) );
};


//----------------------------------------------------------
proc genericity(matrix M)
"USAGE:  genericity(M); M is a matrix/module
RETURN:  list (of strings)
PURPOSE: determine parametric expressions which have been assumed to be non-zero in the process of computing the Groebner basis
NOTE: the output list consists of strings. The first string contains the variables only, whereas each further string contain a single polynomial in parameters.
$* We strongly recommend to switch on the redSB and redTail options.
@* The procedure is effective with the lift procedure for modules with parameters
EXAMPLE:  example genericity; shows an example
"
{
  // returns "-", if there are no parameters!
  if (npars(basering)==0)
  {
    return("-");
  }
  list RT = evas_genericity(M); // list of strings
  if ((size(RT)==1) && (RT[1] == ""))
  {
    return("-");
  }
  return(RT);
}
example
{  // TwoPendula
  "EXAMPLE:"; echo = 2;
  ring r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  module RR =
    [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
    [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
    [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];
  module R = transpose(RR);
  module SR = std(R);
  matrix T = lift(R,SR);
  genericity(T);
  //-- The result might be different when computing reduced bases:
  matrix T2;
  option(redSB);
  option(redTail);
  module SR2 = std(R);
  T2 =  lift(R,SR2);
  genericity(T2);
}
//---------------------------------------------------------------
static proc victors_genericity(matrix M)
{
  // returns "-", if there are no parameters!
  if (npars(basering)==0)
  {
    return("-");
  }
  int plevel = printlevel-voice+2;
  // M is a matrix over a ring with params and vars;
  ideal I = ideal(M); // a list of entries
  I = simplify(I,2); // delete 0's
  // decompose every coeff in every poly
  int i;
  int s = size(I);
  ideal NM;
  poly p;
  number num;
  int  cl=1;
  intvec ZeroVec; ZeroVec[nvars(basering)] = 0;
  intvec W;
  ideal Numero, Denomiro;
  int cNu=0; int cDe=0;
  for (i=1; i<=s; i++)
  {
    // remove contents and add them as polys
    p   = I[i];
    W   = leadexp(p);
    if (W == ZeroVec) // i.e. just a coef
    {
      num    = denominator(leadcoef(p)); // from poly.lib
      NM[cl] = numerator(leadcoef(p));
      dbprint(p,"numerator:");
      dbprint(p, string(NM[cl]));
      cNu++; Numero[cNu]= NM[cl];
      cl++;
      NM[cl] = num; // denominator
      dbprint(p,"denominator:");
      dbprint(p, string(NM[cl]));
      cDe++; Denomiro[cDe]= NM[cl];
      cl++;
      p = p - lead(p); // for the next cycle
    }
    if ( p!= 0)
    {
      num = content(p);
      p   = p/num;
      NM[cl] = denominator(num);
      dbprint(p,"content denominator:");
      dbprint(p, string(NM[cl]));
      cNu++; Numero[cNu]= NM[cl];
      cl++;
      NM[cl] = numerator(num);
      dbprint(p,"content numerator:");
      dbprint(p, string(NM[cl]));
      cDe++; Denomiro[cDe]= NM[cl];
      cl++;
    }
    // it seems that the next elements will not have real influence
    while( p != 0)
    {
      NM[cl] = leadcoef(p); // should be all integer, i.e. non-rational
      dbprint(p,"coef:");
      dbprint(p, string(NM[cl]));
      cl++;
      p = p - lead(p);
    }
  }
  NM = simplify(NM,4); // delete identical
  string newvars = parstr(basering);
  def save = basering;
  string NewRing = "ring @NR =" +string(char(basering))+",("+newvars+"),Dp;";
  execute(NewRing);
  // get params as variables
  // create a list of non-monomials
  ideal @L;
  ideal F;
  ideal NM = imap(save,NM);
  NM = simplify(NM,8); //delete multiples
  poly p,q;
  cl = 1;
  int j, cf;
  for(i=1; i<=size(NM);i++)
  {
    p = NM[i] - lead(NM[i]);
    if (p!=0)
    {
      //      L[cl] = p;
      F  = factorize(NM[i],1); //non-constant factors only
      cf = 1;
      // factorize every polynomial
      // throw away every monomial from factorization (also constants from above ring)
      for (j=1; j<=size(F);j++)
      {
        q = F[j]-lead(F[j]);
        if (q!=0)
        {
          @L[cl] = F[j];
          cl++;
        }
      }
    }
  }
  // return the result [in string-format]
  @L = simplify(@L,2+4+8); // skip zeroes, doubled and entries, diff. by a constant
  list SL;
  for (j=1; j<=size(@L);j++)
  {
    SL[j] = string(@L[j]);
  }
  setring save;
  return(SL);
}
//---------------------------------------------------------------
static proc evas_genericity(matrix M)
{
  // called from the main genericity proc
  ideal I = ideal(M);
  I = simplify(I,2+4);
  int s = ncols(I);
  ideal Den;
  poly p;
  int i;
  for (i=1; i<=s; i++)
  {
    p = I[i];
    while (p !=0)
    {
      Den = Den, denominator(leadcoef(p));
      p   = p-lead(p);
    }
  }
  Den = simplify(Den,2+4);
  string newvars = parstr(basering);
  def save = basering;
  string NewRing = "ring @NR =(" +string(char(basering))+"),("+newvars+"),Dp;";
  execute(NewRing);
  ideal F;
  ideal Den = imap(save,Den);
  Den = simplify(Den,2);
  int s1 = ncols(Den);
  for (i=1; i<=s1; i++)
  {
    if (normalize(Den[i]) !=1)
    {
      F= F, factorize(Den[i],1);
    }
  }
  F = simplify(F, 2+4+8);
  ideal @L = F;
  @L = simplify(@L,2);
  list SL;
  int c,j;
  string Mono;
  c = 1;
  for (j=1; j<= ncols(@L);j++)
  {
    if (leadcoef(@L[j]) <0)
    {
      @L[j] = -1*@L[j];
    }
    if (((@L[j] - lead(@L[j]))==0 ) && (@L[j]!=0) ) //@L[j] is a monomial
    {
      Mono = Mono + string(@L[j])+ ","; // concatenation
    }
    else
    {
      c++;
      SL[c] = string(@L[j]);
    }
  }
  if (Mono!="")
  {
    Mono = Mono[1..size(Mono)-1]; // delete the last colon
  }
  SL[1] = Mono;
  setring save;
  return(SL);
}

//---------------------------------------------------------------
proc canonize(list L)
"USAGE:  canonize(L); L a list
RETURN:  list
PURPOSE: modules in the list are canonized by computing their reduced minimal (= unique up to constant factor w.r.t. the given ordering) Groebner bases
ASSUME:  L is the output of control/autonomy procedures
EXAMPLE:  example canonize; shows an example
"
{
  list M = L;
  intvec v=Opt_Our();
  int s = size(L);
  int i;
  for (i=2; i<=s; i=i+2)
  {
    if (typeof(M[i])=="module")
    {
      M[i] = std(M[i]);
      //      M[i] = prune(M[i]); // mimimal embedding: no need yet
      //      M[i] = std(M[i]);
    }
  }
  option(set, v); //set old values back
  return(M);
}
example
{  // TwoPendula with L1=L2=L
  "EXAMPLE:"; echo = 2;
  ring r=(0,m1,m2,M,g,L),Dt,dp;
  module RR =
    [m1*L*Dt^2, m2*L*Dt^2, -1, (M+m1+m2)*Dt^2],
    [m1*L^2*Dt^2-m1*L*g, 0, 0, m1*L*Dt^2],
    [0, m2*L^2*Dt^2-m2*L*g, 0, m2*L*Dt^2];
  module R = transpose(RR);
  list C = control(R);
  list CC = canonize(C);
  view(CC);
}

//----------------------------------------------------------------

static proc elementof (int i, intvec v)
{
  int b=0;
  for(int j=1;j<=nrows(v);j++)
    {
      if(v[j]==i)
        {
          b=1;
          return (b);
        }
    }
  return (b);
}
//-----------------------------------------------------------------
proc iostruct(module R)
"USAGE: iostruct( R ); R a module
RETURN:  list L with entries: string s, intvec v, module P and module Q
PURPOSE:  if R is the kernel-representation-matrix of some system, then we output a input-ouput representation Py=Qu of the system, the components that have been chosen as outputs(intvec v) and a comment s
NOTE:  the procedure uses Bareiss algorithm
EXAMPLE:  example iostruct; shows an example
"
{
  // NOTE cont'd
  //which might not terminate in some cases
  list L = bareiss(R);
  int R_rank = ncols(L[1]);
  int NCols=ncols(R);
  intvec v=L[2];
  int temp;
  int NRows=nrows(v);
  int i,j;
  int b=1;
  module P;
  module Q;
  int n=0;

  while(b==1)               //sort v through bubblesort
    {
      b=0;
      for(i=1;i<NRows;i++)
        {
          if(v[i]>v[i+1])
          {
            temp=v[i];
            v[i]=v[i+1];
            v[i+1]=temp;
            b=1;
          }
        }
    }
  P=R[v];                     //generate P
  for(i=1;i<=NCols;i++)       //generate Q
    {
      if(elementof(i,v)==1)
        {
          i++;
          continue;
        }
      Q=Q,R[i];
    }
  Q=simplify(Q,2);
  string s="The following components have been chosen as outputs: ";
  return (list(s,v,P,Q));
}
example
{"EXAMPLE:";echo = 2;
 //Example Antenna
 ring r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), (c,dp);
 module RR;
 RR =
   [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
   [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
   [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
   [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
   [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
   [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];
 module R = transpose(RR);
 view(iostruct(R));
};

//---------------------------------------------------------------
static proc smdeg(matrix N)
"USAGE: smdeg( N ); N a matrix
RETURN:  intvec
PURPOSE: returns an intvec of length 2 with the index of an element of N with smallest degree
"
{
  int n = nrows(N);
  int m = ncols(N);
  int d,d_temp;
  intvec v;
  int i,j;            // counter

  if (N==0)
  {
    v = 1,1;
    return(v);
  }

  for (i=1; i<=n; i++)
// hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet
  {
    for (j=1; j<=m; j++)
    {
      if( deg(N[i,j])!=-1 )
      {
        d=deg(N[i,j]);
        break;
      }
    }
    if (d != -1)
    {
      break;
    }
  }
  for(i=1; i<=n; i++)
  {
    for(j=1; j<=m; j++)
    {
      d_temp = deg(N[i,j]);
      if ( (d_temp < d) && (N[i,j]!=0) )
      {
        d=d_temp;
      }
    }
  }
  for (i=1; i<=n; i++)
  {
    for (j=1; j<=m;j++)
    {
      if ( (deg(N[i,j]) == d) && (N[i,j]!=0) )
      {
        v = i,j;
        return(v);
      }
    }
  }
}
//---------------------------------------------------------------
static proc NoNon0Pol(vector v)
"USAGE: NoNon0Pol(v), v a vector
RETURN:  int
PURPOSE: returns 1, if there is only one non-zero element in v and 0 else
"{
  int i,j;
  int n = nrows(v);
  for( j=1; j<=n; j++)
  {
    if (v[j] != 0)
    {
      i++;
    }
  }
  if ( i!=1 )
  {
    i=0;
  }
  return(i);
}
//---------------------------------------------------------------
static proc extgcd_Our(poly p, poly q)
{
  ideal J;   //for extgcd-computations
  matrix T; //----------"------------
  list L;
  // the extgcd-command has a bug in versions before 2-0-7
  if ( system("version")<=2006 )
  {
    J = p,q; // J = N[k-1,k-1],N[k,k]; //J is of type ideal
    L[1] = liftstd(J,T);  //T is of type matrix
    if(J[1]==p) //this is just for the case the SINGULAR swaps the
    //      two elements due to ordering
    {
      L[2] = T[1,1];
      L[3] = T[2,1];
    }
    else
    {
      L[2] = T[2,1];
      L[3] = T[1,1];
    }
  }
  else
  {
    L=extgcd(p,q);
    //    L=extgcd(N[k-1,k-1],N[k,k]);
    //one can use this line if extgcd-bug is fixed
  }
  return(L);
}
static proc normalize_Our(matrix N, matrix Q)
"USAGE: normalize_Our(N,Q), N, Q are two matrices
PURPOSE: normalizes N and divides the columns of Q through the leading coefficients of the columns of N
RETURN: normalized matrix N and altered Q(according to the scheme mentioned in purpose). If number of columns of N and Q do not coincide, N and Q are returned unchanged
NOTE: number of columns of N and Q must coincide.
"
{
  if(ncols(N) != ncols(Q))
    {
      return (N,Q);
    }
  module M = module(N);
  module S = module(Q);
  int NCols = ncols(N);
  number n;
  for(int i=1;i<=NCols;i++)
    {
      n = leadcoef(M[i]);
      if( n != 0 )
        {
          M[i]=M[i]/n;
          S[i]=S[i]/n;
        }
     }
   N = matrix(M);
   Q = matrix(S);
   return (N,Q);
}

//---------------------------------------------------------------
proc smith( module M )
"USAGE: smith(M); M a module/matrix
PURPOSE: computes the Smith normal form of a matrix
RETURN: a list of length 4 with the following entries:
@*      [1]: the Smith normal form S of M,
@*      [2]: the rank of M,
@*      [3]: a unimodular matrix U,
@*      [4]: a unimodular matrix V,
such that U*M*V=S. An warning is returned when no Smith form exists.
NOTE: The Smith form only exists over PIDs (principal ideal domains). Use global ordering for computations!
"
{
  if (nvars(basering)>1) //if more than one variable, return empty list
  {
    string s="The Smith-Form only exists for principal ideal domains";
    return (s);
  }
  matrix N = matrix(M);         //Typecasting
  int n = nrows(N);
  int m = ncols(N);
  matrix P = unitmat(n);       //left transformation matrix
  matrix Q = unitmat(m);       //right transformation matrix
  int k, i, j, deg_temp;
  poly tmp;
  vector v;
  list L;                      //for extgcd-computation
  intmat f[1][n];              //to save degrees
  matrix lambda[1][n];         //to save leadcoefficients
  intmat g[1][m];              //to save degrees
  matrix mu[1][m];             //to save leadcoefficients
  int ii;                       //counter

  while ((k!=n) && (k!=m) )
  {
    k++;
    while ((k<=n) && (k<=m))  //outer while-loop for column-operations
    {
      while(k<=m )        //inner while-loop for row-operations
      {
        if( (n>m) && (k < n) && (k<m))
        {
          if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0)
          {
            return(N,k-1,P,Q);
          }
        }
        i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix
        i=i+(k-1);                        //indices adjusted to the whole matrix
        j=j+(k-1);
        if(i!=k)                    //take the element with smallest degree in the first position
        {
          N=permrow(N,i,k);
          P=permrow(P,i,k);
        }
        if(j!=k)
        {
          N=permcol(N,j,k);
          Q=permcol(Q,j,k);
        }
        if(NoNon0Pol(N[k])==1)
        {
          break;
        }
        tmp=leadcoef(N[k,k]);
        deg_temp=ord(N[k,k]);             //ord outputs the leading degree of N[k,k]
        for(ii=k+1;ii<=n;ii++)
        {
          lambda[1,ii]=leadcoef(N[ii,k])/tmp;
          f[1,ii]=deg(N[ii,k])-deg_temp;
        }
        for(ii=k+1;ii<=n;ii++)
        {
          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
          P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
          N,Q=normalize_Our(N,Q);
        }
      }
      if (k>n)
      {
        break;
      }
      if(NoNon0Pol(transpose(N)[k])==1)
      {
        break;
      }
      tmp=leadcoef(N[k,k]);
      deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k][k]

      for(ii=k+1;ii<=m;ii++)
      {
        mu[1,ii]=leadcoef(N[k,ii])/tmp;
        g[1,ii]=deg(N[k,ii])-deg_temp;
      }
      for(ii=k+1;ii<=m;ii++)
      {
        N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii);
        Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii);
        N,Q=normalize_Our(N,Q);
      }
    }
    if( (k!=1) && (k<n) && (k<m) )
    {
      L = extgcd_Our(N[k-1,k-1],N[k,k]);
      if ( N[k-1,k-1]!=L[1] )  //means that N[k-1,k-1] is not a divisor of N[k,k]
      {
        N=addrow(N,k-1,L[2],k);
        P=addrow(P,k-1,L[2],k);
        N,Q=normalize_Our(N,Q);

        N=addcol(N,k,-L[3],k-1);
        Q=addcol(Q,k,-L[3],k-1);
        N,Q=normalize_Our(N,Q);
        k=k-2;
      }
    }
  }
  if( (k<=n) && (k<=m) )
  {
    if( N[k,k]==0)
    {
      return(N,k-1,P,Q);
    }
  }
  return(N,k,P,Q);
}
example
{
  "EXAMPLE:";echo = 2;
  option(redSB);
  option(redTail);
  ring r   = 0,x,dp;
  module M = [x2,x,3x3-4], [2x2-1,4x,5x2], [2x5,3x,4x];
  print(M);
  list P = smith(M);
  print(P[1]);
  matrix N = matrix(M);
  matrix B = P[3]*N*P[4];
  print(B);
}
// see what happens when the matrix is already in Smith-Form
//  module M = [x,0,0],[0,x2,0],[0,0,x3];
//  list L = smith(M);
// print(L[1]);
//matrix N=matrix(M);
//matrix B=L[3]*N*L[4];
//print(B);
//---------------------------------------------------------------
static proc list_tex(L, string name,link l,int nr_loop)
"USAGE: list_tex(L,name,l), where L is a list, name a string, l a link
         writes the content of list L in a tex-file 'name'
RETURN: nothing
"
{
  if(typeof(L)!="list")  //in case L is not a list
  {
    texobj(name,L);
  }
  if(size(L)==0)
  {
  }
  else
  {
    string t;
    for (int i=1;i<=size(L);i++)
    {
      while(1)
      {
        if(typeof(L[i])=="string")  //Fehler hier fuer normalen output->nur wenn string in liste dann verbatim
        {
          t=L[i];
          if(nr_loop==1)
          {
            write(l,"\\begin\{center\}");
            write(l,"\\begin\{verbatim\}");
          }
          write(l,t);
          if(nr_loop==0)
          {
            write(l,"\\par");
          }
          if(nr_loop==1)
          {
            write(l,"\\end\{verbatim\}");
            write(l,"\\end\{center\}");
          }
          break;
        }
        if(typeof(L[i])=="module")
        {
          texobj(name,matrix(L[i]));
          break;
        }
        if(typeof(L[i])=="list")
        {
          list_tex(L[i],name,l,1);
          break;
        }
        write(l,"\\begin\{center\}");
        texobj(name,L[i]);
        write(l,"\\end\{center\}");
        write(l,"\\par");
        break;
      }
    }
  }
}
example
{
  "EXAMPLE:";echo = 2;
}
//---------------------------------------------------------------
proc verbatim_tex(string s, link l)
"USAGE: verbatim_tex(s,l), where s is a string and l a link
PURPOSE: writes the content of s in verbatim-environment in the file
         specified by link
RETURN: nothing
"
{
  write(l,"\\begin{verbatim}");
  write(l,s);
  write(l,"\\end{verbatim}");
  write(l,"\\par");
}
example
{
  "EXAMPLE:";echo = 2;
}
//---------------------------------------------------------------
proc findTorsion(module R, ideal TAnn)
"USAGE:  findTorsion(R, I);   R an ideal/matrix/module, I an ideal
RETURN:  module
PURPOSE: computes the Groebner basis of the submodule of R, annihilated by I
NOTE: especially helpful, when I is the annihilator of the t(R) - the torsion submodule of R. In this case, the result is the explicit presentation of t(R) as
the submodule of R
EXAMPLE: example findTorsion; shows an example
"
{
  // motivation: let R be a module,
  // TAnn is the annihilator of t(R)\subset R
  // compute the generators of t(R) explicitly
  ideal AS = TAnn;
  module S = R;
  if (attrib(S,"isSB")<>1)
  {
    S = std(S);
  }
  if (attrib(AS,"isSB")<>1)
  {
    AS = std(AS);
  }
  int nc  = ncols(S);
  module To = quotient(S,AS);
  To = std(NF(To,S));
  return(To);
}
example
{
  "EXAMPLE:";echo = 2;
  // Flexible Rod
  ring A = 0,(D1, D2), (c,dp);
  module R= [D1, -D1*D2, -1], [2*D1*D2, -D1-D1*D2^2, 0];
  module RR = transpose(R);
  list L = control(RR);
  // here, we have the annihilator:
  ideal LAnn = D1; // = L[10]
  module Tr  = findTorsion(RR,LAnn);
  print(RR);  // the module itself
  print(Tr); // generators of the torsion submodule
}


proc controlExample(string s)
"USAGE:  controlExample(s);   s a string
RETURN:  ring
PURPOSE: set up an example from the mini database by initalizing a ring and a module in a ring
NOTE: in order to see the list of available examples, execute @code{controlExample(\"show\");}
@* To use ab example, one has to do the following. Suppose one calls the ring, where the example will be activated, A. Then, by executing
@*  @code{def A = controlExample(\"Antenna\");} and @code{setring A;},
@* A will become a basering from the example \"Antenna\" with
the predefined system module R (transposed).
After that one can just execute @code{control(R);} respectively
@code{autonom(R);} to perform the control resp. autonomy analysis of R.
EXAMPLE: example controlExample; shows an example
"{
  list E, S, D; // E=official name, S=synonym, D=description
  E[1] = "Cauchy1";  S[1] = "cauchy1";  D[1] = "1-dimensional Cauchy equation";
  E[2] = "Cauchy2";  S[2] = "cauchy2";  D[2] = "2-dimensional Cauchy equation";
  E[3] = "Control1"; S[3] = "control1"; D[3] = "example of a simple noncontrollable system";
  E[4] = "Control2"; S[4] = "control2"; D[4] = "example of a simple controllable system";
  E[5] = "Antenna";  S[5] = "antenna";  D[5] = "antenna";
  E[6] = "Einstein"; S[6] = "einstein"; D[6] = "Einstein equations in vacuum";
  E[7] = "FlexibleRod"; S[7] = "flexible rod"; D[7] = "flexible rod";
  E[8] = "TwoPendula";  S[8] = "two pendula";  D[8] = "two pendula mounted on a cart";
  E[9] = "WindTunnel";  S[9] = "wind tunnel";D[9] = "wind tunnel";
  E[10] = "Zerz1";      S[10] = "zerz1"; D[10] = "example from the lecture of Eva Zerz";
  // all the examples so far
  int i;
  if ( (s=="show") || (s=="Show") )
  {
    print("The list of examples:");
    for (i=1; i<=size(E); i++)
    {
      printf("name: %s,  desc: %s", E[i],D[i]);
    }
    return();
  }
  string t;
  for (i=1; i<=size(E); i++)
  {
    if ( (s==E[i]) || (s==S[i]) )
    {
      t = "def @A = ex"+E[i]+"();";
      execute(t);
      return(@A);
    }
  }
  "No example found";
  return();
}
example
{
  "EXAMPLE:";echo = 2;
  controlExample("show");   // let us see all available examples:
  def B = controlExample("TwoPendula"); // let us set up a particular example
  setring B;
  print(R);
}

//----------------------------------------------------------
//
//Some example rings with defined systems
//----------------------------------------------------------
//autonomy:
//
//----------------------------------------------------------
static proc exCauchy1()
{
  ring @r=0,(s1,s2),dp;
  module R= [s1,-s2],
            [s2, s1];
  R=transpose(R);
  export R;
  return(@r);
}
//----------------------------------------------------------
static proc exCauchy2()
{
  ring @r=0,(s1,s2,s3,s4),dp;
  module R= [s1,-s2],
            [s2, s1],
            [s3,-s4],
            [s4, s3];
  R=transpose(R);
  export R;
  return(@r);
}
//----------------------------------------------------------
static proc exZerz1()
{
  ring @r=0,(d1,d2),dp;
  module R=[d1^2-d2],
           [d2^2-1];
  R=transpose(R);
  export R;
  return(@r);
}
//----------------------------------------------------------
//control
//----------------------------------------------------------
static proc exControl1()
{
  ring @r=0,(s1,s2,s3),dp;
  module R=[0,-s3,s2],
           [s3,0,-s1];
  R=transpose(R);
  export R;
  return(@r);
}
//----------------------------------------------------------
static proc exControl2()
{
  ring @r=0,(s1,s2,s3),dp;
  module R=[0,-s3,s2],
           [s3,0,-s1],
           [-s2,s1,0];
  R=transpose(R);
  export R;
  return(@r);
}
//----------------------------------------------------------
static proc exAntenna()
{
  ring @r = (0, K1, K2, Te, Kp, Kc),(Dt, delta), dp;
  module R;
  R = [Dt, -K1, 0, 0, 0, 0, 0, 0, 0],
      [0, Dt+K2/Te, 0, 0, 0, 0, -Kp/Te*delta, -Kc/Te*delta, -Kc/Te*delta],
      [0, 0, Dt, -K1, 0, 0, 0, 0, 0],
      [0, 0, 0, Dt+K2/Te, 0, 0, -Kc/Te*delta, -Kp/Te*delta, -Kc/Te*delta],
      [0, 0, 0, 0, Dt, -K1, 0, 0, 0],
      [0, 0, 0, 0, 0, Dt+K2/Te, -Kc/Te*delta, -Kc/Te*delta, -Kp/Te*delta];

  R=transpose(R);
  export R;
  return(@r);
}

//----------------------------------------------------------

static proc exEinstein()
{
  ring @r = 0,(D(1..4)),dp;
  module R =
  [D(2)^2+D(3)^2-D(4)^2, D(1)^2, D(1)^2, -D(1)^2, -2*D(1)*D(2), 0, 0, -2*D(1)*D(3), 0, 2*D(1)*D(4)],
  [D(2)^2, D(1)^2+D(3)^2-D(4)^2, D(2)^2, -D(2)^2, -2*D(1)*D(2), -2*D(2)*D(3), 0, 0, 2*D(2)*D(4), 0],
  [D(3)^2, D(3)^2, D(1)^2+D(2)^2-D(4)^2, -D(3)^2, 0, -2*D(2)*D(3), 2*D(3)*D(4), -2*D(1)*D(3), 0, 0],
  [D(4)^2, D(4)^2, D(4)^2, D(1)^2+D(2)^2+D(3)^2, 0, 0, -2*D(3)*D(4), 0, -2*D(2)*D(4), -2*D(1)*D(4)],
  [0, 0, D(1)*D(2), -D(1)*D(2), D(3)^2-D(4)^2, -D(1)*D(3), 0, -D(2)*D(3), D(1)*D(4), D(2)*D(4)],
  [D(2)*D(3), 0, 0, -D(2)*D(3),-D(1)*D(3), D(1)^2-D(4)^2, D(2)*D(4), -D(1)*D(2), D(3)*D(4), 0],
  [D(3)*D(4), D(3)*D(4), 0, 0, 0, -D(2)*D(4), D(1)^2+D(2)^2, -D(1)*D(4), -D(2)*D(3), -D(1)*D(3)],
  [0, D(1)*D(3), 0, -D(1)*D(3), -D(2)*D(3), -D(1)*D(2), D(1)*D(4), D(2)^2-D(4)^2, 0, D(3)*D(4)],
  [D(2)*D(4), 0, D(2)*D(4), 0, -D(1)*D(4), -D(3)*D(4), -D(2)*D(3), 0, D(1)^2+D(3)^2, -D(1)*D(2)],
  [0, D(1)*D(4), D(1)*D(4), 0, -D(2)*D(4), 0, -D(1)*D(3), -D(3)*D(4), -D(1)*D(2), D(2)^2+D(3)^2];

  R=transpose(R);
  export R;
  return(@r);
}

//----------------------------------------------------------
static proc exFlexibleRod()
{
  ring @r = 0,(D1, delta), dp;
  module R;
  R = [D1, -D1*delta, -1], [2*D1*delta, -D1-D1*delta^2, 0];

  R=transpose(R);
  export R;
  return(@r);
}

//----------------------------------------------------------
static proc exTwoPendula()
{
  ring @r=(0,m1,m2,M,g,L1,L2),Dt,dp;
  module R = [m1*L1*Dt^2, m2*L2*Dt^2, -1, (M+m1+m2)*Dt^2],
             [m1*L1^2*Dt^2-m1*L1*g, 0, 0, m1*L1*Dt^2],
             [0, m2*L2^2*Dt^2-m2*L2*g, 0, m2*L2*Dt^2];

  R=transpose(R);
  export R;
  return(@r);
}
//----------------------------------------------------------
static proc exWindTunnel()
{
  ring @r = (0,a, omega, zeta, k),(D1, delta),dp;
  module R = [D1+a, -k*a*delta, 0, 0],
             [0, D1, -1, 0],
             [0, omega^2, D1+2*zeta*omega, -omega^2];

  R=transpose(R);
  export R;
  return(@r);
}