Changeset 501f7c5 in git


Ignore:
Timestamp:
Dec 16, 2004, 5:28:53 PM (20 years ago)
Author:
Viktor Levandovskyy <levandov@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
175c561dc37255410efb4b041fbb41a6e1b0d8b8
Parents:
1d5418b9bdcc715becec9d153116793ed3147e8f
Message:
*levandov: LeftInverse and others rewritten in order to be nc-capable, Smith form routines added


git-svn-id: file:///usr/local/Singular/svn/trunk@7635 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/control.lib

    r1d5418 r501f7c5  
    1 version="$Id: control.lib,v 1.17 2004-12-14 21:05:09 levandov Exp $";
     1version="$Id: control.lib,v 1.18 2004-12-16 16:28:53 levandov Exp $";
    22category="Applications";
    33info="
     
    1818LeftKernel(module R);         a left kernel of R,
    1919RightKernel(module R);        a right kernel of R,
    20 LeftInverse(module R)         a left inverse of matrix (module).
     20LeftInverse(module R)         a left inverse of matrix (module),
     21
     22smith(module M);              a Smith form of a module M.
    2123
    2224AUXILIARY PROCEDURES:
     
    6163LIB "poly.lib";
    6264LIB "primdec.lib";
     65LIB "matrix.lib";
    6366// -- noncommutative procs:
    6467LIB "ncalg.lib";
     68LIB "nctools.lib";
    6569LIB "nchomolog.lib";
    6670LIB "gkdim.lib";
     
    270274};
    271275//------------------------------------------------------------------------
    272 proc LeftInverse(matrix M)
    273 "USAGE:  LeftInverse(M);
    274            M:  matrix
    275 RETURN:  left inverse of M if exists, i.e., matrix L such that LM == id;
     276proc LeftInverse(module M)
     277"USAGE:  LeftInverse(M);  M a module
     278RETURN:  matrix L such that LM == Id;
    276279EXAMPLE:  example LeftInverse; shows an example
     280NOTE: exists only in the case when Id \subseteq M!
    277281"
    278282{
     283  // now it works also for the NC case;
    279284  int NCols=ncols(M);
    280   M=transpose(M);
    281   //  matrix I[NCols][NCols];
    282   //  I=I+1;
    283   //  module Id=I;
    284285  module Id = freemodule(NCols);
    285   return( transpose( lift( module(M),Id ) )  );
     286  M = transpose(M);
     287  Id = std(Id);
     288  matrix T;
     289  option(redSB);
     290  option(redTail);
     291  // check the correctness:
     292  // Id \subseteq M
     293  module N = liftstd(M, T);
     294  ideal TT = ideal(NF(N,Id));
     295  TT = simplify(TT,2);
     296
     297  if (TT[1]!=0)
     298  {
     299    "No left inverse exists";
     300    return(matrix(0));
     301  }
     302  matrix T2 = lift(N, Id);
     303  T2 = transpose(T2)*transpose(T);
     304  T2 =  transpose(T2);
     305  return(T2);
    286306}
    287307example
     
    295315static proc dim_Our(module R)
    296316{
    297   return(GKdim(R));
     317  int d;
     318  if (Is_NC())
     319  {
     320    d = GKdim(R);
     321  }
     322  else { d = dim(R); }
     323  return(d);
    298324}
    299325//-----------------------------------------------------------------------
     
    311337}
    312338//-----------------------------------------------------------------------
    313 static proc Ext_Our(int i, module R,list #)
     339static proc Ext_Our(int i, module R, list #)
    314340// just a copy of 'Ext_R' from "homolog.lib" in commutative case
    315341{
    316   if ( size(#)==0 )
    317   {
    318     return( ncExt_R(i,R) );
     342  int ISNC = Is_NC();
     343  int ExtraArg = ( size(#)>0 );
     344  if (ISNC)
     345  {
     346    // ExtraArg not yet implemented and thus ignored
     347    return( ncExt_R(i,R) );
     348  }
     349  // the commutative case:
     350  if (ExtraArg)
     351  {
     352    return( Ext_R(i,R,#[1]) );
    319353  }
    320354  else
    321355  {
    322     return( Ext_R(i,R,#[1]) );
    323   };
     356    return( Ext_R(i,R) );
     357  }
    324358}
    325359//------------------------------------------------------------------------
     
    360394                  "obstruction to controllability",
    361395                   Ext_1,
    362                   "annihilator of torsion module(of obstruction to controllability)",
     396                  "annihilator of torsion module (of obstruction to controllability)",
    363397                   Ann_Our(Ext_1),
    364398                   DofS,
     
    871905      }
    872906      //      if( (d-d_lead != 0) || (LExp > 1) )
    873       if ( ( (d-d_lead) != 0) || (LExp > 1) || ( (LExp==0) && !( (d_lead==1) || (d_lead==-1)) ) )
    874       {
    875         "i,j";
    876         i,j;
    877         return( "wrong input" );
     907if ( ( (d-d_lead) != 0) || (LExp > 1) || ( (LExp==0) && ((d_lead>1) || (d_lead<-1)) ) )
     908      {
     909        return(theta);
    878910      }
    879911
     
    10311063    }
    10321064  }
    1033   // return the result [in string=format]
     1065  // return the result [in string-format]
    10341066  L = simplify(L,2+4); // skip zeroes and double entries
    10351067  list SL;
     
    10711103    {
    10721104      M[i] = std(M[i]);
    1073       M[i] = prune(M[i]); // mimimal embedding
    1074       M[i] = std(M[i]);
     1105      //      M[i] = prune(M[i]); // mimimal embedding
     1106      //      M[i] = std(M[i]);
    10751107    }
    10761108  }
     
    10901122  view(CC);
    10911123}
     1124
     1125//---------------------------------------------------------------
     1126
     1127static proc smdeg(matrix N) 
     1128// returns an intvec of length 2 with the index of an element of N with smallest degree
     1129{
     1130  int n = nrows(N);
     1131  int m = ncols(N);
     1132  int d,d_temp;
     1133  intvec v;
     1134  int i,j;            // counter
     1135 
     1136  for(i=1;i<=n;i++) // hier wird ein Element ausgewaehlt(!=0) und mit dessen Grad gestartet
     1137  {
     1138    for(j=1;j<=m;j++)
     1139    {
     1140      if( deg(N[i,j])!=-1 )
     1141      {
     1142    d=deg(N[i,j]);
     1143    break;
     1144      }
     1145    }
     1146    if(d!=-1)
     1147    {
     1148      break;
     1149    } 
     1150  }
     1151  for(i =1; i<=n; i++)
     1152  {
     1153    for(j=1;j<=m;j++)
     1154    {
     1155      d_temp = deg(N[i,j]);
     1156      if ((d_temp < d) && (N[i,j]!=0) )
     1157      {
     1158    d=d_temp;
     1159      }
     1160    }
     1161  }
     1162  for (i=1; i<=n; i++)
     1163  {
     1164    for (j =1; j<=m;j++)
     1165    {
     1166      if ( (deg(N[i,j]) == d) && (N[i,j]!=0) )
     1167      {
     1168    v = i,j;
     1169    return(v);
     1170      }
     1171    }
     1172  }
     1173}
     1174
     1175static proc NoNon0Pol(vector v)
     1176// returns 1, if there is only one non-zero element in v and 0 else
     1177{
     1178  int i,j;
     1179  int n = nrows(v);
     1180  for(j=1; j<=n;j++)
     1181  {
     1182    if(v[j]!=0)
     1183    {
     1184      i++;
     1185    }
     1186  }
     1187  if ( i!=1 )
     1188  {
     1189    i=0;
     1190  }
     1191  return(i);
     1192}
     1193
     1194
     1195proc smith( module M )
     1196"USAGE: smith(M), M a module or a matrix,
     1197RETURN: a list of length 4 with the following entries:
     1198[1]: The Smith-Form S of M,
     1199[2]: the rank of M,
     1200[3]: a unimodular matrix U,
     1201[4]: a unimodular matrix V,
     1202such that U*M*V=S. An empty list is returned when no Smith Form exists.
     1203NOTE: The Smith form only exists over PIDs (principal ideal domains).
     1204"
     1205{
     1206  if (nvars(basering)>1)
     1207    {
     1208      list @L;
     1209      return (@L);
     1210    }
     1211  matrix N = matrix(M);         //Typecasting
     1212  int n = nrows(N);
     1213  int m = ncols(N);
     1214  matrix P = unitmat(n);       //left transformation matrix
     1215  matrix Q = unitmat(m);       //right transformation matrix
     1216  int k, i, j, deg_temp;
     1217  poly tmp;
     1218  vector v;
     1219  list L;                      //for extgcd-computation
     1220  intmat f[1][n];              //to save degrees
     1221  matrix lambda[1][n];         //to save leadcoefficients
     1222  intmat g[1][m];              //to save degrees
     1223  matrix mu[1][m];             //to save leadcoefficients
     1224  int ii;                       //counter
     1225  ideal J;                      //for extgcd-computations
     1226  matrix T;                     //----------"------------
     1227 
     1228  while ((k!=n) && (k!=m) )   
     1229  {
     1230    k++;
     1231    while ((k<=n) && (k<=m))  //outer while-loop for column-operations
     1232    {
     1233      while(k<=m )        //inner while-loop for row-operations
     1234      {
     1235    if( (n>m) && (k < n) && (k<m))
     1236    {     
     1237      if( simplify((ideal(submat(N,k+1..n,k+1..m))),2)== 0)
     1238      {
     1239        return(N,k-1,P,Q);
     1240      }
     1241    }
     1242        i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix
     1243        i=i+(k-1);                        //indices adjusted to the whole matrix
     1244        j=j+(k-1);
     1245        if(i!=k)                    //take the element with smallest degree in the first position
     1246        {
     1247          N=permrow(N,i,k);
     1248          P=permrow(P,i,k);
     1249        }
     1250        if(j!=k)
     1251        {
     1252          N=permcol(N,j,k);
     1253          Q=permcol(Q,j,k);
     1254        }
     1255        if(NoNon0Pol(N[k])==1)
     1256        {
     1257          break;
     1258        }
     1259        tmp=leadcoef(N[k,k]);
     1260        deg_temp=ord(N[k,k]);             //ord outputs the leading degree of N[k,k]
     1261        for(ii=k+1;ii<=n;ii++)
     1262        {
     1263          lambda[1,ii]=leadcoef(N[ii,k])/tmp;     
     1264          f[1,ii]=deg(N[ii,k])-deg_temp;
     1265        }
     1266        for(ii=k+1;ii<=n;ii++)
     1267        {
     1268          N = addrow(N,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
     1269      N = normalize(N);
     1270          P = addrow(P,k,-lambda[1,ii]*var(1)^f[1,ii],ii);
     1271    }
     1272      }
     1273      if (k>n)
     1274      {
     1275    break;
     1276      }
     1277      if(NoNon0Pol(transpose(N)[k])==1)
     1278      {
     1279        break;
     1280      }
     1281      tmp=leadcoef(N[k,k]);
     1282      deg_temp=ord(N[k,k]); //ord outputs the leading degree of N[k][k]
     1283     
     1284      for(ii=k+1;ii<=m;ii++)
     1285      {
     1286        mu[1,ii]=leadcoef(N[k,ii])/tmp;     
     1287        g[1,ii]=deg(N[k,ii])-deg_temp;
     1288      }
     1289      for(ii=k+1;ii<=m;ii++)
     1290      {
     1291        N=addcol(N,k,-mu[1,ii]*var(1)^g[1,ii],ii);
     1292        N=normalize(N);
     1293        Q=addcol(Q,k,-mu[1,ii]*var(1)^g[1,ii],ii);
     1294      }
     1295    }
     1296    if( (k!=1) && (k<n) && (k<m) )
     1297    {
     1298      L=extgcd(N[k-1,k-1],N[k,k]);  //one can use this line if extgcd-bug is fixed
     1299     
     1300      //the next lines are needed for the divisibility condition of the entries of the Smith Form
     1301     
     1302      if(!(N[k-1,k-1]==L[1]))  //means that N[k-1,k-1] is not a divisor of N[k,k]
     1303        {
     1304          N=addrow(N,k-1,L[2],k);
     1305          normalize(N);
     1306          P=addrow(P,k-1,L[2],k);
     1307         
     1308          N=addcol(N,k,-L[3],k-1);
     1309          N=normalize(N);
     1310          Q=addcol(Q,k,-L[3],k-1);
     1311          k=k-2;
     1312        }
     1313    }
     1314  } 
     1315  if( (k<=n) && (k<=m) )
     1316    {
     1317      if( N[k,k]==0)
     1318      {
     1319        return(N,k-1,P,Q);
     1320      }
     1321    }
     1322  return(N,k,P,Q);
     1323}
     1324example
     1325{
     1326  "EXAMPLE:";echo = 2;
     1327  option(redSB);
     1328  option(redTail);
     1329  ring r=0,x,dp;
     1330  // This is example is trivial, but we want to see,
     1331  // that nothing is changed, when the matrix is already in Smith-Form!
     1332  module M = [x,0,0],[0,x2,0],[0,0,x3];
     1333  print(M);
     1334  list L = smith(M);
     1335  print(L[1]);
     1336  matrix N=matrix(M);
     1337  matrix B=L[3]*N*L[4];
     1338  B=normalize(std(B));
     1339  print(B);
     1340  //------- and yet another example --------------
     1341  module M2=[x2,x,3x3-4],[2x2-1,4x,5x2],[2x5,3x,4x];
     1342  print(M2);
     1343  list P=smith(M2);
     1344  print(P[1]);
     1345  matrix N2=matrix(M2);
     1346  matrix B2=P[3]*N2*P[4];
     1347  print(B2);
     1348  B2=normalize(std(B2));
     1349  print(B2);
     1350}
Note: See TracChangeset for help on using the changeset viewer.