My Project  debian-1:4.1.1-p2+ds-4build4
Functions
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 evColElim (leftv res, leftv h)
 
BOOLEAN evHessenberg (leftv res, leftv h)
 
lists evEigenvals (matrix M)
 
BOOLEAN evEigenvals (leftv res, leftv h)
 

Function Documentation

◆ evColElim()

BOOLEAN evColElim ( leftv  res,
leftv  h 
)

Definition at line 75 of file eigenval_ip.cc.

76 {
77  if(currRing)
78  {
79  const short t[]={4,MATRIX_CMD,INT_CMD,INT_CMD,INT_CMD};
80  if (iiCheckTypes(h,t,1))
81  {
82  matrix M=(matrix)h->Data();
83  h=h->next;
84  int i=(int)(long)h->Data();
85  h=h->next;
86  int j=(int)(long)h->Data();
87  h=h->next;
88  int k=(int)(long)h->Data();
89  res->rtyp=MATRIX_CMD;
90  res->data=(void *)evColElim(mp_Copy(M, currRing),i,j,k);
91  return FALSE;
92  }
93  return TRUE;
94  }
95  WerrorS("no ring active");
96  return TRUE;
97 }
#define TRUE
Definition: auxiliary.h:98
#define FALSE
Definition: auxiliary.h:94
int i
Definition: cfEzgcd.cc:125
int k
Definition: cfEzgcd.cc:92
BOOLEAN evColElim(leftv res, leftv h)
Definition: eigenval_ip.cc:75
CanonicalForm res
Definition: facAbsFact.cc:64
int j
Definition: facHensel.cc:105
void WerrorS(const char *s)
Definition: feFopen.cc:24
@ MATRIX_CMD
Definition: grammar.cc:285
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:6503
static Poly * h
Definition: janet.cc:972
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:63
ip_smatrix * matrix
Definition: matpol.h:31
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:24
@ INT_CMD
Definition: tok.h:96

◆ evEigenvals() [1/2]

BOOLEAN evEigenvals ( leftv  res,
leftv  h 
)

Definition at line 267 of file eigenval_ip.cc.

268 {
269  if(currRing)
270  {
271  if(h&&h->Typ()==MATRIX_CMD)
272  {
273  matrix M=(matrix)h->CopyD();
274  res->rtyp=LIST_CMD;
275  res->data=(void *)evEigenvals(M);
276  return FALSE;
277  }
278  WerrorS("<matrix> expected");
279  return TRUE;
280  }
281  WerrorS("no ring active");
282  return TRUE;
283 }
lists evEigenvals(matrix M)
Definition: eigenval_ip.cc:118
@ LIST_CMD
Definition: tok.h:118

◆ evEigenvals() [2/2]

lists evEigenvals ( matrix  M)

Definition at line 118 of file eigenval_ip.cc.

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

◆ evHessenberg()

BOOLEAN evHessenberg ( leftv  res,
leftv  h 
)

Definition at line 99 of file eigenval_ip.cc.

100 {
101  if(currRing)
102  {
103  if(h&&h->Typ()==MATRIX_CMD)
104  {
105  matrix M=(matrix)h->Data();
106  res->rtyp=MATRIX_CMD;
107  res->data=(void *)evHessenberg(mp_Copy(M, currRing));
108  return FALSE;
109  }
110  WerrorS("<matrix> expected");
111  return TRUE;
112  }
113  WerrorS("no ring active");
114  return TRUE;
115 }

◆ 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

◆ 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