Loading...
Searching...
No Matches
sispasm.cc
Go to the documentation of this file.
/*
* provides (after defintion of a ring with coeffs in Z/p)
* - type spasm
* - assignment smatrix,module ->spasm, matrix ->spasm
* - printing/string(spasm)
* - transpose(spasm) -> spasm
* - nrows(spasm) -> int
* - ncols(spasm) -> int
* - to_matrix(spams) -> matrix
* - to_smatrix(spasm) -> smatrix
* - spasm_kernel(spasm)->spasm
* - spasm_rref(spasm) -> spasm
* - <spasm>[<int>,<int>] -> number: reading an entry (get_spasm_entry)
*/
#include "singularconfig.h"
#include "libpolys/polys/monomials/monomials.h"
#include "kernel/ideals.h"
#include "Singular/ipid.h"
#include "Singular/ipshell.h"
#include "Singular/blackbox.h"
#include "Singular/mod_lib.h"
#ifdef HAVE_SPASM_H
extern "C"
{
#include "spasm.h"
}
static spasm* conv_matrix2spasm(matrix M, const ring R)
{
int i=MATROWS(M);
int j=MATCOLS(M);
spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
for (int ii=1;ii<=i;ii++)
{
for(int jj=1;jj<=j;jj++)
{
poly p;
if ((p=MATELEM(M,ii,jj))!=NULL)
{
spasm_add_entry(T,ii-1,jj-1,(spasm_GFp)(long)pGetCoeff(p));
}
}
}
spasm* A=spasm_compress(T);
spasm_triplet_free(T);
return A;
}
static spasm* conv_smatrix2spasm(ideal M, const ring R)
{
int i=MATROWS((matrix)M);
int j=MATCOLS((matrix)M);
spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
for(int jj=0;jj<j;jj++)
{
poly p=M->m[jj];
while (p!=NULL)
{
int ii=p_GetComp(p,R);
spasm_add_entry(T,ii-1,jj,(spasm_GFp)(long)pGetCoeff(p));
pIter(p);
}
}
spasm* A=spasm_compress(T);
spasm_triplet_free(T);
return A;
}
static matrix conv_spasm2matrix(spasm *A, const ring R)
{
matrix M=mpNew(A->n,A->m);
int n=A->n;
int *Aj = A->j;
int *Ap = A->p;
spasm_GFp *Ax = A->x;
for (int i = 0; i < n; i++)
{
for (int px = Ap[i]; px < Ap[i + 1]; px++)
{
spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
MATELEM(M,i+1,Aj[px] + 1)=p_ISet(x,R);
}
}
return M;
}
static ideal conv_spasm2smatrix(spasm *A, const ring R)
{
ideal M=idInit(A->m,A->n);
int n=A->n;
int *Aj = A->j;
int *Ap = A->p;
spasm_GFp *Ax = A->x;
for (int i = 0; i < n; i++)
{
for (int px = Ap[i]; px < Ap[i + 1]; px++)
{
spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
poly p=p_ISet(x,R);
p_SetComp(p,i+1,R);p_SetmComp(p,R);
M->m[Aj[px]]=p_Add_q(M->m[Aj[px]],p,R);
}
}
return M;
}
static number get_spasm_entry(spasm *A, int i, int j, const ring R)
{
matrix M=mpNew(A->n,A->m);
int n=A->n;
int *Aj = A->j;
int *Ap = A->p;
i--;j--;
spasm_GFp *Ax = A->x;
if (i<n)
{
for (int px = Ap[i]; px < Ap[i + 1]; px++)
{
spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
if (j==Aj[px])
return n_Init(x,R->cf);
}
}
return n_Init(0,R->cf);;
}
static spasm* sp_kernel(spasm* A, const ring R)
{
int n = A->n;
int m = A->m;
int* p = (int*)spasm_malloc(n * sizeof(int));
int * qinv = (int*)spasm_malloc(m * sizeof(int));
spasm_find_pivots(A, p, qinv); /* this does some useless stuff, but
* pushes zero rows to the bottom */
#if 0
/*from kernel.c*/
spasm* A_clean = spasm_permute(A, p, SPASM_IDENTITY_PERMUTATION, SPASM_WITH_NUMERICAL_VALUES);
A = A_clean;
for (int i = 0; i < n; i++)
{
if (spasm_row_weight(A, i) == 0)
{
//fprintf(stderr, "[kernel] ignoring %d empty rows\n", n - i);
A->n = i;
n = i;
break;
}
}
spasm* A_t = spasm_transpose(A, SPASM_WITH_NUMERICAL_VALUES);
spasm_find_pivots(A_t, qinv, p);
spasm* K = spasm_kernel(A_t, qinv);
spasm_csr_free(A_t);
#else
spasm* K = spasm_kernel(A, p);
#endif
free(p);
free(qinv);
return K;
}
static spasm* sp_rref(spasm* A)
{ /* from rref_gplu.c: compute an echelonized form, WITHOUT COLUMN PERMUTATION */
spasm_lu *LU = spasm_LU(A, SPASM_IDENTITY_PERMUTATION, 1);
spasm *U = spasm_transpose(LU->L, 1);
spasm_make_pivots_unitary(U, SPASM_IDENTITY_PERMUTATION, U->n);
spasm_free_LU(LU);
return U;
}
static spasm* sp_Mult_v(spasm* A, int *v)
{
int *y=(int*)omAlloc0(A->n*sizeof(int));
spasm *AA=spasm_submatrix(A,0,A->n,0,A->m,1); /*copy A*/
spasm_gaxpy(AA,v,y);
return AA;
}
/*----------------------------------------------------------------*/
VAR int SPASM_CMD;
static void* sp_Init(blackbox* /*b*/)
{
if ((currRing!=NULL)&&(rField_is_Zp(currRing)))
{
spasm_triplet *T = spasm_triplet_alloc(0, 0, 1, currRing->cf->ch, 1);
spasm* A=spasm_compress(T);
spasm_triplet_free(T);
return (void*)A;
}
else
{
WerrorS("ring with Z/p coeffs required");
return NULL;
}
}
static void sp_destroy(blackbox* /*b*/, void *d)
{
if (d!=NULL) spasm_csr_free((spasm*)d);
d=NULL;
}
static char* sp_String(blackbox* /*b*/, void *d)
{
char buf[30];
spasm* A=(spasm*)d;
snprintf(buf,30,"spasm matrix %dx%d",A->n,A->m);
return omStrDup(buf);
}
static void* sp_Copy(blackbox* /*b*/, void *d)
{
if (d!=NULL)
{
spasm* A=(spasm*)d;
spasm* B=spasm_submatrix(A,0,A->n,0,A->m,1);
return (void*)B;
}
return NULL;
}
static BOOLEAN sp_Assign(leftv l, leftv r)
{
spasm* A;
void*d=l->Data();
if (d!=NULL) spasm_csr_free((spasm*)d);
if (l->rtyp==IDHDL)
{
IDDATA((idhdl)l->data) = (char*)NULL;
}
else
{
l->data = (void*)NULL;
}
int rt=r->Typ();
if (rt==l->Typ())
{
A=(spasm*)sp_Copy(NULL,r->Data());
}
else if ((rt==SMATRIX_CMD)||(rt==MODUL_CMD))
{
A=conv_smatrix2spasm((ideal)r->Data(),currRing);
}
else if (rt==MATRIX_CMD)
{
A=conv_matrix2spasm((matrix)r->Data(),currRing);
}
else
return TRUE;
if (l->rtyp==IDHDL)
{
IDDATA((idhdl)l->data) = (char*)A;
}
else
{
l->data = (void*)A;
}
return FALSE;
}
static BOOLEAN to_smatrix(leftv res, leftv args)
{
leftv u = args;
if ((u!=NULL) && (u->Typ()==SPASM_CMD))
{
res->rtyp=SMATRIX_CMD;
res->data=(void*)conv_spasm2smatrix((spasm*)u->Data(),currRing);
return FALSE;
}
return TRUE;
}
static BOOLEAN to_matrix(leftv res, leftv args)
{
leftv u = args;
if ((u!=NULL) && (u->Typ()==SPASM_CMD))
{
res->rtyp=MATRIX_CMD;
res->data=(void*)conv_spasm2matrix((spasm*)u->Data(),currRing);
return FALSE;
}
return TRUE;
}
static BOOLEAN kernel(leftv res, leftv args)
{
leftv u = args;
if ((u!=NULL) && (u->Typ()==SPASM_CMD))
{
res->rtyp=SPASM_CMD;
res->data=(void*)sp_kernel((spasm*)u->Data(),currRing);
return FALSE;
}
return TRUE;
}
static BOOLEAN rref(leftv res, leftv args)
{
leftv u = args;
if ((u!=NULL) && (u->Typ()==SPASM_CMD))
{
res->rtyp=SPASM_CMD;
res->data=(void*)sp_rref((spasm*)u->Data());
return FALSE;
}
return TRUE;
}
static BOOLEAN sp_Op1(int op,leftv res, leftv arg)
{
if(op==TRANSPOSE_CMD)
{
res->rtyp=arg->Typ();
res->data=(void*)spasm_transpose((spasm*)arg->Data(),SPASM_WITH_NUMERICAL_VALUES);
return FALSE;
}
else if (op==COLS_CMD)
{
spasm* A=(spasm*)arg->Data();
res->rtyp=INT_CMD;
res->data=(void*)(long)A->m;
return FALSE;
}
else if (op==ROWS_CMD)
{
spasm* A=(spasm*)arg->Data();
res->rtyp=INT_CMD;
res->data=(void*)(long)A->n;
return FALSE;
}
return blackboxDefaultOp1(op,res,arg);
}
static BOOLEAN sp_Op3(int op,leftv res, leftv a1, leftv a2, leftv a3)
{
if ((op=='[')
&& (a2->Typ()==INT_CMD)
&& (a3->Typ()==INT_CMD))
{
res->rtyp=NUMBER_CMD;
res->data=(char*)get_spasm_entry((spasm*)a1->Data(),(int)(long)a2->Data(),
(int)(long)a3->Data(),currRing);
return FALSE;
}
return blackboxDefaultOp3(op,res,a1,a2,a3);
}
/*----------------------------------------------------------------*/
// initialisation of the module
extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* p)
{
blackbox *b=(blackbox*)omAlloc0(sizeof(blackbox));
b->blackbox_destroy=sp_destroy;
b->blackbox_String=sp_String;
b->blackbox_Init=sp_Init;
b->blackbox_Copy=sp_Copy;
b->blackbox_Assign=sp_Assign;
b->blackbox_Op1=sp_Op1;
b->blackbox_Op3=sp_Op3;
SPASM_CMD=setBlackboxStuff(b,"spasm");
p->iiAddCproc("spasm.so","spasm_kernel",FALSE,kernel);
p->iiAddCproc("spasm.so","spasm_rref",FALSE,rref);
p->iiAddCproc("spasm.so","to_smatrix",FALSE,to_smatrix);
p->iiAddCproc("spasm.so","to_matrix",FALSE,to_matrix);
return (MAX_TOK);
}
#else
extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* psModulFunctions)
{
PrintS("no spasm support\n");
return MAX_TOK;
}
#endif
BOOLEAN
int BOOLEAN
Definition auxiliary.h:88
TRUE
#define TRUE
Definition auxiliary.h:101
FALSE
#define FALSE
Definition auxiliary.h:97
setBlackboxStuff
int setBlackboxStuff(blackbox *bb, const char *n)
define a new type
Definition blackbox.cc:143
blackboxDefaultOp3
BOOLEAN blackboxDefaultOp3(int, leftv, leftv, leftv, leftv)
default procedure blackboxDefaultOp3, to be called as "default:" branch
Definition blackbox.cc:102
blackboxDefaultOp1
BOOLEAN blackboxDefaultOp1(int op, leftv l, leftv r)
default procedure blackboxDefaultOp1, to be called as "default:" branch
Definition blackbox.cc:78
blackbox.h
l
int l
Definition cfEzgcd.cc:100
m
int m
Definition cfEzgcd.cc:128
i
int i
Definition cfEzgcd.cc:132
x
Variable x
Definition cfModGcd.cc:4090
p
int p
Definition cfModGcd.cc:4086
b
CanonicalForm b
Definition cfModGcd.cc:4111
sleftv::Typ
int Typ()
Definition subexpr.cc:1048
sleftv::Data
void * Data()
Definition subexpr.cc:1192
n_Init
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition coeffs.h:541
y
const CanonicalForm int const CFList const Variable & y
Definition facAbsFact.cc:53
res
CanonicalForm res
Definition facAbsFact.cc:60
B
b *CanonicalForm B
Definition facBivar.cc:52
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
j
int j
Definition facHensel.cc:110
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
VAR
#define VAR
Definition globaldefs.h:5
MATRIX_CMD
@ MATRIX_CMD
Definition grammar.cc:287
MODUL_CMD
@ MODUL_CMD
Definition grammar.cc:288
SMATRIX_CMD
@ SMATRIX_CMD
Definition grammar.cc:292
NUMBER_CMD
@ NUMBER_CMD
Definition grammar.cc:289
ideals.h
ipid.h
IDDATA
#define IDDATA(a)
Definition ipid.h:126
SModulFunctions
Definition ipid.h:69
ipshell.h
T
STATIC_VAR jList * T
Definition janet.cc:30
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition matpol.cc:37
MATELEM
#define MATELEM(mat, i, j)
1-based access to matrix
Definition matpol.h:29
matrix
ip_smatrix * matrix
Definition matpol.h:43
MATROWS
#define MATROWS(i)
Definition matpol.h:26
MATCOLS
#define MATCOLS(i)
Definition matpol.h:27
mod_lib.h
monomials.h
p_GetComp
#define p_GetComp(p, r)
Definition monomials.h:64
pIter
#define pIter(p)
Definition monomials.h:37
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition monomials.h:44
omStrDup
#define omStrDup(s)
Definition omAllocDecl.h:263
omAlloc0
#define omAlloc0(size)
Definition omAllocDecl.h:211
free
#define free
Definition omAllocFunc.c:14
NULL
#define NULL
Definition omList.c:12
p_ISet
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition p_polys.cc:1298
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition p_polys.h:938
p_SetComp
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition p_polys.h:249
p_SetmComp
#define p_SetmComp
Definition p_polys.h:246
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
PrintS
void PrintS(const char *s)
Definition reporter.cc:288
rField_is_Zp
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
idhdl
idrec * idhdl
Definition ring.h:22
buf
int status int void * buf
Definition si_signals.h:69
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition simpleideals.cc:35
R
#define R
Definition sirandom.c:27
A
#define A
Definition sirandom.c:24
M
#define M
Definition sirandom.c:25
sispasm
int SI_MOD_INIT sispasm(SModulFunctions *psModulFunctions)
Definition sispasm.cc:361
leftv
sleftv * leftv
Definition structs.h:53
IDHDL
#define IDHDL
Definition tok.h:31
TRANSPOSE_CMD
@ TRANSPOSE_CMD
Definition tok.h:193
INT_CMD
@ INT_CMD
Definition tok.h:96
ROWS_CMD
@ ROWS_CMD
Definition tok.h:177
MAX_TOK
@ MAX_TOK
Definition tok.h:220
COLS_CMD
@ COLS_CMD
Definition tok.h:52