source: git/Singular/eigenval_ip.cc @ 599326

spielwiese
Last change on this file since 599326 was 599326, checked in by Kai Krüger <krueger@…>, 14 years ago
Anne, Kai, Frank: - changes to #include "..." statements to allow cleaner build structure - affected directories: omalloc, kernel, Singular - not yet done: IntergerProgramming git-svn-id: file:///usr/local/Singular/svn/trunk@13032 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 6.1 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: eigenvalues of constant square matrices
7*/
8
9#include <Singular/mod2.h>
10
11#ifdef HAVE_EIGENVAL
12
13#include <kernel/febase.h>
14#include <Singular/tok.h>
15#include <Singular/ipid.h>
16#include <kernel/intvec.h>
17#include <kernel/numbers.h>
18#include <kernel/polys.h>
19#include <kernel/ideals.h>
20#include <Singular/lists.h>
21#include <kernel/matpol.h>
22#include <Singular/clapsing.h>
23#include <Singular/eigenval.h>
24#include <Singular/eigenval_ip.h>
25
26
27BOOLEAN evSwap(leftv res,leftv h)
28{
29  if(currRingHdl)
30  {
31    if(h&&h->Typ()==MATRIX_CMD)
32    {
33      matrix M=(matrix)h->Data();
34      h=h->next;
35      if(h&&h->Typ()==INT_CMD)
36      {
37        int i=(int)(long)h->Data();
38        h=h->next;
39        if(h&&h->Typ()==INT_CMD)
40        {
41          int j=(int)(long)h->Data();
42          res->rtyp=MATRIX_CMD;
43          res->data=(void *)evSwap(mpCopy(M),i,j);
44          return FALSE;
45        }
46      }
47    }
48    WerrorS("<matrix>,<int>,<int> expected");
49    return TRUE;
50  }
51  WerrorS("no ring active");
52  return TRUE;
53}
54
55BOOLEAN evRowElim(leftv res,leftv h)
56{
57  if(currRingHdl)
58  {
59    if(h&&h->Typ()==MATRIX_CMD)
60    {
61      matrix M=(matrix)h->Data();
62      h=h->next;
63      if(h&&h->Typ()==INT_CMD)
64      {
65        int i=(int)(long)h->Data();
66        h=h->next;
67        if(h&&h->Typ()==INT_CMD)
68        {
69          int j=(int)(long)h->Data();
70          h=h->next;
71          if(h&&h->Typ()==INT_CMD)
72          {
73            int k=(int)(long)h->Data();
74            res->rtyp=MATRIX_CMD;
75            res->data=(void *)evRowElim(mpCopy(M),i,j,k);
76            return FALSE;
77          }
78        }
79      }
80    }
81    WerrorS("<matrix>,<int>,<int>,<int> expected");
82    return TRUE;
83  }
84  WerrorS("no ring active");
85  return TRUE;
86}
87
88BOOLEAN evColElim(leftv res,leftv h)
89{
90  if(currRingHdl)
91  {
92    if(h&&h->Typ()==MATRIX_CMD)
93    {
94      matrix M=(matrix)h->Data();
95      h=h->next;
96      if(h&&h->Typ()==INT_CMD)
97      {
98        int i=(int)(long)h->Data();
99        h=h->next;
100        if(h&&h->Typ()==INT_CMD)
101        {
102          int j=(int)(long)h->Data();
103          h=h->next;
104          if(h&&h->Typ()==INT_CMD)
105          {
106            int k=(int)(long)h->Data();
107            res->rtyp=MATRIX_CMD;
108            res->data=(void *)evColElim(mpCopy(M),i,j,k);
109            return FALSE;
110          }
111        }
112      }
113    }
114    WerrorS("<matrix>,<int>,<int>,<int> expected");
115    return TRUE;
116  }
117  WerrorS("no ring active");
118  return TRUE;
119}
120
121BOOLEAN evHessenberg(leftv res,leftv h)
122{
123  if(currRingHdl)
124  {
125    if(h&&h->Typ()==MATRIX_CMD)
126    {
127      matrix M=(matrix)h->Data();
128      res->rtyp=MATRIX_CMD;
129      res->data=(void *)evHessenberg(mpCopy(M));
130      return FALSE;
131    }
132    WerrorS("<matrix> expected");
133    return TRUE;
134  }
135  WerrorS("no ring active");
136  return TRUE;
137}
138
139
140lists evEigenvals(matrix M)
141{
142  lists l=(lists)omAllocBin(slists_bin);
143  if(MATROWS(M)!=MATCOLS(M))
144  {
145    l->Init(0);
146    return(l);
147  }
148
149  M=evHessenberg((matrix)idJet((ideal)M,0));
150
151  int n=MATROWS(M);
152  ideal e=idInit(n,1);
153  intvec *m=new intvec(n);
154
155  poly t=pOne();
156  pSetExp(t,1,1);
157  pSetm(t);
158
159  for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
160  {
161    while(j<=n&&MATELEM(M,j,j-1)!=NULL)
162      j++;
163    if(j==j0+1)
164    {
165      e->m[k]=pHead(MATELEM(M,j0,j0));
166      (*m)[k]=1;
167      k++;
168    }
169    else
170    {
171      int n0=j-j0;
172      matrix M0=mpNew(n0,n0);
173
174      j0--;
175      for(int i=1;i<=n0;i++)
176        for(int j=1;j<=n0;j++)
177          MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
178      for(int i=1;i<=n0;i++)
179        MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
180
181      intvec *m0;
182      ideal e0=singclap_factorize(mpDetBareiss(M0),&m0,2);
183      if (e0==NULL)
184      {
185        l->Init(0);
186        return(l);
187      }
188
189      for(int i=0;i<IDELEMS(e0);i++)
190      {
191        if(pNext(e0->m[i])==NULL)
192        {
193          (*m)[k]=(*m0)[i];
194          k++;
195        }
196        else
197        if(pGetExp(e0->m[i],1)<2&&pGetExp(pNext(e0->m[i]),1)<2&&
198           pNext(pNext(e0->m[i]))==NULL)
199        {
200          number e1=nNeg(nCopy(pGetCoeff(e0->m[i])));
201          if(pGetExp(pNext(e0->m[i]),1)==0)
202            e->m[k]=pNSet(nDiv(pGetCoeff(pNext(e0->m[i])),e1));
203          else
204            e->m[k]=pNSet(nDiv(e1,pGetCoeff(pNext(e0->m[i]))));
205          nDelete(&e1);
206          pNormalize(e->m[k]);
207          (*m)[k]=(*m0)[i];
208          k++;
209        }
210        else
211        {
212          e->m[k]=e0->m[i];
213          pNormalize(e->m[k]);
214          e0->m[i]=NULL;
215          (*m)[k]=(*m0)[i];
216          k++;
217        }
218      }
219
220      delete(m0);
221      idDelete(&e0);
222    }
223  }
224
225  pDelete(&t);
226  idDelete((ideal *)&M);
227
228  for(int i0=0;i0<n-1;i0++)
229  {
230    for(int i1=i0+1;i1<n;i1++)
231    {
232      if(pEqualPolys(e->m[i0],e->m[i1]))
233      {
234        (*m)[i0]+=(*m)[i1];
235        (*m)[i1]=0;
236      }
237      else
238      {
239        if(e->m[i0]==NULL&&!nGreaterZero(pGetCoeff(e->m[i1]))||
240           e->m[i1]==NULL&&
241          (nGreaterZero(pGetCoeff(e->m[i0]))||pNext(e->m[i0])!=NULL)||
242           e->m[i0]!=NULL&&e->m[i1]!=NULL&&
243          (pNext(e->m[i0])!=NULL&&pNext(e->m[i1])==NULL||
244           pNext(e->m[i0])==NULL&&pNext(e->m[i1])==NULL&&
245           nGreater(pGetCoeff(e->m[i0]),pGetCoeff(e->m[i1]))))
246        {
247          poly e1=e->m[i0];
248          e->m[i0]=e->m[i1];
249          e->m[i1]=e1;
250          int m1=(*m)[i0];
251          (*m)[i0]=(*m)[i1];
252          (*m)[i1]=m1;
253        }
254      }
255    }
256  }
257
258  int n0=0;
259  for(int i=0;i<n;i++)
260    if((*m)[i]>0)
261      n0++;
262
263  ideal e0=idInit(n0,1);
264  intvec *m0=new intvec(n0);
265
266  for(int i=0,i0=0;i<n;i++)
267    if((*m)[i]>0)
268    {
269      e0->m[i0]=e->m[i];
270      e->m[i]=NULL;
271      (*m0)[i0]=(*m)[i];
272      i0++;
273    }
274
275  idDelete(&e);
276  delete(m);
277
278  l->Init(2);
279  l->m[0].rtyp=IDEAL_CMD;
280  l->m[0].data=e0;
281  l->m[1].rtyp=INTVEC_CMD;
282  l->m[1].data=m0;
283
284  return(l);
285}
286
287
288BOOLEAN evEigenvals(leftv res,leftv h)
289{
290  if(currRingHdl)
291  {
292    if(h&&h->Typ()==MATRIX_CMD)
293    {
294      matrix M=(matrix)h->Data();
295      res->rtyp=LIST_CMD;
296      res->data=(void *)evEigenvals(mpCopy(M));
297      return FALSE;
298    }
299    WerrorS("<matrix> expected");
300    return TRUE;
301  }
302  WerrorS("no ring active");
303  return TRUE;
304}
305
306#endif /* HAVE_EIGENVAL */
Note: See TracBrowser for help on using the repository browser.