Changeset a7f2a04 in git


Ignore:
Timestamp:
Feb 16, 2002, 11:56:19 AM (21 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
Children:
03495f1856f489d63ec19f337bd5240315ef5371
Parents:
502ea23f40b3a43fb8b3bcbc306d20b0c02ab50b
Message:
*mschulze: new procedures for eigenvalues


git-svn-id: file:///usr/local/Singular/svn/trunk@5877 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/eigenval.cc

    r502ea2 ra7f2a04  
    1 /****************************************
    2 *  Computer Algebra System SINGULAR     *
    3 ****************************************/
    4 /* $Id: eigenval.cc,v 1.5 2002-01-19 14:48:14 obachman Exp $ */
     1/*****************************************
     2*  Computer Algebra System SINGULAR      *
     3*****************************************/
     4/* $Id: eigenval.cc,v 1.6 2002-02-16 10:56:14 mschulze Exp $ */
    55/*
    66* ABSTRACT: eigenvalues of constant square matrices
    … …  
    1111#ifdef HAVE_EIGENVAL
    1212
     13#include "febase.h"
    1314#include "tok.h"
    14 #include "febase.h"
     15#include "ipid.h"
     16#include "intvec.h"
    1517#include "numbers.h"
    1618#include "polys.h"
    1719#include "ideals.h"
     20#include "lists.h"
    1821#include "matpol.h"
    1922#include "clapsing.h"
    20 
    21 matrix swap(matrix M,int i,int j)
     23#include "eigenval.h"
     24
     25
     26matrix evSwap(matrix M,int i,int j)
    2227{
    2328  if(i==j)
    24     return M;
    25   poly p;
     29    return(M);
     30
    2631  for(int k=1;k<=MATROWS(M);k++)
    2732  {
    28     p=MATELEM(M,i,k);
     33    poly p=MATELEM(M,i,k);
    2934    MATELEM(M,i,k)=MATELEM(M,j,k);
    3035    MATELEM(M,j,k)=p;
    3136  }
     37
    3238  for(int k=1;k<=MATCOLS(M);k++)
    3339  {
    34     p=MATELEM(M,k,i);
     40    poly p=MATELEM(M,k,i);
    3541    MATELEM(M,k,i)=MATELEM(M,k,j);
    3642    MATELEM(M,k,j)=p;
    3743  }
     44
    3845  return(M);
    3946}
    4047
    41 matrix rowelim(matrix M,int i, int j, int k)
     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)
    4279{
    4380  if(MATELEM(M,i,k)==NULL||MATELEM(M,j,k)==NULL)
    44     return M;
    45   number n=nDiv(pGetCoeff(MATELEM(M,i,k)),
    46                 pGetCoeff(MATELEM(M,j,k)));
     81    return(M);
     82
     83  poly p=pNSet(nDiv(pGetCoeff(MATELEM(M,i,k)),pGetCoeff(MATELEM(M,j,k))));
     84
    4785  for(int l=1;l<=MATCOLS(M);l++)
    48   {
    49     MATELEM(M,i,l)=pSub(MATELEM(M,i,l),pMult_nn(pCopy(MATELEM(M,j,l)),n));
    50   }
     86    MATELEM(M,i,l)=pSub(MATELEM(M,i,l),pMult(pCopy(p),pCopy(MATELEM(M,j,l))));
     87
    5188  for(int l=1;l<=MATROWS(M);l++)
    52   {
    53     MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),pMult_nn(pCopy(MATELEM(M,l,i)),n));
    54   }
    55   nDelete(&n);
    56   return M;
    57 }
    58 
    59 matrix colelim(matrix M,int i, int j, int k)
    60 {
    61   if(MATELEM(M,k,i)==NULL||MATELEM(M,k,j)==NULL)
    62     return M;
    63   number n=nDiv(pGetCoeff(MATELEM(M,k,i)),
    64                 pGetCoeff(MATELEM(M,k,j)));
     89    MATELEM(M,l,j)=pAdd(MATELEM(M,l,j),pMult(pCopy(p),pCopy(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
    65138  for(int l=1;l<=MATROWS(M);l++)
    66   {
    67     MATELEM(M,l,i)=pSub(MATELEM(M,l,i),pMult_nn(pCopy(MATELEM(M,l,j)),n));
    68   }
     139    MATELEM(M,l,i)=pSub(MATELEM(M,l,i),pMult(pCopy(p),pCopy(MATELEM(M,l,j))));
     140
    69141  for(int l=1;l<=MATCOLS(M);l++)
    70   {
    71     MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),pMult_nn(pCopy(MATELEM(M,i,l)),n));
    72   }
    73   nDelete(&n);
    74   return M;
    75 }
    76 
    77 matrix tridiag(matrix M)
    78 {
    79   int n=MATCOLS(M);
    80   for(int k=1;k<=n-2;k++)
    81   {
    82     int j=k+1;
    83     while(j<=n&&MATELEM(M,j,k)==NULL) j++;
    84     if(j<=n)
    85     {
    86       for(int i=j+1;i<=n;i++)
    87       {
    88         M=rowelim(M,i,j,k);
    89       }
    90       M=swap(M,j,k+1);
    91     }
     142    MATELEM(M,j,l)=pAdd(MATELEM(M,j,l),pMult(pCopy(p),pCopy(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  int i,j;
     188
     189  for(int k=1;k<n-1;k++)
     190  {
    92191    j=k+1;
    93     while(j<=n&&MATELEM(M,k,j)==NULL) j++;
    94     if(j<=n)
    95     {
    96       for(int i=j+1;i<=n;i++)
    97       {
    98         M=colelim(M,i,j,k);
    99       }
    100       M=swap(M,j,k+1);
    101     }
    102   }
     192    while(j<n&&MATELEM(M,j,k)==0)
     193      j++;
     194
     195    if(MATELEM(M,j,k)!=0)
     196    {
     197      M=evSwap(M,j,k+1);
     198
     199      for(i=j+1;i<=n;i++)
     200        M=evRowElim(M,i,k+1,k);
     201    }
     202  }
     203
    103204  return(M);
    104205}
    105206
    106 lists addval(lists l,poly e0,int m0)
    107 {
    108   ideal ee=(ideal)l->m[0].data;
    109   intvec *mm=(intvec*)l->m[1].data;
    110   int n=0;
    111   if(ee!=NULL)
    112     n=IDELEMS(ee);
    113   for(int i=n-1;i>=0;i--)
    114   {
    115     if(pEqualPolys(ee->m[i],e0))
    116     {
    117       (*mm)[i]+=m0;
    118       return l;
    119     }
    120   }
    121   ideal e=idInit(n+1,1);
    122   for(int i=n-1;i>=0;i--)
    123   {
    124     e->m[i]=ee->m[i];
    125     ee->m[i]=NULL;
    126   }
    127   e->m[n]=e0;
    128   l->m[0].data=e;
    129   if(ee!=NULL)
    130     idDelete(&ee);
    131   mm->resize(n+1);
    132   (*mm)[n]=m0;
    133   return l;
    134 }
    135 
    136 lists sortval(lists l)
    137 {
    138   ideal ee=(ideal)l->m[0].data;
    139   intvec *mm=(intvec*)l->m[1].data;
    140   int n=IDELEMS(ee);
    141   for(int i=n-1;i>=1;i--)
    142     for(int j=i-1;j>=0;j--)
    143       if(nGreater(pGetCoeff(ee->m[j]),pGetCoeff(ee->m[i])))
    144       {
    145         poly e=ee->m[i];
    146         ee->m[i]=ee->m[j];
    147         ee->m[j]=e;
    148         int m=(*mm)[i];
    149         (*mm)[i]=(*mm)[j];
    150         (*mm)[j]=m;
    151       }
    152   return l;
    153 }
    154 
    155 lists eigenval(matrix M)
    156 {
    157   M=tridiag(M);
    158   int n=MATCOLS(M);
     207
     208BOOLEAN evHessenberg(leftv res,leftv h)
     209{
     210  if(currRingHdl)
     211  {
     212    if(h&&h->Typ()==MATRIX_CMD)
     213    {
     214      matrix M=(matrix)h->Data();
     215      res->rtyp=MATRIX_CMD;
     216      res->data=(void*)evHessenberg(mpCopy(M));
     217      return FALSE;
     218    }
     219    WerrorS("<matrix> expected");
     220    return TRUE;
     221  }
     222  WerrorS("no ring active");
     223  return TRUE;
     224}
     225
     226
     227lists evEigenvalue(matrix M)
     228{
    159229  lists l=(lists)omAllocBin(slists_bin);
     230  if(MATROWS(M)!=MATCOLS(M))
     231  {
     232    l->Init(0);
     233    return(l);
     234  }
     235
     236  M=evHessenberg(M);
     237
     238  int n=MATROWS(M);
     239  ideal e=idInit(n,1);
     240  intvec *m=new intvec(n);
     241
     242  poly t=pOne();
     243  pSetExp(t,1,1);
     244
     245  for(int j0=1,j=2,k=0;j<=n+1;j0=j,j++)
     246  {
     247    while(j<=n&&MATELEM(M,j,j-1)!=NULL)
     248      j++;
     249    if(j==j0+1)
     250    {
     251      e->m[k]=pHead(MATELEM(M,j0,j0));
     252      (*m)[k]=1;
     253      k++;
     254    }
     255    else
     256    {
     257      int n0=j-j0;
     258      matrix M0=mpNew(n0,n0);
     259
     260      j0--;
     261      for(int i=1;i<=n0;i++)
     262        for(int j=1;j<=n0;j++)
     263          MATELEM(M0,i,j)=pCopy(MATELEM(M,j0+i,j0+j));
     264      for(int i=1;i<=n0;i++)
     265        MATELEM(M0,i,i)=pSub(MATELEM(M0,i,i),pCopy(t));
     266
     267      intvec *m0;
     268      ideal e0=singclap_factorize(mpDetBareiss(M0),&m0,2);
     269
     270      for(int i=0;i<IDELEMS(e0);i++)
     271      {
     272        number e1=nNeg(pGetCoeff(e0->m[i]));
     273        pDeleteLm(&e0->m[i]);
     274        if(pGetExp(e0->m[i],1)==0)
     275          e->m[k]=pNSet(nDiv(pGetCoeff(e0->m[i]),e1));
     276        else
     277          e->m[k]=pNSet(nDiv(e1,pGetCoeff(e0->m[i])));
     278        nDelete(&e1);
     279        (*m)[k]=(*m0)[i];
     280        k++;
     281      }
     282
     283      delete(m0);
     284      idDelete(&e0);
     285    }
     286  }
     287
     288  pDelete(&t);
     289  idDelete((ideal *)&M);
     290
     291  for(int i=0;i<n-1;i++)
     292  {
     293    if(e->m[i]!=NULL)
     294    for(int j=i+1;j<n;j++)
     295    {
     296      if(e->m[j]!=NULL)
     297      if(nEqual(pGetCoeff(e->m[i]),pGetCoeff(e->m[j])))
     298      {
     299        (*m)[i]+=(*m)[j];
     300        (*m)[j]=0;
     301      }
     302      else
     303      if(nGreater(pGetCoeff(e->m[i]),pGetCoeff(e->m[j])))
     304      {
     305        poly p=e->m[i];
     306        e->m[i]=e->m[j];
     307        e->m[j]=p;
     308        int k=(*m)[i];
     309        (*m)[i]=(*m)[j];
     310        (*m)[j]=k;
     311      }
     312    }
     313  }
     314
     315  int n0=0;
     316  for(int i=0;i<n;i++)
     317    if((*m)[i]>0)
     318      n0++;
     319
     320  ideal e0=idInit(n0,1);
     321  intvec *m0=new intvec(n0);
     322
     323  for(int i=0,i0=0;i<n;i++)
     324    if((*m)[i]>0)
     325    {
     326      e0->m[i0]=e->m[i];
     327      e->m[i]=NULL;
     328      (*m0)[i0]=(*m)[i];
     329      i0++;
     330    }
     331
     332  idDelete(&e);
     333  delete(m);
     334
    160335  l->Init(2);
    161336  l->m[0].rtyp=IDEAL_CMD;
    162   l->m[0].data=NULL;
     337  l->m[0].data=e0;
    163338  l->m[1].rtyp=INTVEC_CMD;
    164   l->m[1].data=new intvec;
    165   int j=1;
    166   while(j<=n)
    167   {
    168     while(j<n&&(MATELEM(M,j,j+1)==NULL||MATELEM(M,j+1,j)==NULL)||j==n)
    169     {
    170       l=addval(l,MATELEM(M,j,j),1);
    171       MATELEM(M,j,j)=NULL;
    172       j++;
    173     }
    174     if(j<n)
    175     {
    176       poly t=pOne();
    177       pSetExp(t,1,1);
    178       pSetm(t);
    179       poly d0=pSub(MATELEM(M,j,j),t);
    180       MATELEM(M,j,j)=NULL;
    181       j++;
    182       poly d1=pOne();
    183       poly d2=NULL;
    184       while(j<=n&&MATELEM(M,j,j-1)!=NULL&&MATELEM(M,j-1,j)!=NULL)
    185       {
    186         d2=d1;
    187         d1=pCopy(d0);
    188         poly t=pOne();
    189         pSetExp(t,1,1);
    190         pSetm(t);
    191         d0=pSub(pMult(d0,pSub(MATELEM(M,j,j),t)),
    192                 pMult(pMult(d2,MATELEM(M,j-1,j)),MATELEM(M,j,j-1)));
    193         MATELEM(M,j,j)=NULL;
    194         MATELEM(M,j-1,j)=NULL;
    195         MATELEM(M,j,j-1)=NULL;
    196         j++;
    197       }
    198       pDelete(&d1);
    199       intvec *m0=NULL;
    200 #ifdef HAVE_FACTORY
    201       ideal e0=singclap_factorize(d0,&m0,0);
    202 #else
    203       ideal e0 = NULL;
    204 #endif
    205       pDelete(&d0);
    206       for(int i=IDELEMS(e0)-1;i>=1;i--)
    207       {
    208         poly p=e0->m[i];
    209         e0->m[i]=NULL;
    210         poly p0,p1;
    211         poly pp=p;
    212         while(pp!=NULL&&pGetExp(pp,1)<=1)
    213         {
    214           if(pGetExp(pp,1)==0)
    215             p0=pp;
    216           else
    217           if(pGetExp(pp,1)==1)
    218             p1=pp;
    219           pp=pNext(pp);
    220         }
    221         if(pp==NULL)
    222         {
    223           pp=p;
    224           p=pNSet(nNeg(nDiv(pGetCoeff(p0),pGetCoeff(p1))));
    225           pDelete(&pp);
    226         }
    227         else
    228         {
    229           p=pMult_nn(p,pGetCoeff(e0->m[0]));
    230         }
    231         l=addval(l,p,(*m0)[i]);
    232       }
    233       delete m0;
    234       idDelete(&e0);
    235     }
    236   }
    237   idDelete((ideal*)&M);
    238   return sortval(l);
    239 }
    240 
    241 BOOLEAN tridiag(leftv res,leftv h)
    242 {
    243   if((h!=NULL) && (h->Typ()==MATRIX_CMD))
    244   {
    245     matrix M=(matrix)h->Data();
    246     if(MATCOLS(M)!=MATROWS(M))
    247     {
    248       WerrorS("square matrix expected");
    249     }
    250     res->rtyp=MATRIX_CMD;
    251     res->data=(void*)tridiag(mpCopy(M));
    252     return FALSE;
    253   }
    254   WerrorS("<matrix> expected");
    255   return TRUE;
    256 }
    257 
    258 BOOLEAN eigenval(leftv res,leftv h)
    259 {
    260   if((h!=NULL) && (h->Typ()==MATRIX_CMD))
    261   {
    262     matrix M=(matrix)h->Data();
    263     if(MATCOLS(M)!=MATROWS(M))
    264     {
    265       WerrorS("square matrix expected");
    266     }
    267     res->rtyp=LIST_CMD;
    268     res->data=(void*)eigenval(mpCopy(M));
    269     return FALSE;
    270   }
    271   WerrorS("<matrix> expected");
    272   return TRUE;
    273 }
    274 #endif
     339  l->m[1].data=m0;
     340
     341  return(l);
     342}
     343
     344
     345BOOLEAN evEigenvalue(leftv res,leftv h)
     346{
     347  if(currRingHdl)
     348  {
     349    if(h&&h->Typ()==MATRIX_CMD)
     350    {
     351      matrix M=(matrix)h->Data();
     352      res->rtyp=LIST_CMD;
     353      res->data=(void*)evEigenvalue(mpCopy(M));
     354      return FALSE;
     355    }
     356    WerrorS("<matrix> expected");
     357    return TRUE;
     358  }
     359  WerrorS("no ring active");
     360  return TRUE;
     361}
     362
     363#endif /* HAVE_EIGENVAL */

  • Singular/eigenval.h

    r502ea2 ra7f2a04  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: eigenval.h,v 1.1 2001-03-05 17:31:40 mschulze Exp $ */
     4/* $Id: eigenval.h,v 1.2 2002-02-16 10:56:19 mschulze Exp $ */
    55/*
    66* ABSTRACT: eigenvalues of constant square matrices
    … …  
    1010#define EIGENVAL_H
    1111
    12 matrix tridiag(matrix M);
    13 lists eigenval(matrix M);
    14 BOOLEAN tridiag(leftv res,leftv h);
    15 BOOLEAN eigenval(leftv res,leftv h);
     12matrix evSwap(matrix M,int i,int j);
     13BOOLEAN evSwap(leftv res,leftv h);
     14matrix evRowElim(matrix M,int i,int j,int k);
     15BOOLEAN evRowElim(leftv res,leftv h);
     16matrix evColElim(matrix M,int i,int j,int k);
     17BOOLEAN evColElim(leftv res,leftv h);
     18matrix evHessenberg(matrix M);
     19BOOLEAN evHessenberg(leftv res,leftv h);
     20lists evEigenvalue(matrix M);
     21BOOLEAN evEigenvalue(leftv res,leftv h);
    1622
    17 #endif
     23#endif /* EIGENVAL_H */
Note: See TracChangeset for help on using the changeset viewer.