Changeset c33a45 in git


Ignore:
Timestamp:
Mar 2, 2017, 4:12:22 PM (7 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '648d28f488f6ff08f5607ff229b9ad9e4a5b93c2')
Children:
77ba486a2bb537dad7c328a5710ec099c240c06f
Parents:
dc29160a567b9e04c172ea21b48f03a654afd0a8
Message:
add: mp_Tensor for Kronecker product of 2 matrices given as modules
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • Singular/extra.cc

    rdc2916 rc33a45  
    38343834    }
    38353835    else
     3836/*==================== tensor =================*/
     3837    if(strcmp(sys_cmd,"tensor")==0)
     3838    {
     3839      if ((h->Typ()==MODUL_CMD)
     3840      && (h->next!=NULL)
     3841      && (h->next->Typ()==MODUL_CMD))
     3842      {
     3843        res->data=(void*)mp_Tensor((ideal)h->Data(),(ideal)h->next->Data(),currRing);
     3844        res->rtyp=MODUL_CMD;
     3845        return FALSE;
     3846      }
     3847      else
     3848        return TRUE;
     3849    }
     3850    else
    38363851/*==================== power* ==================================*/
    38373852    #if 0
  • libpolys/polys/matpol.cc

    rdc2916 rc33a45  
    99#include <stdio.h>
    1010#include <math.h>
    11 
    12 
    13 
    1411
    1512#include <misc/auxiliary.h>
     
    17101707  return (result);
    17111708}
     1709
     1710static void p_DecomposeComp(poly p, poly *a, int l, const ring r)
     1711{
     1712  poly h=p;
     1713  while(h!=NULL)
     1714  {
     1715    poly hh=pNext(h);
     1716    pNext(h)=a[p_GetComp(h,r)-1];
     1717    a[p_GetComp(h,r)-1]=h;
     1718    p_SetComp(h,0,r);
     1719    p_SetmComp(h,r);
     1720    h=hh;
     1721  }
     1722  for(int i=0;i<l;i++)
     1723  {
     1724    if(a[i]!=NULL) a[i]=pReverse(a[i]);
     1725  }
     1726}
     1727
     1728static ideal mp_MultAndShift(poly f, ideal B, int s, const ring r)
     1729{
     1730  ideal res=idInit(IDELEMS(B),B->rank+s);
     1731  if (f!=NULL)
     1732  {
     1733    int q=IDELEMS(B); // p x q
     1734    for(int j=0;j<q;j++)
     1735    {
     1736      poly h=pp_Mult_qq(f,B->m[j],r);
     1737      if (h!=NULL)
     1738      {
     1739        if (s>0) p_Shift(&h,s,r);
     1740        res->m[j]=h;
     1741      }
     1742    }
     1743  }
     1744  return res;
     1745}
     1746static ideal mp_AddSubMat(ideal res, ideal sm, int col, const ring r)
     1747{
     1748  for(int i=0;i<IDELEMS(sm);i++)
     1749  {
     1750    res->m[col+i]=p_Add_q(res->m[col+i],sm->m[i],r);
     1751    sm->m[i]=NULL;
     1752  }
     1753  return res;
     1754}
     1755
     1756ideal mp_Tensor(ideal A, ideal B, const ring r)
     1757{
     1758  // size of the result n*q x m*p
     1759  int n=IDELEMS(A); // m x n
     1760  int m=A->rank;
     1761  int q=IDELEMS(B); // p x q
     1762  int p=B->rank;
     1763  ideal res=idInit(n*q,m*p);
     1764  poly *a=(poly*)omAlloc0(m*sizeof(poly));
     1765  for(int i=0; i<n; i++)
     1766  {
     1767    p_DecomposeComp(p_Copy(A->m[i],r),a,m,r);
     1768    for(int j=0;j<m;j++)
     1769    {
     1770      ideal sm=mp_MultAndShift(a[j], // A_i_j
     1771                               B,
     1772                               j*p, // shift j*p down
     1773                               r);
     1774      res=mp_AddSubMat(res,sm,i*q,r); // add this colums to col i*q ff
     1775      id_Delete(&sm,r);
     1776    }
     1777    memset(a,0,m*sizeof(poly));
     1778  }
     1779  omFreeSize(a,m*sizeof(poly));
     1780  return res;
     1781}
     1782
  • libpolys/polys/matpol.h

    rdc2916 rc33a45  
    9494
    9595int mp_Compare(matrix a, matrix b, const ring r);
     96
     97ideal mp_Tensor(ideal A, ideal B, const ring r);
    9698#endif/* MATPOL_H */
Note: See TracChangeset for help on using the changeset viewer.