My Project
Loading...
Searching...
No Matches
sispasm.cc
Go to the documentation of this file.
1/*
2 * provides (after defintion of a ring with coeffs in Z/p)
3 * - type spasm
4 * - assignment smatrix ->spasm, matrix ->spasm
5 * - printing/string(spasm)
6 * - transpose(spasm) -> spasm
7 * - to_matrix(spams) -> matrix
8 * - to_smatrix(spasm) -> smatrix
9 * - spasm_kernel(spasm)->spasm
10 * - spasm_rref(spasm) -> spasm
11*/
12#include "singularconfig.h"
14#include "kernel/ideals.h"
15#include "Singular/ipid.h"
16#include "Singular/ipshell.h"
17#include "Singular/blackbox.h"
18#include "Singular/mod_lib.h"
19#ifdef HAVE_SPASM_H
20extern "C"
21{
22#include "spasm.h"
23}
24
25spasm* conv_matrix2spasm(matrix M, const ring R)
26{
27 int i=MATROWS(M);
28 int j=MATCOLS(M);
29 spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
30 for (int ii=1;ii<=i;ii++)
31 {
32 for(int jj=1;jj<=j;jj++)
33 {
34 poly p;
35 if ((p=MATELEM(M,ii,jj))!=NULL)
36 {
37 spasm_add_entry(T,ii-1,jj-1,(spasm_GFp)(long)pGetCoeff(p));
38 }
39 }
40 }
41 spasm* A=spasm_compress(T);
42 spasm_triplet_free(T);
43 return A;
44}
45
46spasm* conv_smatrix2spasm(ideal M, const ring R)
47{
48 int i=MATROWS((matrix)M);
49 int j=MATCOLS((matrix)M);
50 spasm_triplet *T = spasm_triplet_alloc(i, j, 1, R->cf->ch, 1);
51 for(int jj=0;jj<j;jj++)
52 {
53 poly p=M->m[jj];
54 while (p!=NULL)
55 {
56 int ii=p_GetComp(p,R);
57 spasm_add_entry(T,ii-1,jj,(spasm_GFp)(long)pGetCoeff(p));
58 pIter(p);
59 }
60 }
61 spasm* A=spasm_compress(T);
62 spasm_triplet_free(T);
63 return A;
64}
65
66matrix conv_spasm2matrix(spasm *A, const ring R)
67{
68 matrix M=mpNew(A->n,A->m);
69 int n=A->n;
70 int *Aj = A->j;
71 int *Ap = A->p;
72 spasm_GFp *Ax = A->x;
73 for (int i = 0; i < n; i++)
74 {
75 for (int px = Ap[i]; px < Ap[i + 1]; px++)
76 {
77 spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
78 MATELEM(M,i+1,Aj[px] + 1)=p_ISet(x,R);
79 }
80 }
81 return M;
82}
83
84ideal conv_spasm2smatrix(spasm *A, const ring R)
85{
86 ideal M=idInit(A->m,A->n);
87 int n=A->n;
88 int *Aj = A->j;
89 int *Ap = A->p;
90 spasm_GFp *Ax = A->x;
91 for (int i = 0; i < n; i++)
92 {
93 for (int px = Ap[i]; px < Ap[i + 1]; px++)
94 {
95 spasm_GFp x = (Ax != NULL) ? Ax[px] : 1;
96 poly p=p_ISet(x,R);
98 M->m[Aj[px]]=p_Add_q(M->m[Aj[px]],p,R);
99 }
100 }
101 return M;
102}
103
104spasm* sp_kernel(spasm* A, const ring R)
105{ /*from kernel.c*/
106 int n = A->n;
107 int m = A->m;
108 int* p = (int*)spasm_malloc(n * sizeof(int));
109 int * qinv = (int*)spasm_malloc(m * sizeof(int));
110 spasm_find_pivots(A, p, qinv); /* this does some useless stuff, but
111 * pushes zero rows to the bottom */
112 spasm* A_clean = spasm_permute(A, p, SPASM_IDENTITY_PERMUTATION, SPASM_WITH_NUMERICAL_VALUES);
113 spasm_csr_free(A);
114 A = A_clean;
115 for (int i = 0; i < n; i++)
116 {
117 if (spasm_row_weight(A, i) == 0)
118 {
119 //fprintf(stderr, "[kernel] ignoring %d empty rows\n", n - i);
120 A->n = i;
121 n = i;
122 break;
123 }
124 }
125
126 spasm* A_t = spasm_transpose(A, SPASM_WITH_NUMERICAL_VALUES);
127 spasm_find_pivots(A_t, qinv, p);
128
129 spasm* K = spasm_kernel(A_t, qinv);
130 free(p);
131 free(qinv);
132 spasm_csr_free(A_t);
133 return K;
134}
135
136spasm* sp_rref(spasm* A)
137{ /* from rref_gplu.c: compute an echelonized form, WITHOUT COLUMN PERMUTATION */
138 spasm_lu *LU = spasm_LU(A, SPASM_IDENTITY_PERMUTATION, 1);
139 spasm *U = spasm_transpose(LU->L, 1);
140 spasm_make_pivots_unitary(U, SPASM_IDENTITY_PERMUTATION, U->n);
141 spasm_free_LU(LU);
142 return U;
143}
144
145spasm* sp_Mult_v(spasm* A, int *v)
146{
147 int *y=(int*)omAlloc0(A->n*sizeof(int));
148 spasm *AA=spasm_submatrix(A,0,A->n,0,A->m,1); /*copy A*/
149 spasm_gaxpy(AA,v,y);
150 return AA;
151}
152/*----------------------------------------------------------------*/
153VAR int SPASM_CMD;
154
155static void* sp_Init(blackbox* /*b*/)
156{
158 {
159 spasm_triplet *T = spasm_triplet_alloc(0, 0, 1, currRing->cf->ch, 1);
160 spasm* A=spasm_compress(T);
161 spasm_triplet_free(T);
162 return (void*)A;
163 }
164 else
165 {
166 WerrorS("ring with Z/p coeffs required");
167 return NULL;
168 }
169}
170static void sp_destroy(blackbox* /*b*/, void *d)
171{
172 if (d!=NULL) spasm_csr_free((spasm*)d);
173 d=NULL;
174}
175static char* sp_String(blackbox* /*b*/, void *d)
176{
177 char buf[30];
178 spasm* A=(spasm*)d;
179 sprintf(buf,"spasm matrix %dx%d",A->n,A->m);
180 return omStrDup(buf);
181}
182static void* sp_Copy(blackbox* /*b*/, void *d)
183{
184 if (d!=NULL)
185 {
186 spasm* A=(spasm*)d;
187 spasm* B=spasm_submatrix(A,0,A->n,0,A->m,1);
188 return (void*)B;
189 }
190 return NULL;
191}
192static BOOLEAN sp_Assign(leftv l, leftv r)
193{
194 spasm* A;
195 void*d=l->Data();
196 if (d!=NULL) spasm_csr_free((spasm*)d);
197 if (l->rtyp==IDHDL)
198 {
199 IDDATA((idhdl)l->data) = (char*)NULL;
200 }
201 else
202 {
203 l->data = (void*)NULL;
204 }
205
206 if (r->Typ()==l->Typ())
207 {
208 A=(spasm*)sp_Copy(NULL,r->Data());
209 }
210 else if (r->Typ()==SMATRIX_CMD)
211 {
212 A=conv_smatrix2spasm((ideal)r->Data(),currRing);
213 }
214 else if (r->Typ()==MATRIX_CMD)
215 {
216 A=conv_matrix2spasm((matrix)r->Data(),currRing);
217 }
218 else
219 return TRUE;
220
221 if (l->rtyp==IDHDL)
222 {
223 IDDATA((idhdl)l->data) = (char*)A;
224 }
225 else
226 {
227 l->data = (void*)A;
228 }
229 return FALSE;
230}
231
232static BOOLEAN to_smatrix(leftv res, leftv args)
233{
234 leftv u = args;
235 if ((u!=NULL) && (u->Typ()==SPASM_CMD))
236 {
237 res->rtyp=SMATRIX_CMD;
238 res->data=(void*)conv_spasm2smatrix((spasm*)u->Data(),currRing);
239 return FALSE;
240 }
241 return TRUE;
242}
243static BOOLEAN to_matrix(leftv res, leftv args)
244{
245 leftv u = args;
246 if ((u!=NULL) && (u->Typ()==SPASM_CMD))
247 {
248 res->rtyp=MATRIX_CMD;
249 res->data=(void*)conv_spasm2matrix((spasm*)u->Data(),currRing);
250 return FALSE;
251 }
252 return TRUE;
253}
254static BOOLEAN kernel(leftv res, leftv args)
255{
256 leftv u = args;
257 if ((u!=NULL) && (u->Typ()==SPASM_CMD))
258 {
259 res->rtyp=SPASM_CMD;
260 res->data=(void*)sp_kernel((spasm*)u->Data(),currRing);
261 return FALSE;
262 }
263 return TRUE;
264}
265static BOOLEAN rref(leftv res, leftv args)
266{
267 leftv u = args;
268 if ((u!=NULL) && (u->Typ()==SPASM_CMD))
269 {
270 res->rtyp=SPASM_CMD;
271 res->data=(void*)sp_rref((spasm*)u->Data());
272 return FALSE;
273 }
274 return TRUE;
275}
276static BOOLEAN sp_Op1(int op,leftv l, leftv r)
277{
278 if(op==TRANSPOSE_CMD)
279 {
280 l->rtyp=r->Typ();
281 l->data=(void*)spasm_transpose((spasm*)r->Data(),SPASM_WITH_NUMERICAL_VALUES);
282 return FALSE;
283 }
284 return blackboxDefaultOp1(op,l,r);
285}
286/*----------------------------------------------------------------*/
287// initialisation of the module
288extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* p)
289{
290 blackbox *b=(blackbox*)omAlloc0(sizeof(blackbox));
291 b->blackbox_destroy=sp_destroy;
292 b->blackbox_String=sp_String;
293 b->blackbox_Init=sp_Init;
294 b->blackbox_Copy=sp_Copy;
295 b->blackbox_Assign=sp_Assign;
296 b->blackbox_Op1=sp_Op1;
297 SPASM_CMD=setBlackboxStuff(b,"spasm");
298 p->iiAddCproc("spasm.so","spasm_kernel",FALSE,kernel);
299 p->iiAddCproc("spasm.so","spasm_rref",FALSE,rref);
300 p->iiAddCproc("spasm.so","to_smatrix",FALSE,to_smatrix);
301 p->iiAddCproc("spasm.so","to_matrix",FALSE,to_matrix);
302 return (MAX_TOK);
303}
304#else
305extern "C" int SI_MOD_INIT(sispasm)(SModulFunctions* psModulFunctions)
306{
307 PrintS("no spasm support\n");
308 return MAX_TOK;
309}
310#endif
int BOOLEAN
Definition: auxiliary.h:87
#define TRUE
Definition: auxiliary.h:100
#define FALSE
Definition: auxiliary.h:96
int setBlackboxStuff(blackbox *bb, const char *n)
define a new type
Definition: blackbox.cc:142
BOOLEAN blackboxDefaultOp1(int op, leftv l, leftv r)
default procedure blackboxDefaultOp1, to be called as "default:" branch
Definition: blackbox.cc:78
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4081
int p
Definition: cfModGcd.cc:4077
CanonicalForm b
Definition: cfModGcd.cc:4102
Definition: idrec.h:35
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
int Typ()
Definition: subexpr.cc:1011
void * Data()
Definition: subexpr.cc:1154
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
b *CanonicalForm B
Definition: facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define VAR
Definition: globaldefs.h:5
@ MATRIX_CMD
Definition: grammar.cc:286
@ SMATRIX_CMD
Definition: grammar.cc:291
#define IDDATA(a)
Definition: ipid.h:126
STATIC_VAR jList * T
Definition: janet.cc:30
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition: matpol.cc:37
#define MATELEM(mat, i, j)
1-based access to matrix
Definition: matpol.h:29
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
#define p_GetComp(p, r)
Definition: monomials.h:64
#define pIter(p)
Definition: monomials.h:37
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
#define omStrDup(s)
Definition: omAllocDecl.h:263
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define free
Definition: omAllocFunc.c:14
#define NULL
Definition: omList.c:12
poly p_ISet(long i, const ring r)
returns the poly representing the integer i
Definition: p_polys.cc:1297
static poly p_Add_q(poly p, poly q, const ring r)
Definition: p_polys.h:938
static unsigned long p_SetComp(poly p, unsigned long c, ring r)
Definition: p_polys.h:249
#define p_SetmComp
Definition: p_polys.h:246
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
void PrintS(const char *s)
Definition: reporter.cc:284
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:501
int status int void * buf
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
int SI_MOD_INIT() sispasm(SModulFunctions *psModulFunctions)
Definition: sispasm.cc:305
#define IDHDL
Definition: tok.h:31
@ TRANSPOSE_CMD
Definition: tok.h:191
@ MAX_TOK
Definition: tok.h:218