source: git/modules/eigenval.mod @ 2dc9fea

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