My Project  debian-1:4.1.1-p2+ds-4build3
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 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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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 }
FALSE
#define FALSE
Definition: auxiliary.h:94
matrix
ip_smatrix * matrix
Definition: matpol.h:31
ip_smatrix
Definition: matpol.h:15
j
int j
Definition: facHensel.cc:105
k
int k
Definition: cfEzgcd.cc:92
idDelete
#define idDelete(H)
delete an ideal
Definition: ideals.h:29
MATELEM
#define MATELEM(mat, i, j)
Definition: matpol.h:30
nGreater
#define nGreater(a, b)
Definition: numbers.h:29
LIST_CMD
@ LIST_CMD
Definition: tok.h:118
pGetExp
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
singclap_factorize
ideal singclap_factorize(poly f, intvec **v, int with_exps, const ring r)
Definition: clapsing.cc:854
omAllocBin
#define omAllocBin(bin)
Definition: omAllocDecl.h:205
pDelete
#define pDelete(p_ptr)
Definition: polys.h:173
slists_bin
omBin slists_bin
Definition: lists.cc:23
evSwap
BOOLEAN evSwap(leftv res, leftv h)
Definition: eigenval_ip.cc:29
MATRIX_CMD
@ MATRIX_CMD
Definition: grammar.cc:285
nGreaterZero
#define nGreaterZero(n)
Definition: numbers.h:28
nInpNeg
#define nInpNeg(n)
Definition: numbers.h:22
currRing
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
TRUE
#define TRUE
Definition: auxiliary.h:98
i
int i
Definition: cfEzgcd.cc:125
res
CanonicalForm res
Definition: facAbsFact.cc:64
INT_CMD
@ INT_CMD
Definition: tok.h:96
M
#define M
Definition: sirandom.c:24
IDEAL_CMD
@ IDEAL_CMD
Definition: grammar.cc:283
h
static Poly * h
Definition: janet.cc:972
evColElim
BOOLEAN evColElim(leftv res, leftv h)
Definition: eigenval_ip.cc:75
pOne
#define pOne()
Definition: polys.h:301
intvec
Definition: intvec.h:21
nDiv
#define nDiv(a, b)
Definition: numbers.h:33
pNSet
#define pNSet(n)
Definition: polys.h:299
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:36
slists
Definition: lists.h:23
INTVEC_CMD
@ INTVEC_CMD
Definition: tok.h:101
evHessenberg
BOOLEAN evHessenberg(leftv res, leftv h)
Definition: eigenval_ip.cc:99
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:37
WerrorS
void WerrorS(const char *s)
Definition: feFopen.cc:24
pEqualPolys
#define pEqualPolys(p1, p2)
Definition: polys.h:386
m
int m
Definition: cfEzgcd.cc:121
MATCOLS
#define MATCOLS(i)
Definition: matpol.h:29
NULL
#define NULL
Definition: omList.c:10
lists
slists * lists
Definition: mpr_numeric.h:146
pSetm
#define pSetm(p)
Definition: polys.h:257
evEigenvals
lists evEigenvals(matrix M)
Definition: eigenval_ip.cc:118
l
int l
Definition: cfEzgcd.cc:93
nDelete
#define nDelete(n)
Definition: numbers.h:17
mp_DetBareiss
poly mp_DetBareiss(matrix a, const ring r)
returns the determinant of the matrix m; uses Bareiss algorithm
Definition: matpol.cc:1576
pSetExp
#define pSetExp(p, i, v)
Definition: polys.h:42
iiCheckTypes
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
pCopy
#define pCopy(p)
return a copy of the poly
Definition: polys.h:172
IDELEMS
#define IDELEMS(i)
Definition: simpleideals.h:26
pNormalize
#define pNormalize(p)
Definition: polys.h:303
pHead
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
pGetCoeff
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
mp_Copy
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition: matpol.cc:63
MATROWS
#define MATROWS(i)
Definition: matpol.h:28
pSub
#define pSub(a, b)
Definition: polys.h:273
nCopy
#define nCopy(n)
Definition: numbers.h:16
evRowElim
BOOLEAN evRowElim(leftv res, leftv h)
Definition: eigenval_ip.cc:51
pNext
#define pNext(p)
Definition: monomials.h:43