source: git/Singular/eigenval.cc @ a7c2db

spielwiese
Last change on this file since a7c2db was a7c2db, checked in by Mathias Schulze <mschulze@…>, 22 years ago
*** empty log message *** git-svn-id: file:///usr/local/Singular/svn/trunk@5902 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 7.2 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* $Id: eigenval.cc,v 1.10 2002-02-20 15:50:34 mschulze Exp $ */
5/*
6* ABSTRACT: eigenvalues of constant square matrices
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_EIGENVAL
12
13#include "febase.h"
14#include "tok.h"
15#include "ipid.h"
16#include "intvec.h"
17#include "numbers.h"
18#include "polys.h"
19#include "ideals.h"
20#include "lists.h"
21#include "matpol.h"
22#include "clapsing.h"
23#include "eigenval.h"
24
25
26matrix evSwap(matrix M,int i,int j)
27{
28  if(i==j)
29    return(M);
30
31  for(int k=1;k<=MATROWS(M);k++)
32  {
33    poly p=MATELEM(M,i,k);
34    MATELEM(M,i,k)=MATELEM(M,j,k);
35    MATELEM(M,j,k)=p;
36  }
37
38  for(int k=1;k<=MATCOLS(M);k++)
39  {
40    poly p=MATELEM(M,k,i);
41    MATELEM(M,k,i)=MATELEM(M,k,j);
42    MATELEM(M,k,j)=p;
43  }
44
45  return(M);
46}
47
48
49BOOLEAN evSwap(leftv res,leftv h)
50{
51  if(currRingHdl)
52  {
53    if(h&&h->Typ()==MATRIX_CMD)
54    {
55      matrix M=(matrix)h->Data();
56      h=h->next;
57      if(h&&h->Typ()==INT_CMD)
58      {
59        int i=(int)h->Data();
60        h=h->next;
61        if(h&&h->Typ()==INT_CMD)
62        {
63          int j=(int)h->Data();
64          res->rtyp=MATRIX_CMD;
65          res->data=(void *)evSwap(mpCopy(M),i,j);
66          return FALSE;
67        }
68      }
69    }
70    WerrorS("<matrix>,<int>,<int> expected");
71    return TRUE;
72  }
73  WerrorS("no ring active");
74  return TRUE;
75}
76
77
78matrix evRowElim(matrix M,int i,int j,int k)
79{
80  if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
81    return(M);
82
83  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k))));
84
85  for(int l=1;l<=MATCOLS(M);l++)
86    MATELEM(M,i,l)=pSub(MATELEM(M,i,l),ppMult_qq(p,MATELEM(M,j,l)));
87
88  for(int l=1;l<=MATROWS(M);l++)
89    MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),ppMult_qq(p,MATELEM(M,l,i)));
90
91  pDelete(&p);
92
93  return(M);
94}
95
96
97BOOLEAN evRowElim(leftv res,leftv h)
98{
99  if(currRingHdl)
100  {
101    if(h&&h->Typ()==MATRIX_CMD)
102    {
103      matrix M=(matrix)h->Data();
104      h=h->next;
105      if(h&&h->Typ()==INT_CMD)
106      {
107        int i=(int)h->Data();
108        h=h->next;
109        if(h&&h->Typ()==INT_CMD)
110        {
111          int j=(int)h->Data();
112          h=h->next;
113          if(h&&h->Typ()==INT_CMD)
114          {
115            int k=(int)h->Data();
116            res->rtyp=MATRIX_CMD;
117            res->data=(void *)evRowElim(mpCopy(M),i,j,k);
118            return FALSE;
119          }
120        }
121      }
122    }
123    WerrorS("<matrix>,<int>,<int>,<int> expected");
124    return TRUE;
125  }
126  WerrorS("no ring active");
127  return TRUE;
128}
129
130
131matrix evColElim(matrix M,int i,int j,int k)
132{
133  if(MATELEM(M,k,i)==0||MATELEM(M,k,j)==0)
134    return(M);
135
136  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j))));
137
138  for(int l=1;l<=MATROWS(M);l++)
139    MATELEM(M,l,i)=pSub(MATELEM(M,l,i),ppMult_qq(p,MATELEM(M,l,j)));
140
141  for(int l=1;l<=MATCOLS(M);l++)
142    MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),ppMult_qq(p,MATELEM(M,i,l)));
143
144  pDelete(&p);
145
146  return(M);
147}
148
149
150BOOLEAN evColElim(leftv res,leftv h)
151{
152  if(currRingHdl)
153  {
154    if(h&&h->Typ()==MATRIX_CMD)
155    {
156      matrix M=(matrix)h->Data();
157      h=h->next;
158      if(h&&h->Typ()==INT_CMD)
159      {
160        int i=(int)h->Data();
161        h=h->next;
162        if(h&&h->Typ()==INT_CMD)
163        {
164          int j=(int)h->Data();
165          h=h->next;
166          if(h&&h->Typ()==INT_CMD)
167          {
168            int k=(int)h->Data();
169            res->rtyp=MATRIX_CMD;
170            res->data=(void *)evColElim(mpCopy(M),i,j,k);
171            return FALSE;
172          }
173        }
174      }
175    }
176    WerrorS("<matrix>,<int>,<int>,<int> expected");
177    return TRUE;
178  }
179  WerrorS("no ring active");
180  return TRUE;
181}
182
183
184matrix evHessenberg(matrix M)
185{
186  int n=MATROWS(M);
187
188  for(int k=1,j=2;k<n-1;k++,j=k+1)
189  {
190    while(j<=n&&MATELEM(M,j,k)==0)
191      j++;
192
193    if(j<=n)
194    {
195      M=evSwap(M,j,k+1);
196
197      for(int i=j+1;i<=n;i++)
198        M=evRowElim(M,i,k+1,k);
199    }
200  }
201
202  return(M);
203}
204
205
206BOOLEAN evHessenberg(leftv res,leftv h)
207{
208  if(currRingHdl)
209  {
210    if(h&&h->Typ()==MATRIX_CMD)
211    {
212      matrix M=(matrix)h->Data();
213      res->rtyp=MATRIX_CMD;
214      res->data=(void *)evHessenberg(mpCopy(M));
215      return FALSE;
216    }
217    WerrorS("<matrix> expected");
218    return TRUE;
219  }
220  WerrorS("no ring active");
221  return TRUE;
222}
223
224
225lists evEigenvals(matrix M)
226{
227  lists l=(lists)omAllocBin(slists_bin);
228  if(MATROWS(M)!=MATCOLS(M))
229  {
230    l->Init(0);
231    return(l);
232  }
233
234  M=evHessenberg((matrix)idJet((ideal)M,0));
235
236  int n=MATROWS(M);
237  ideal e=idInit(n,1);
238  intvec *m=new intvec(n);
239
240  poly t=pOne();
241  pSetExp(t,1,1);
242  pSetm(t);
243
244  for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
245  {
246    while(j<=n&&MATELEM(M,j,j-1)!=NULL)
247      j++;
248    if(j==j0+1)
249    {
250      e->m[k]=pHead(MATELEM(M,j0,j0));
251      (*m)[k]=1;
252      k++;
253    }
254    else
255    {
256      int n0=j-j0;
257      matrix M0=mpNew(n0,n0);
258
259      j0--;
260      for(int i=1;i<=n0;i++)
261        for(int j=1;j<=n0;j++)
262          MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
263      for(int i=1;i<=n0;i++)
264        MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
265
266      intvec *m0;
267      ideal e0=singclap_factorize(mpDetBareiss(M0),&m0,2);
268
269      for(int i=0;i<IDELEMS(e0);i++)
270      {
271        if(pNext(e0->m[i])==NULL)
272        {
273          (*m)[k]=(*m0)[i];
274          k++;
275        }
276        else
277        if(pGetExp(e0->m[i],1)<2&&pGetExp(pNext(e0->m[i]),1)<2&&
278           pNext(pNext(e0->m[i]))==NULL)
279        {
280          number e1=nNeg(pGetCoeff(e0->m[i]));
281          if(pGetExp(pNext(e0->m[i]),1)==0)
282            e->m[k]=pNSet(nDiv(pGetCoeff(pNext(e0->m[i])),e1));
283          else
284            e->m[k]=pNSet(nDiv(e1,pGetCoeff(pNext(e0->m[i]))));
285          nDelete(&e1);
286          (*m)[k]=(*m0)[i];
287          k++;
288        }
289        else
290        {
291          e->m[k]=e0->m[i];
292          e0->m[i]=NULL;
293          (*m)[k]=(*m0)[i];
294          k++;
295        }
296      }
297
298      delete(m0);
299      idDelete(&e0);
300    }
301  }
302
303  pDelete(&t);
304  idDelete((ideal *)&M);
305
306  for(int i0=0;i0<n-1;i0++)
307  {
308    for(int i1=i0+1;i1<n;i1++)
309    {
310      if(pEqualPolys(e->m[i0],e->m[i1]))
311      {
312        (*m)[i0]+=(*m)[i1];
313        (*m)[i1]=0;
314      }
315      else
316      {
317        if(e->m[i0]==NULL&&!nGreaterZero(pGetCoeff(e->m[i1]))||
318           e->m[i1]==NULL&&
319          (nGreaterZero(pGetCoeff(e->m[i0]))||pNext(e->m[i0])!=NULL)||
320           e->m[i0]!=NULL&&e->m[i1]!=NULL&&
321          (pNext(e->m[i0])!=NULL&&pNext(e->m[i1])==NULL||
322           pNext(e->m[i0])==NULL&&pNext(e->m[i1])==NULL&&
323           nGreater(pGetCoeff(e->m[i0]),pGetCoeff(e->m[i1]))))
324        {
325          poly e1=e->m[i0];
326          e->m[i0]=e->m[i1];
327          e->m[i1]=e1;
328          int m1=(*m)[i0];
329          (*m)[i0]=(*m)[i1];
330          (*m)[i1]=m1;
331        }
332      }
333    }
334  }
335
336  int n0=0;
337  for(int i=0;i<n;i++)
338    if((*m)[i]>0)
339      n0++;
340
341  ideal e0=idInit(n0,1);
342  intvec *m0=new intvec(n0);
343
344  for(int i=0,i0=0;i<n;i++)
345    if((*m)[i]>0)
346    {
347      e0->m[i0]=e->m[i];
348      e->m[i]=NULL;
349      (*m0)[i0]=(*m)[i];
350      i0++;
351    }
352
353  idDelete(&e);
354  delete(m);
355
356  l->Init(2);
357  l->m[0].rtyp=IDEAL_CMD;
358  l->m[0].data=e0;
359  l->m[1].rtyp=INTVEC_CMD;
360  l->m[1].data=m0;
361
362  return(l);
363}
364
365
366BOOLEAN evEigenvals(leftv res,leftv h)
367{
368  if(currRingHdl)
369  {
370    if(h&&h->Typ()==MATRIX_CMD)
371    {
372      matrix M=(matrix)h->Data();
373      res->rtyp=LIST_CMD;
374      res->data=(void *)evEigenvals(mpCopy(M));
375      return FALSE;
376    }
377    WerrorS("<matrix> expected");
378    return TRUE;
379  }
380  WerrorS("no ring active");
381  return TRUE;
382}
383
384#endif /* HAVE_EIGENVAL */
Note: See TracBrowser for help on using the repository browser.