My Project
Loading...
Searching...
No Matches
eigenval_ip.cc File Reference
#include "kernel/mod2.h"
#include "Singular/tok.h"
#include "Singular/ipid.h"
#include "misc/intvec.h"
#include "coeffs/numbers.h"
#include "kernel/polys.h"
#include "kernel/ideals.h"
#include "Singular/lists.h"
#include "polys/matpol.h"
#include "polys/clapsing.h"
#include "kernel/linear_algebra/eigenval.h"
#include "Singular/ipshell.h"
#include "Singular/eigenval_ip.h"

Go to the source code of this file.

Functions

BOOLEAN evSwap (leftv res, leftv h)
 
BOOLEAN evRowElim (leftv res, leftv h)
 
BOOLEAN evHessenberg (leftv res, leftv h)
 
lists evEigenvals (matrix M)
 
BOOLEAN evEigenvals (leftv res, leftv h)
 

Function Documentation

◆ evEigenvals() [1/2]

BOOLEAN evEigenvals ( leftv res,
leftv h )

Definition at line 269 of file eigenval_ip.cc.

270{
271 if(currRing)
272 {
273 if(h&&h->Typ()==MATRIX_CMD)
274 {
275 matrix M=(matrix)h->CopyD();
276 res->rtyp=LIST_CMD;
277 res->data=(void *)evEigenvals(M);
278 return FALSE;
279 }
280 WerrorS("<matrix> expected");
281 return TRUE;
282 }
283 WerrorS("no ring active");
284 return TRUE;
285}
#define TRUE
Definition auxiliary.h:101
#define FALSE
Definition auxiliary.h:97
lists evEigenvals(matrix M)
CanonicalForm res
Definition facAbsFact.cc:60
void WerrorS(const char *s)
Definition feFopen.cc:24
@ MATRIX_CMD
Definition grammar.cc:287
STATIC_VAR Poly * h
Definition janet.cc:971
ip_smatrix * matrix
Definition matpol.h:43
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
#define M
Definition sirandom.c:25
@ LIST_CMD
Definition tok.h:118

◆ evEigenvals() [2/2]

lists evEigenvals ( matrix M)

Definition at line 120 of file eigenval_ip.cc.

121{
123 if(MATROWS(M)!=MATCOLS(M))
124 {
125 l->Init(0);
126 return(l);
127 }
128
130
131 int n=MATCOLS(M);
132 ideal e=idInit(n,1);
133 intvec *m=new intvec(n);
134
135 poly t=pOne();
136 pSetExp(t,1,1);
137 pSetm(t);
138
139 for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
140 {
141 while(j<=n&&MATELEM(M,j,j-1)!=NULL)
142 j++;
143 if(j==j0+1)
144 {
145 e->m[k]=pHead(MATELEM(M,j0,j0));
146 (*m)[k]=1;
147 k++;
148 }
149 else
150 {
151 int n0=j-j0;
152 matrix M0=mpNew(n0,n0);
153
154 j0--;
155 for(int i=1;i<=n0;i++)
156 for(int j=1;j<=n0;j++)
157 MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
158 for(int i=1;i<=n0;i++)
159 MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
160
161 intvec *m0;
163 if (e0==NULL)
164 {
165 l->Init(0);
166 return(l);
167 }
168
169 for(int i=0;i<IDELEMS(e0);i++)
170 {
171 if(pNext(e0->m[i])==NULL)
172 {
173 (*m)[k]=(*m0)[i];
174 k++;
175 }
176 else
177 if(pGetExp(e0->m[i],1)<2&&pGetExp(pNext(e0->m[i]),1)<2&&
178 pNext(pNext(e0->m[i]))==NULL)
179 {
180 number e1=nCopy(pGetCoeff(e0->m[i]));
181 e1=nInpNeg(e1);
182 if(pGetExp(pNext(e0->m[i]),1)==0)
183 e->m[k]=pNSet(nDiv(pGetCoeff(pNext(e0->m[i])),e1));
184 else
185 e->m[k]=pNSet(nDiv(e1,pGetCoeff(pNext(e0->m[i]))));
186 nDelete(&e1);
187 pNormalize(e->m[k]);
188 (*m)[k]=(*m0)[i];
189 k++;
190 }
191 else
192 {
193 e->m[k]=e0->m[i];
194 pNormalize(e->m[k]);
195 e0->m[i]=NULL;
196 (*m)[k]=(*m0)[i];
197 k++;
198 }
199 }
200
201 delete(m0);
202 idDelete(&e0);
203 }
204 }
205
206 pDelete(&t);
207 idDelete((ideal *)&M);
208
209 for(int i0=0;i0<n-1;i0++)
210 {
211 for(int i1=i0+1;i1<n;i1++)
212 {
213 if(pEqualPolys(e->m[i0],e->m[i1]))
214 {
215 (*m)[i0]+=(*m)[i1];
216 (*m)[i1]=0;
217 }
218 else
219 {
220 if((e->m[i0]==NULL&&!nGreaterZero(pGetCoeff(e->m[i1])))||
221 (e->m[i1]==NULL&&
222 (nGreaterZero(pGetCoeff(e->m[i0]))||pNext(e->m[i0])!=NULL))||
223 e->m[i0]!=NULL&&e->m[i1]!=NULL&&
224 ((pNext(e->m[i0])!=NULL&&pNext(e->m[i1])==NULL)||
225 (pNext(e->m[i0])==NULL&&pNext(e->m[i1])==NULL&&
226 nGreater(pGetCoeff(e->m[i0]),pGetCoeff(e->m[i1])))))
227 {
228 poly e1=e->m[i0];
229 e->m[i0]=e->m[i1];
230 e->m[i1]=e1;
231 int m1=(*m)[i0];
232 (*m)[i0]=(*m)[i1];
233 (*m)[i1]=m1;
234 }
235 }
236 }
237 }
238
239 int n0=0;
240 for(int i=0;i<n;i++)
241 if((*m)[i]>0)
242 n0++;
243
244 ideal e0=idInit(n0,1);
245 intvec *m0=new intvec(n0);
246
247 for(int i=0,i0=0;i<n;i++)
248 if((*m)[i]>0)
249 {
250 e0->m[i0]=e->m[i];
251 e->m[i]=NULL;
252 (*m0)[i0]=(*m)[i];
253 i0++;
254 }
255
256 idDelete(&e);
257 delete(m);
258
259 l->Init(2);
260 l->m[0].rtyp=IDEAL_CMD;
261 l->m[0].data=e0;
262 l->m[1].rtyp=INTVEC_CMD;
263 l->m[1].data=m0;
264
265 return(l);
266}
int l
Definition cfEzgcd.cc:100
int m
Definition cfEzgcd.cc:128
int i
Definition cfEzgcd.cc:132
int k
Definition cfEzgcd.cc:99
ideal singclap_factorize(poly f, intvec **v, int with_exps, const ring r)
Definition clapsing.cc:948
BOOLEAN evHessenberg(leftv res, leftv h)
int j
Definition facHensel.cc:110
@ IDEAL_CMD
Definition grammar.cc:285
#define idDelete(H)
delete an ideal
Definition ideals.h:29
VAR omBin slists_bin
Definition lists.cc:23
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition matpol.cc:37
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition matpol.cc:1670
#define MATELEM(mat, i, j)
1-based access to matrix
Definition matpol.h:29
#define MATROWS(i)
Definition matpol.h:26
#define MATCOLS(i)
Definition matpol.h:27
#define pNext(p)
Definition monomials.h:36
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition monomials.h:44
slists * lists
#define nDiv(a, b)
Definition numbers.h:32
#define nDelete(n)
Definition numbers.h:16
#define nInpNeg(n)
Definition numbers.h:21
#define nCopy(n)
Definition numbers.h:15
#define nGreater(a, b)
Definition numbers.h:28
#define nGreaterZero(n)
Definition numbers.h:27
#define omAllocBin(bin)
#define NULL
Definition omList.c:12
#define pDelete(p_ptr)
Definition polys.h:187
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition polys.h:68
#define pSetm(p)
Definition polys.h:272
#define pNSet(n)
Definition polys.h:314
#define pSub(a, b)
Definition polys.h:288
#define pGetExp(p, i)
Exponent.
Definition polys.h:42
#define pNormalize(p)
Definition polys.h:318
#define pEqualPolys(p1, p2)
Definition polys.h:400
#define pSetExp(p, i, v)
Definition polys.h:43
#define pCopy(p)
return a copy of the poly
Definition polys.h:186
#define pOne()
Definition polys.h:316
ideal idInit(int idsize, int rank)
initialise an ideal / module
#define IDELEMS(i)
@ INTVEC_CMD
Definition tok.h:101

