My Project
Functions
eigenval_ip.h File Reference
#include "kernel/linear_algebra/eigenval.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:100
#define FALSE
Definition: auxiliary.h:96
lists evEigenvals(matrix M)
Definition: eigenval_ip.cc:120
CanonicalForm res
Definition: facAbsFact.cc:60
void WerrorS(const char *s)
Definition: feFopen.cc:24
@ MATRIX_CMD
Definition: grammar.cc:286
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 
129  M=evHessenberg(M);
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;
162  ideal e0=singclap_factorize(mp_DetBareiss(M0,currRing),&m0,2, currRing);
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
Definition: intvec.h:23
Definition: lists.h:24
BOOLEAN evHessenberg(leftv res, leftv h)
Definition: eigenval_ip.cc:101
int j
Definition: facHensel.cc:110
@ IDEAL_CMD
Definition: grammar.cc:284
#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:1676
#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
Definition: mpr_numeric.h:146
#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)
Definition: omAllocDecl.h:205
#define NULL
Definition: omList.c:12
#define pDelete(p_ptr)
Definition: polys.h:186
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:271
#define pNSet(n)
Definition: polys.h:313
#define pSub(a, b)
Definition: polys.h:287
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
#define pNormalize(p)
Definition: polys.h:317
#define pEqualPolys(p1, p2)
Definition: polys.h:399
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pCopy(p)
return a copy of the poly
Definition: polys.h:185
#define pOne()
Definition: polys.h:315
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define IDELEMS(i)
Definition: simpleideals.h:23
@ 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:64

◆ 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)
Definition: eigenval_ip.cc:51
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:6572
@ 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)
Definition: eigenval_ip.cc:29