source: git/kernel/eigenval.cc @ 6d77472

spielwiese
Last change on this file since 6d77472 was ba5e9e, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Changed configure-scripts to generate individual public config files for each package: resources, libpolys, singular (main) fix: sources should include correct corresponding config headers.
  • Property mode set to 100644
File size: 2.3 KB
Line 
1/*****************************************
2*  Computer Algebra System SINGULAR      *
3*****************************************/
4/*
5* ABSTRACT: eigenvalues of constant square matrices
6*/
7
8#ifdef HAVE_CONFIG_H
9#include "singularconfig.h"
10#endif /* HAVE_CONFIG_H */
11#include <kernel/mod2.h>
12
13#ifdef HAVE_EIGENVAL
14
15#include <kernel/febase.h>
16#include <kernel/structs.h>
17//#include "ipid.h"
18#include <misc/intvec.h>
19#include <coeffs/numbers.h>
20#include <kernel/polys.h>
21#include <kernel/ideals.h>
22//#include "lists.h"
23#include <polys/matpol.h>
24#include <polys/clapsing.h>
25#include <kernel/eigenval.h>
26
27
28matrix evSwap(matrix M,int i,int j)
29{
30  if(i==j)
31    return(M);
32
33  for(int k=1;k<=MATROWS(M);k++)
34  {
35    poly p=MATELEM(M,i,k);
36    MATELEM(M,i,k)=MATELEM(M,j,k);
37    MATELEM(M,j,k)=p;
38  }
39
40  for(int k=1;k<=MATCOLS(M);k++)
41  {
42    poly p=MATELEM(M,k,i);
43    MATELEM(M,k,i)=MATELEM(M,k,j);
44    MATELEM(M,k,j)=p;
45  }
46
47  return(M);
48}
49
50matrix evRowElim(matrix M,int i,int j,int k)
51{
52  if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
53    return(M);
54
55  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k))));
56  pNormalize(p);
57
58  for(int l=1;l<=MATCOLS(M);l++)
59  {
60    MATELEM(M,i,l)=pSub(MATELEM(M,i,l),ppMult_qq(p,MATELEM(M,j,l)));
61    pNormalize(MATELEM(M,i,l));
62  }
63  for(int l=1;l<=MATROWS(M);l++)
64  {
65    MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),ppMult_qq(p,MATELEM(M,l,i)));
66    pNormalize(MATELEM(M,l,j));
67  }
68
69  pDelete(&p);
70
71  return(M);
72}
73
74matrix evColElim(matrix M,int i,int j,int k)
75{
76  if(MATELEM(M,k,i)==0||MATELEM(M,k,j)==0)
77    return(M);
78
79  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,k,i)),pGetCoeff(MATELEM(M,k,j))));
80  pNormalize(p);
81
82  for(int l=1;l<=MATROWS(M);l++)
83  {
84    MATELEM(M,l,i)=pSub(MATELEM(M,l,i),ppMult_qq(p,MATELEM(M,l,j)));
85    pNormalize(MATELEM(M,l,i));
86  }
87  for(int l=1;l<=MATCOLS(M);l++)
88  {
89    MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),ppMult_qq(p,MATELEM(M,i,l)));
90    pNormalize(MATELEM(M,j,l));
91  }
92
93  pDelete(&p);
94
95  return(M);
96}
97
98matrix evHessenberg(matrix M)
99{
100  int n=MATROWS(M);
101  if(n!=MATCOLS(M))
102    return(M);
103
104  for(int k=1,j=2;k<n-1;k++,j=k+1)
105  {
106    while(j<=n&&MATELEM(M,j,k)==0)
107      j++;
108
109    if(j<=n)
110    {
111      M=evSwap(M,j,k+1);
112
113      for(int i=j+1;i<=n;i++)
114        M=evRowElim(M,i,k+1,k);
115    }
116  }
117
118  return(M);
119}
120
121#endif /* HAVE_EIGENVAL */
Note: See TracBrowser for help on using the repository browser.