Changeset c33a45 in git

Timestamp:

Mar 2, 2017, 4:12:22 PM (7 years ago)

Author:

Hans Schoenemann <hannes@…>

Branches:

(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')

Children:

77ba486a2bb537dad7c328a5710ec099c240c06f

Parents:

dc29160a567b9e04c172ea21b48f03a654afd0a8

Message:

add: mp_Tensor for Kronecker product of 2 matrices given as modules

Files:

• Singular/extra.cc (modified) (1 diff)
• libpolys/polys/matpol.cc (modified) (2 diffs)
• libpolys/polys/matpol.h (modified) (1 diff)

Singular/extra.cc

@@ -3834 +3834 @@
     }
     else
+/*==================== tensor =================*/
+    if(strcmp(sys_cmd,"tensor")==0)
+    {
+      if ((h->Typ()==MODUL_CMD)
+      && (h->next!=NULL)
+      && (h->next->Typ()==MODUL_CMD))
+      {
+        res->data=(void*)mp_Tensor((ideal)h->Data(),(ideal)h->next->Data(),currRing);
+        res->rtyp=MODUL_CMD;
+        return FALSE;
+      }
+      else
+        return TRUE;
+    }
+    else
 /*==================== power* ==================================*/
     #if 0

libpolys/polys/matpol.cc

@@ -9 +9 @@
 #include <stdio.h>
 #include <math.h>
-
-
-
 
 #include <misc/auxiliary.h>
@@ -1710 +1707 @@
   return (result);
 }
+
+static void p_DecomposeComp(poly p, poly *a, int l, const ring r)
+{
+  poly h=p;
+  while(h!=NULL)
+  {
+    poly hh=pNext(h);
+    pNext(h)=a[p_GetComp(h,r)-1];
+    a[p_GetComp(h,r)-1]=h;
+    p_SetComp(h,0,r);
+    p_SetmComp(h,r);
+    h=hh;
+  }
+  for(int i=0;i<l;i++)
+  {
+    if(a[i]!=NULL) a[i]=pReverse(a[i]);
+  }
+}
+
+static ideal mp_MultAndShift(poly f, ideal B, int s, const ring r)
+{
+  ideal res=idInit(IDELEMS(B),B->rank+s);
+  if (f!=NULL)
+  {
+    int q=IDELEMS(B); // p x q
+    for(int j=0;j<q;j++)
+    {
+      poly h=pp_Mult_qq(f,B->m[j],r);
+      if (h!=NULL)
+      {
+        if (s>0) p_Shift(&h,s,r);
+        res->m[j]=h;
+      }
+    }
+  }
+  return res;
+}
+static ideal mp_AddSubMat(ideal res, ideal sm, int col, const ring r)
+{
+  for(int i=0;i<IDELEMS(sm);i++)
+  {
+    res->m[col+i]=p_Add_q(res->m[col+i],sm->m[i],r);
+    sm->m[i]=NULL;
+  }
+  return res;
+}
+
+ideal mp_Tensor(ideal A, ideal B, const ring r)
+{
+  // size of the result n*q x m*p
+  int n=IDELEMS(A); // m x n
+  int m=A->rank;
+  int q=IDELEMS(B); // p x q
+  int p=B->rank;
+  ideal res=idInit(n*q,m*p);
+  poly *a=(poly*)omAlloc0(m*sizeof(poly));
+  for(int i=0; i<n; i++)
+  {
+    p_DecomposeComp(p_Copy(A->m[i],r),a,m,r);
+    for(int j=0;j<m;j++)
+    {
+      ideal sm=mp_MultAndShift(a[j], // A_i_j
+                               B,
+                               j*p, // shift j*p down
+                               r);
+      res=mp_AddSubMat(res,sm,i*q,r); // add this colums to col i*q ff
+      id_Delete(&sm,r);
+    }
+    memset(a,0,m*sizeof(poly));
+  }
+  omFreeSize(a,m*sizeof(poly));
+  return res;
+}
+

libpolys/polys/matpol.h

@@ -94 +94 @@
 
 int mp_Compare(matrix a, matrix b, const ring r);
+
+ideal mp_Tensor(ideal A, ideal B, const ring r);
 #endif/* MATPOL_H */