◆ evHessenberg()

BOOLEAN evHessenberg ( leftv res,
leftv h )

Definition at line 101 of file eigenval_ip.cc.

102{
103 if(currRing)
104 {
105 if(h&&h->Typ()==MATRIX_CMD)
106 {
107 matrix M=(matrix)h->Data();
108 res->rtyp=MATRIX_CMD;
109 res->data=(void *)evHessenberg(mp_Copy(M, currRing));
110 return FALSE;
111 }
112 WerrorS("<matrix> expected");
113 return TRUE;
114 }
115 WerrorS("no ring active");
116 return TRUE;
117}
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition matpol.cc:57

◆ evRowElim()

BOOLEAN evRowElim ( leftv res,
leftv h )

Definition at line 51 of file eigenval_ip.cc.

52{
53 if(currRing)
54 {
55 const short t[]={4,MATRIX_CMD,INT_CMD,INT_CMD,INT_CMD};
56 if (iiCheckTypes(h,t,1))
57 {
58 matrix M=(matrix)h->CopyD();
59 h=h->next;
60 int i=(int)(long)h->Data();
61 h=h->next;
62 int j=(int)(long)h->Data();
63 h=h->next;
64 int k=(int)(long)h->Data();
65 res->rtyp=MATRIX_CMD;
66 res->data=(void *)evRowElim(M,i,j,k);
67 return FALSE;
68 }
69 return TRUE;
70 }
71 WerrorS("no ring active");
72 return TRUE;
73}
BOOLEAN evRowElim(leftv res, leftv h)
BOOLEAN iiCheckTypes(leftv args, const short *type_list, int report)
check a list of arguemys against a given field of types return TRUE if the types match return FALSE (...
Definition ipshell.cc:6569
@ INT_CMD
Definition tok.h:96

◆ evSwap()

BOOLEAN evSwap ( leftv res,
leftv h )

Definition at line 29 of file eigenval_ip.cc.

30{
31 if(currRing)
32 {
33 const short t[]={3,MATRIX_CMD,INT_CMD,INT_CMD};
34 if (iiCheckTypes(h,t,1))
35 {
36 matrix M=(matrix)h->Data();
37 h=h->next;
38 int i=(int)(long)h->Data();
39 h=h->next;
40 int j=(int)(long)h->Data();
41 res->rtyp=MATRIX_CMD;
42 res->data=(void *)evSwap(mp_Copy(M, currRing),i,j);
43 return FALSE;
44 }
45 return TRUE;
46 }
47 WerrorS("no ring active");
48 return TRUE;
49}
BOOLEAN evSwap(leftv res, leftv h)