source: git/kernel/eigenval.cc @ 76e501

spielwiese
Last change on this file since 76e501 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: 2.2 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/* $Id$ */
5/*
6* ABSTRACT: eigenvalues of constant square matrices
7*/
8
9#include <kernel/mod2.h>
10
11#ifdef HAVE_EIGENVAL
12
13#include <kernel/febase.h>
14#include <kernel/structs.h>
15//#include "ipid.h"
16#include <kernel/intvec.h>
17#include <kernel/numbers.h>
18#include <kernel/polys.h>
19#include <kernel/ideals.h>
20//#include "lists.h"
21#include <kernel/matpol.h>
22#include <kernel/clapsing.h>
23#include <kernel/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
48matrix evRowElim(matrix M,int i,int j,int k)
49{
50  if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
51    return(M);
52
53  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k))));
54  pNormalize(p);
55
56  for(int l=1;l<=MATCOLS(M);l++)
57  {
58    MATELEM(M,i,l)=pSub(MATELEM(M,i,l),ppMult_qq(p,MATELEM(M,j,l)));
59    pNormalize(MATELEM(M,i,l));
60  }
61  for(int l=1;l<=MATROWS(M);l++)
62  {
63    MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),ppMult_qq(p,MATELEM(M,l,i)));
64    pNormalize(MATELEM(M,l,j));
65  }
66
67  pDelete(&p);
68
69  return(M);
70}
71
72matrix evColElim(matrix M,int i,int j,int k)
73{
74  if(MATELEM(M,k,i)==0||MATELEM(M,k,j)==0)
75    return(M);
76
77  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j))));
78  pNormalize(p);
79
80  for(int l=1;l<=MATROWS(M);l++)
81  {
82    MATELEM(M,l,i)=pSub(MATELEM(M,l,i),ppMult_qq(p,MATELEM(M,l,j)));
83    pNormalize(MATELEM(M,l,i));
84  }
85  for(int l=1;l<=MATCOLS(M);l++)
86  {
87    MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),ppMult_qq(p,MATELEM(M,i,l)));
88    pNormalize(MATELEM(M,j,l));
89  }
90
91  pDelete(&p);
92
93  return(M);
94}
95
96matrix evHessenberg(matrix M)
97{
98  int n=MATROWS(M);
99  if(n!=MATCOLS(M))
100    return(M);
101
102  for(int k=1,j=2;k<n-1;k++,j=k+1)
103  {
104    while(j<=n&&MATELEM(M,j,k)==0)
105      j++;
106
107    if(j<=n)
108    {
109      M=evSwap(M,j,k+1);
110
111      for(int i=j+1;i<=n;i++)
112        M=evRowElim(M,i,k+1,k);
113    }
114  }
115
116  return(M);
117}
118
119#endif /* HAVE_EIGENVAL */
Note: See TracBrowser for help on using the repository browser.