Changeset a7f2a04 in git
 Timestamp:
 Feb 16, 2002, 11:56:19 AM (21 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
 Children:
 03495f1856f489d63ec19f337bd5240315ef5371
 Parents:
 502ea23f40b3a43fb8b3bcbc306d20b0c02ab50b
 Location:
 Singular
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

Singular/eigenval.cc
r502ea2 ra7f2a04 1 /**************************************** 2 * Computer Algebra System SINGULAR *3 **************************************** /4 /* $Id: eigenval.cc,v 1. 5 20020119 14:48:14 obachmanExp $ */1 /***************************************** 2 * Computer Algebra System SINGULAR * 3 *****************************************/ 4 /* $Id: eigenval.cc,v 1.6 20020216 10:56:14 mschulze Exp $ */ 5 5 /* 6 6 * ABSTRACT: eigenvalues of constant square matrices … … 11 11 #ifdef HAVE_EIGENVAL 12 12 13 #include "febase.h" 13 14 #include "tok.h" 14 #include "febase.h" 15 #include "ipid.h" 16 #include "intvec.h" 15 17 #include "numbers.h" 16 18 #include "polys.h" 17 19 #include "ideals.h" 20 #include "lists.h" 18 21 #include "matpol.h" 19 22 #include "clapsing.h" 20 21 matrix swap(matrix M,int i,int j) 23 #include "eigenval.h" 24 25 26 matrix evSwap(matrix M,int i,int j) 22 27 { 23 28 if(i==j) 24 return M;25 poly p; 29 return(M); 30 26 31 for(int k=1;k<=MATROWS(M);k++) 27 32 { 28 p =MATELEM(M,i,k);33 poly p=MATELEM(M,i,k); 29 34 MATELEM(M,i,k)=MATELEM(M,j,k); 30 35 MATELEM(M,j,k)=p; 31 36 } 37 32 38 for(int k=1;k<=MATCOLS(M);k++) 33 39 { 34 p =MATELEM(M,k,i);40 poly p=MATELEM(M,k,i); 35 41 MATELEM(M,k,i)=MATELEM(M,k,j); 36 42 MATELEM(M,k,j)=p; 37 43 } 44 38 45 return(M); 39 46 } 40 47 41 matrix rowelim(matrix M,int i, int j, int k) 48 49 BOOLEAN evSwap(leftv res,leftv h) 50 { 51 if(currRingHdl) 52 { 53 if(h&&h>Typ()==MATRIX_CMD) 54 { 55 matrix M=(matrix)h>Data(); 56 h=h>next; 57 if(h&&h>Typ()==INT_CMD) 58 { 59 int i=(int)h>Data(); 60 h=h>next; 61 if(h&&h>Typ()==INT_CMD) 62 { 63 int j=(int)h>Data(); 64 res>rtyp=MATRIX_CMD; 65 res>data=(void*)evSwap(mpCopy(M),i,j); 66 return FALSE; 67 } 68 } 69 } 70 WerrorS("<matrix>,<int>,<int> expected"); 71 return TRUE; 72 } 73 WerrorS("no ring active"); 74 return TRUE; 75 } 76 77 78 matrix evRowElim(matrix M,int i,int j,int k) 42 79 { 43 80 if(MATELEM(M,i,k)==NULLMATELEM(M,j,k)==NULL) 44 return M; 45 number n=nDiv(pGetCoeff(MATELEM(M,i,k)), 46 pGetCoeff(MATELEM(M,j,k))); 81 return(M); 82 83 poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k)))); 84 47 85 for(int l=1;l<=MATCOLS(M);l++) 48 { 49 MATELEM(M,i,l)=pSub(MATELEM(M,i,l),pMult_nn(pCopy(MATELEM(M,j,l)),n)); 50 } 86 MATELEM(M,i,l)=pSub(MATELEM(M,i,l),pMult(pCopy(p),pCopy(MATELEM(M,j,l)))); 87 51 88 for(int l=1;l<=MATROWS(M);l++) 52 { 53 MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),pMult_nn(pCopy(MATELEM(M,l,i)),n)); 54 } 55 nDelete(&n); 56 return M; 57 } 58 59 matrix colelim(matrix M,int i, int j, int k) 60 { 61 if(MATELEM(M,k,i)==NULLMATELEM(M,k,j)==NULL) 62 return M; 63 number n=nDiv(pGetCoeff(MATELEM(M,k,i)), 64 pGetCoeff(MATELEM(M,k,j))); 89 MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),pMult(pCopy(p),pCopy(MATELEM(M,l,i)))); 90 91 pDelete(&p); 92 93 return(M); 94 } 95 96 97 BOOLEAN evRowElim(leftv res,leftv h) 98 { 99 if(currRingHdl) 100 { 101 if(h&&h>Typ()==MATRIX_CMD) 102 { 103 matrix M=(matrix)h>Data(); 104 h=h>next; 105 if(h&&h>Typ()==INT_CMD) 106 { 107 int i=(int)h>Data(); 108 h=h>next; 109 if(h&&h>Typ()==INT_CMD) 110 { 111 int j=(int)h>Data(); 112 h=h>next; 113 if(h&&h>Typ()==INT_CMD) 114 { 115 int k=(int)h>Data(); 116 res>rtyp=MATRIX_CMD; 117 res>data=(void*)evRowElim(mpCopy(M),i,j,k); 118 return FALSE; 119 } 120 } 121 } 122 } 123 WerrorS("<matrix>,<int>,<int>,<int> expected"); 124 return TRUE; 125 } 126 WerrorS("no ring active"); 127 return TRUE; 128 } 129 130 131 matrix evColElim(matrix M,int i,int j,int k) 132 { 133 if(MATELEM(M,k,i)==0MATELEM(M,k,j)==0) 134 return(M); 135 136 poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j)))); 137 65 138 for(int l=1;l<=MATROWS(M);l++) 66 { 67 MATELEM(M,l,i)=pSub(MATELEM(M,l,i),pMult_nn(pCopy(MATELEM(M,l,j)),n)); 68 } 139 MATELEM(M,l,i)=pSub(MATELEM(M,l,i),pMult(pCopy(p),pCopy(MATELEM(M,l,j)))); 140 69 141 for(int l=1;l<=MATCOLS(M);l++) 70 { 71 MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),pMult_nn(pCopy(MATELEM(M,i,l)),n)); 72 } 73 nDelete(&n); 74 return M; 75 } 76 77 matrix tridiag(matrix M) 78 { 79 int n=MATCOLS(M); 80 for(int k=1;k<=n2;k++) 81 { 82 int j=k+1; 83 while(j<=n&&MATELEM(M,j,k)==NULL) j++; 84 if(j<=n) 85 { 86 for(int i=j+1;i<=n;i++) 87 { 88 M=rowelim(M,i,j,k); 89 } 90 M=swap(M,j,k+1); 91 } 142 MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),pMult(pCopy(p),pCopy(MATELEM(M,i,l)))); 143 144 pDelete(&p); 145 146 return(M); 147 } 148 149 150 BOOLEAN evColElim(leftv res,leftv h) 151 { 152 if(currRingHdl) 153 { 154 if(h&&h>Typ()==MATRIX_CMD) 155 { 156 matrix M=(matrix)h>Data(); 157 h=h>next; 158 if(h&&h>Typ()==INT_CMD) 159 { 160 int i=(int)h>Data(); 161 h=h>next; 162 if(h&&h>Typ()==INT_CMD) 163 { 164 int j=(int)h>Data(); 165 h=h>next; 166 if(h&&h>Typ()==INT_CMD) 167 { 168 int k=(int)h>Data(); 169 res>rtyp=MATRIX_CMD; 170 res>data=(void*)evColElim(mpCopy(M),i,j,k); 171 return FALSE; 172 } 173 } 174 } 175 } 176 WerrorS("<matrix>,<int>,<int>,<int> expected"); 177 return TRUE; 178 } 179 WerrorS("no ring active"); 180 return TRUE; 181 } 182 183 184 matrix evHessenberg(matrix M) 185 { 186 int n=MATROWS(M); 187 int i,j; 188 189 for(int k=1;k<n1;k++) 190 { 92 191 j=k+1; 93 while(j<=n&&MATELEM(M,k,j)==NULL) j++; 94 if(j<=n) 95 { 96 for(int i=j+1;i<=n;i++) 97 { 98 M=colelim(M,i,j,k); 99 } 100 M=swap(M,j,k+1); 101 } 102 } 192 while(j<n&&MATELEM(M,j,k)==0) 193 j++; 194 195 if(MATELEM(M,j,k)!=0) 196 { 197 M=evSwap(M,j,k+1); 198 199 for(i=j+1;i<=n;i++) 200 M=evRowElim(M,i,k+1,k); 201 } 202 } 203 103 204 return(M); 104 205 } 105 206 106 lists addval(lists l,poly e0,int m0) 107 { 108 ideal ee=(ideal)l>m[0].data; 109 intvec *mm=(intvec*)l>m[1].data; 110 int n=0; 111 if(ee!=NULL) 112 n=IDELEMS(ee); 113 for(int i=n1;i>=0;i) 114 { 115 if(pEqualPolys(ee>m[i],e0)) 116 { 117 (*mm)[i]+=m0; 118 return l; 119 } 120 } 121 ideal e=idInit(n+1,1); 122 for(int i=n1;i>=0;i) 123 { 124 e>m[i]=ee>m[i]; 125 ee>m[i]=NULL; 126 } 127 e>m[n]=e0; 128 l>m[0].data=e; 129 if(ee!=NULL) 130 idDelete(&ee); 131 mm>resize(n+1); 132 (*mm)[n]=m0; 133 return l; 134 } 135 136 lists sortval(lists l) 137 { 138 ideal ee=(ideal)l>m[0].data; 139 intvec *mm=(intvec*)l>m[1].data; 140 int n=IDELEMS(ee); 141 for(int i=n1;i>=1;i) 142 for(int j=i1;j>=0;j) 143 if(nGreater(pGetCoeff(ee>m[j]),pGetCoeff(ee>m[i]))) 144 { 145 poly e=ee>m[i]; 146 ee>m[i]=ee>m[j]; 147 ee>m[j]=e; 148 int m=(*mm)[i]; 149 (*mm)[i]=(*mm)[j]; 150 (*mm)[j]=m; 151 } 152 return l; 153 } 154 155 lists eigenval(matrix M) 156 { 157 M=tridiag(M); 158 int n=MATCOLS(M); 207 208 BOOLEAN evHessenberg(leftv res,leftv h) 209 { 210 if(currRingHdl) 211 { 212 if(h&&h>Typ()==MATRIX_CMD) 213 { 214 matrix M=(matrix)h>Data(); 215 res>rtyp=MATRIX_CMD; 216 res>data=(void*)evHessenberg(mpCopy(M)); 217 return FALSE; 218 } 219 WerrorS("<matrix> expected"); 220 return TRUE; 221 } 222 WerrorS("no ring active"); 223 return TRUE; 224 } 225 226 227 lists evEigenvalue(matrix M) 228 { 159 229 lists l=(lists)omAllocBin(slists_bin); 230 if(MATROWS(M)!=MATCOLS(M)) 231 { 232 l>Init(0); 233 return(l); 234 } 235 236 M=evHessenberg(M); 237 238 int n=MATROWS(M); 239 ideal e=idInit(n,1); 240 intvec *m=new intvec(n); 241 242 poly t=pOne(); 243 pSetExp(t,1,1); 244 245 for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++) 246 { 247 while(j<=n&&MATELEM(M,j,j1)!=NULL) 248 j++; 249 if(j==j0+1) 250 { 251 e>m[k]=pHead(MATELEM(M,j0,j0)); 252 (*m)[k]=1; 253 k++; 254 } 255 else 256 { 257 int n0=jj0; 258 matrix M0=mpNew(n0,n0); 259 260 j0; 261 for(int i=1;i<=n0;i++) 262 for(int j=1;j<=n0;j++) 263 MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j)); 264 for(int i=1;i<=n0;i++) 265 MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t)); 266 267 intvec *m0; 268 ideal e0=singclap_factorize(mpDetBareiss(M0),&m0,2); 269 270 for(int i=0;i<IDELEMS(e0);i++) 271 { 272 number e1=nNeg(pGetCoeff(e0>m[i])); 273 pDeleteLm(&e0>m[i]); 274 if(pGetExp(e0>m[i],1)==0) 275 e>m[k]=pNSet(nDiv(pGetCoeff(e0>m[i]),e1)); 276 else 277 e>m[k]=pNSet(nDiv(e1,pGetCoeff(e0>m[i]))); 278 nDelete(&e1); 279 (*m)[k]=(*m0)[i]; 280 k++; 281 } 282 283 delete(m0); 284 idDelete(&e0); 285 } 286 } 287 288 pDelete(&t); 289 idDelete((ideal *)&M); 290 291 for(int i=0;i<n1;i++) 292 { 293 if(e>m[i]!=NULL) 294 for(int j=i+1;j<n;j++) 295 { 296 if(e>m[j]!=NULL) 297 if(nEqual(pGetCoeff(e>m[i]),pGetCoeff(e>m[j]))) 298 { 299 (*m)[i]+=(*m)[j]; 300 (*m)[j]=0; 301 } 302 else 303 if(nGreater(pGetCoeff(e>m[i]),pGetCoeff(e>m[j]))) 304 { 305 poly p=e>m[i]; 306 e>m[i]=e>m[j]; 307 e>m[j]=p; 308 int k=(*m)[i]; 309 (*m)[i]=(*m)[j]; 310 (*m)[j]=k; 311 } 312 } 313 } 314 315 int n0=0; 316 for(int i=0;i<n;i++) 317 if((*m)[i]>0) 318 n0++; 319 320 ideal e0=idInit(n0,1); 321 intvec *m0=new intvec(n0); 322 323 for(int i=0,i0=0;i<n;i++) 324 if((*m)[i]>0) 325 { 326 e0>m[i0]=e>m[i]; 327 e>m[i]=NULL; 328 (*m0)[i0]=(*m)[i]; 329 i0++; 330 } 331 332 idDelete(&e); 333 delete(m); 334 160 335 l>Init(2); 161 336 l>m[0].rtyp=IDEAL_CMD; 162 l>m[0].data= NULL;337 l>m[0].data=e0; 163 338 l>m[1].rtyp=INTVEC_CMD; 164 l>m[1].data=new intvec; 165 int j=1; 166 while(j<=n) 167 { 168 while(j<n&&(MATELEM(M,j,j+1)==NULLMATELEM(M,j+1,j)==NULL)j==n) 169 { 170 l=addval(l,MATELEM(M,j,j),1); 171 MATELEM(M,j,j)=NULL; 172 j++; 173 } 174 if(j<n) 175 { 176 poly t=pOne(); 177 pSetExp(t,1,1); 178 pSetm(t); 179 poly d0=pSub(MATELEM(M,j,j),t); 180 MATELEM(M,j,j)=NULL; 181 j++; 182 poly d1=pOne(); 183 poly d2=NULL; 184 while(j<=n&&MATELEM(M,j,j1)!=NULL&&MATELEM(M,j1,j)!=NULL) 185 { 186 d2=d1; 187 d1=pCopy(d0); 188 poly t=pOne(); 189 pSetExp(t,1,1); 190 pSetm(t); 191 d0=pSub(pMult(d0,pSub(MATELEM(M,j,j),t)), 192 pMult(pMult(d2,MATELEM(M,j1,j)),MATELEM(M,j,j1))); 193 MATELEM(M,j,j)=NULL; 194 MATELEM(M,j1,j)=NULL; 195 MATELEM(M,j,j1)=NULL; 196 j++; 197 } 198 pDelete(&d1); 199 intvec *m0=NULL; 200 #ifdef HAVE_FACTORY 201 ideal e0=singclap_factorize(d0,&m0,0); 202 #else 203 ideal e0 = NULL; 204 #endif 205 pDelete(&d0); 206 for(int i=IDELEMS(e0)1;i>=1;i) 207 { 208 poly p=e0>m[i]; 209 e0>m[i]=NULL; 210 poly p0,p1; 211 poly pp=p; 212 while(pp!=NULL&&pGetExp(pp,1)<=1) 213 { 214 if(pGetExp(pp,1)==0) 215 p0=pp; 216 else 217 if(pGetExp(pp,1)==1) 218 p1=pp; 219 pp=pNext(pp); 220 } 221 if(pp==NULL) 222 { 223 pp=p; 224 p=pNSet(nNeg(nDiv(pGetCoeff(p0),pGetCoeff(p1)))); 225 pDelete(&pp); 226 } 227 else 228 { 229 p=pMult_nn(p,pGetCoeff(e0>m[0])); 230 } 231 l=addval(l,p,(*m0)[i]); 232 } 233 delete m0; 234 idDelete(&e0); 235 } 236 } 237 idDelete((ideal*)&M); 238 return sortval(l); 239 } 240 241 BOOLEAN tridiag(leftv res,leftv h) 242 { 243 if((h!=NULL) && (h>Typ()==MATRIX_CMD)) 244 { 245 matrix M=(matrix)h>Data(); 246 if(MATCOLS(M)!=MATROWS(M)) 247 { 248 WerrorS("square matrix expected"); 249 } 250 res>rtyp=MATRIX_CMD; 251 res>data=(void*)tridiag(mpCopy(M)); 252 return FALSE; 253 } 254 WerrorS("<matrix> expected"); 255 return TRUE; 256 } 257 258 BOOLEAN eigenval(leftv res,leftv h) 259 { 260 if((h!=NULL) && (h>Typ()==MATRIX_CMD)) 261 { 262 matrix M=(matrix)h>Data(); 263 if(MATCOLS(M)!=MATROWS(M)) 264 { 265 WerrorS("square matrix expected"); 266 } 267 res>rtyp=LIST_CMD; 268 res>data=(void*)eigenval(mpCopy(M)); 269 return FALSE; 270 } 271 WerrorS("<matrix> expected"); 272 return TRUE; 273 } 274 #endif 339 l>m[1].data=m0; 340 341 return(l); 342 } 343 344 345 BOOLEAN evEigenvalue(leftv res,leftv h) 346 { 347 if(currRingHdl) 348 { 349 if(h&&h>Typ()==MATRIX_CMD) 350 { 351 matrix M=(matrix)h>Data(); 352 res>rtyp=LIST_CMD; 353 res>data=(void*)evEigenvalue(mpCopy(M)); 354 return FALSE; 355 } 356 WerrorS("<matrix> expected"); 357 return TRUE; 358 } 359 WerrorS("no ring active"); 360 return TRUE; 361 } 362 363 #endif /* HAVE_EIGENVAL */ 
Singular/eigenval.h
r502ea2 ra7f2a04 2 2 * Computer Algebra System SINGULAR * 3 3 ****************************************/ 4 /* $Id: eigenval.h,v 1. 1 20010305 17:31:40mschulze Exp $ */4 /* $Id: eigenval.h,v 1.2 20020216 10:56:19 mschulze Exp $ */ 5 5 /* 6 6 * ABSTRACT: eigenvalues of constant square matrices … … 10 10 #define EIGENVAL_H 11 11 12 matrix tridiag(matrix M); 13 lists eigenval(matrix M); 14 BOOLEAN tridiag(leftv res,leftv h); 15 BOOLEAN eigenval(leftv res,leftv h); 12 matrix evSwap(matrix M,int i,int j); 13 BOOLEAN evSwap(leftv res,leftv h); 14 matrix evRowElim(matrix M,int i,int j,int k); 15 BOOLEAN evRowElim(leftv res,leftv h); 16 matrix evColElim(matrix M,int i,int j,int k); 17 BOOLEAN evColElim(leftv res,leftv h); 18 matrix evHessenberg(matrix M); 19 BOOLEAN evHessenberg(leftv res,leftv h); 20 lists evEigenvalue(matrix M); 21 BOOLEAN evEigenvalue(leftv res,leftv h); 16 22 17 #endif 23 #endif /* EIGENVAL_H */
Note: See TracChangeset
for help on using the changeset viewer.