Changeset 501f7c5 in git
 Timestamp:
 Dec 16, 2004, 5:28:53 PM (20 years ago)
 Branches:
 (u'fiekerDuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
 Children:
 175c561dc37255410efb4b041fbb41a6e1b0d8b8
 Parents:
 1d5418b9bdcc715becec9d153116793ed3147e8f
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/control.lib
r1d5418 r501f7c5 1 version="$Id: control.lib,v 1.1 7 20041214 21:05:09levandov Exp $";1 version="$Id: control.lib,v 1.18 20041216 16:28:53 levandov Exp $"; 2 2 category="Applications"; 3 3 info=" … … 18 18 LeftKernel(module R); a left kernel of R, 19 19 RightKernel(module R); a right kernel of R, 20 LeftInverse(module R) a left inverse of matrix (module). 20 LeftInverse(module R) a left inverse of matrix (module), 21 22 smith(module M); a Smith form of a module M. 21 23 22 24 AUXILIARY PROCEDURES: … … 61 63 LIB "poly.lib"; 62 64 LIB "primdec.lib"; 65 LIB "matrix.lib"; 63 66 //  noncommutative procs: 64 67 LIB "ncalg.lib"; 68 LIB "nctools.lib"; 65 69 LIB "nchomolog.lib"; 66 70 LIB "gkdim.lib"; … … 270 274 }; 271 275 // 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; 276 proc LeftInverse(module M) 277 "USAGE: LeftInverse(M); M a module 278 RETURN: matrix L such that LM == Id; 276 279 EXAMPLE: example LeftInverse; shows an example 280 NOTE: exists only in the case when Id \subseteq M! 277 281 " 278 282 { 283 // now it works also for the NC case; 279 284 int NCols=ncols(M); 280 M=transpose(M);281 // matrix I[NCols][NCols];282 // I=I+1;283 // module Id=I;284 285 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); 286 306 } 287 307 example … … 295 315 static proc dim_Our(module R) 296 316 { 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); 298 324 } 299 325 // … … 311 337 } 312 338 // 313 static proc Ext_Our(int i, module R, list #)339 static proc Ext_Our(int i, module R, list #) 314 340 // just a copy of 'Ext_R' from "homolog.lib" in commutative case 315 341 { 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]) ); 319 353 } 320 354 else 321 355 { 322 return( Ext_R(i,R ,#[1]) );323 } ;356 return( Ext_R(i,R) ); 357 } 324 358 } 325 359 // … … 360 394 "obstruction to controllability", 361 395 Ext_1, 362 "annihilator of torsion module (of obstruction to controllability)",396 "annihilator of torsion module (of obstruction to controllability)", 363 397 Ann_Our(Ext_1), 364 398 DofS, … … 871 905 } 872 906 // if( (dd_lead != 0)  (LExp > 1) ) 873 if ( ( (dd_lead) != 0)  (LExp > 1)  ( (LExp==0) && !( (d_lead==1)  (d_lead==1)) ) ) 874 { 875 "i,j"; 876 i,j; 877 return( "wrong input" ); 907 if ( ( (dd_lead) != 0)  (LExp > 1)  ( (LExp==0) && ((d_lead>1)  (d_lead<1)) ) ) 908 { 909 return(theta); 878 910 } 879 911 … … 1031 1063 } 1032 1064 } 1033 // return the result [in string =format]1065 // return the result [in stringformat] 1034 1066 L = simplify(L,2+4); // skip zeroes and double entries 1035 1067 list SL; … … 1071 1103 { 1072 1104 M[i] = std(M[i]); 1073 M[i] = prune(M[i]); // mimimal embedding1074 M[i] = std(M[i]);1105 // M[i] = prune(M[i]); // mimimal embedding 1106 // M[i] = std(M[i]); 1075 1107 } 1076 1108 } … … 1090 1122 view(CC); 1091 1123 } 1124 1125 // 1126 1127 static 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 1175 static proc NoNon0Pol(vector v) 1176 // returns 1, if there is only one nonzero 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 1195 proc smith( module M ) 1196 "USAGE: smith(M), M a module or a matrix, 1197 RETURN: a list of length 4 with the following entries: 1198 [1]: The SmithForm S of M, 1199 [2]: the rank of M, 1200 [3]: a unimodular matrix U, 1201 [4]: a unimodular matrix V, 1202 such that U*M*V=S. An empty list is returned when no Smith Form exists. 1203 NOTE: 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 extgcdcomputation 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 extgcdcomputations 1226 matrix T; //" 1227 1228 while ((k!=n) && (k!=m) ) 1229 { 1230 k++; 1231 while ((k<=n) && (k<=m)) //outer whileloop for columnoperations 1232 { 1233 while(k<=m ) //inner whileloop for rowoperations 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,k1,P,Q); 1240 } 1241 } 1242 i,j = smdeg(submat(N,k..n,k..m)); //choose smallest degree in the remaining submatrix 1243 i=i+(k1); //indices adjusted to the whole matrix 1244 j=j+(k1); 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[k1,k1],N[k,k]); //one can use this line if extgcdbug is fixed 1299 1300 //the next lines are needed for the divisibility condition of the entries of the Smith Form 1301 1302 if(!(N[k1,k1]==L[1])) //means that N[k1,k1] is not a divisor of N[k,k] 1303 { 1304 N=addrow(N,k1,L[2],k); 1305 normalize(N); 1306 P=addrow(P,k1,L[2],k); 1307 1308 N=addcol(N,k,L[3],k1); 1309 N=normalize(N); 1310 Q=addcol(Q,k,L[3],k1); 1311 k=k2; 1312 } 1313 } 1314 } 1315 if( (k<=n) && (k<=m) ) 1316 { 1317 if( N[k,k]==0) 1318 { 1319 return(N,k1,P,Q); 1320 } 1321 } 1322 return(N,k,P,Q); 1323 } 1324 example 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 SmithForm! 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,3x34],[2x21,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.