1 /************************************************************/
2 /*! \file matrix.cc
3 * \brief C++ function to add matrix support
4 *
5 * \version 1
6 *
7 * \date Created: 19/02/06
8 * \date Last modified: Time-stamp: <2006-06-17 23:10:44 antoine>
9 *
10 * \author A. Lucas
11 *
12 * \note
13 * as usually, matrix x[i,j] (n x p) is represented by a vector
14 * x[i + j x n] (i=0..n-1 ; j=0..p-1)
15 *
16 * \note Licence: GPL (>= 2)
17 */
18
19 #include <vector>
20
21 using namespace std;
22
23 #include "Rgmp.h"
24
25 #include "biginteger.h"
26 #include "bigintegerR.h"
27 #include "matrix.h"
28 // need to call matrix_mul_q()
29 #include "matrixq.h"
30
31 // given that x is "bigz" or "bigq",
32 // return TRUE if x is a bigz/q *matrix*: R's is.matrixZQ(.)
is_matrix_zq(SEXP x)33 SEXP is_matrix_zq(SEXP x) {
34 SEXP nrowSexp = Rf_mkString("nrow");
35 PROTECT(nrowSexp);
36 SEXP attributeRow = Rf_getAttrib(x,nrowSexp );
37 PROTECT(attributeRow);
38 SEXP ans = Rf_ScalarLogical(attributeRow != R_NilValue);
39 UNPROTECT(2);
40 return ans;
41 }
42
43 // C++ side of R function matrix.bigz()
as_matrixz(SEXP x,SEXP nrR,SEXP ncR,SEXP byrowR,SEXP mod)44 SEXP as_matrixz (SEXP x, SEXP nrR, SEXP ncR, SEXP byrowR, SEXP mod)
45 {
46 int
47 nc=INTEGER(ncR)[0],
48 nr=INTEGER(nrR)[0],
49 byrow=INTEGER(byrowR)[0];
50
51 // get "bigz" vector, this makes all conversion int to bigz etc...
52 bigvec mat = bigintegerR::create_bignum(x);
53 int lendat = mat.value.size();
54 // int sizemod = mat.modulus.size();
55 // when modulus specified
56 bigvec modulus = bigintegerR::create_bignum(mod);
57 if(modulus.value.size()>0) // should be allways the case
58 if(!modulus.value[0].isNA()) {
59 mat.modulus.resize(modulus.size());
60 for (unsigned int i = 0; i < modulus.size(); i++)
61 mat.modulus[i] = modulus.value[i];
62 // sizemod = modulus.size();
63 }
64
65 // A copy from R project to get correct dimension
66 // all warnings...
67 //
68 if (nr == NA_INTEGER) // This is < 0
69 error(_("matrix: invalid 'nrow' value (too large or NA)"));
70 if (nr < 0)
71 error(_("matrix: invalid 'nrow' value (< 0)"));
72 if (nc < 0)
73 error(_("matrix: invalid 'ncol' value (< 0)"));
74 if (nc == NA_INTEGER)
75 error(_("matrix: invalid 'ncol' value (too large or NA)"));
76
77 if(lendat > 0 ) {
78 if (lendat > 1 && (nr * nc) % lendat != 0) {
79 if (((lendat > nr) && (lendat / nr) * nr != lendat) ||
80 ((lendat < nr) && (nr / lendat) * lendat != nr))
81 warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d] in matrix"),
82 lendat, nr);
83 else if (((lendat > nc) && (lendat / nc) * nc != lendat) ||
84 ((lendat < nc) && (nc / lendat) * lendat != nc))
85 warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d] in matrix"),
86 lendat, nc);
87 }
88 else if ((lendat > 1) && (nr * nc == 0)){
89 warning(_("data length exceeds size of matrix"));
90 }
91 }
92
93 // update dimension parameters
94 if(nr == 1)
95 nr = (int)ceil(lendat / (double) nc);
96 if(nc == 1)
97 nc = (int)ceil(lendat / (double)nr);
98
99 // when we extend "x"
100 if(nc*nr > lendat)
101 {
102 mat.value.resize(nr*nc);
103 for(int i = lendat; i < nr*nc; i++)
104 mat.value[i] = mat.value[i % lendat];
105 }
106 mat.nrow = nr;
107 if(byrow)
108 {
109 bigvec mat2 = matrixz::bigint_transpose (mat, nc,nr);
110 mat2.nrow = nr;// FIXME - needed ??
111 return( bigintegerR::create_SEXP (mat2));
112 }
113
114 return( bigintegerR::create_SEXP (mat));
115 }
116
117
118 /*
119 * Transposition
120 */
bigint_transposeR(SEXP x)121 SEXP bigint_transposeR(SEXP x)
122 {
123 SEXP dimKey =Rf_mkString("nrow");
124 PROTECT(dimKey);
125 SEXP dimAttr = Rf_getAttrib(x,dimKey );
126 PROTECT(dimAttr);
127 bigvec mat = bigintegerR::create_bignum(x);
128 int nr, n = mat.size();
129
130 if (dimAttr == R_NilValue) { // vector
131 nr = n;
132 } else if (TYPEOF(dimAttr) == INTSXP) {
133 nr = INTEGER(dimAttr)[0];
134 } else { nr = -1;// -Wall
135 error(_("argument must be a matrix of class \"bigz\""));
136 }
137 UNPROTECT(2);
138 int nc = (int) n / nr;
139 // Rprintf(" o bigI_tr(<%d x %d>) ..\n", nr,nc);
140 return( bigintegerR::create_SEXP(matrixz::bigint_transpose(mat, nr,nc)));
141 }
142
143
144 /* \brief matrix cross product
145 *
146 * \param a matrix (n x p)
147 * \param trans if(trans), compute tcrossprod(), else crossprod()
148 * \return crossprod(a) := t(a) %*% a [p x p] or
149 * tcrossprod(a) := a %*% t(a) [n x n]
150 */
matrix_crossp_z(SEXP a,SEXP trans)151 SEXP matrix_crossp_z (SEXP a, SEXP trans)
152 {
153 bool useMod = FALSE,
154 tr = (bool)Rf_asLogical(trans);
155 bigvec mat_a = bigintegerR::create_bignum(a);
156 int sizemod = mat_a.modulus.size(),
157 a_nrow = mat_a.nrow,
158 a_len = mat_a.size();
159
160 // in case of a vector; crossprod() returns scalar product,
161 // whereas tcrossprod() gives n x n matrix.
162 if(a_nrow < 0)
163 a_nrow = a_len;
164 int a_ncol = a_len / a_nrow;
165
166 // Result R is R[1..m, 1..m] -- and R_{ij} = sum_{k=1}^p A.. B..
167 int m, p;
168 if(tr) { // tcrossprod()
169 m= a_nrow; p= a_ncol;
170 } else { // crossprod()
171 m= a_ncol; p= a_nrow;
172 }
173 bigvec res(m*m);
174 res.nrow= m;
175
176 mpz_t R_ij, tt;
177 mpz_init(R_ij);
178 mpz_init(tt);
179 mpz_t common_modulus; mpz_init(common_modulus);
180
181 if(sizemod <= 1) { // maybe 'useMod' i.e., can use common modulus:
182 if(sizemod == 1) {
183 mpz_set(common_modulus, mat_a.modulus[0].getValueTemp());
184 if(!mat_a.modulus[0].isNA())
185 useMod = TRUE;
186 }
187 }
188
189 // here the computation:
190 for(int i=0; i < m; i++)
191 for(int j=0; j < m; j++) {
192 mpz_set_ui(R_ij, 0);
193 bool isna = false;
194 #define K_LOOP \
195 for(int k=0; k < p; k++) { \
196 /* R_ij = \sum_{k=1}^p a_{ik} b_{kj} */ \
197 if( !(A_I_K.isNA() || B_K_J.isNA())) { \
198 mpz_mul(tt, A_I_K.getValueTemp(), B_K_J.getValueTemp()); \
199 mpz_add(R_ij, tt,R_ij); \
200 } \
201 else { \
202 isna = true; break; \
203 } \
204 }
205
206 if(tr) {//------------- tcrossprod ---------------------------
207
208 #define A_I_K mat_a.value [ i + k *a_nrow]
209 #define B_K_J mat_a.value [ j + k *a_nrow]
210 K_LOOP
211 #undef A_I_K
212 #undef B_K_J
213
214 } else {//------------- crossprod ---------------------------
215
216 #define A_I_K mat_a.value [ k + i *a_nrow]
217 #define B_K_J mat_a.value [ k + j *a_nrow]
218 K_LOOP
219 #undef A_I_K
220 #undef B_K_J
221
222 }
223
224 if(isna) {
225 res.value[i + j*m].setValue(0);
226 res.value[i + j*m].NA(true);
227 }
228 else
229 res.value[i + j*m].setValue(R_ij);
230
231 } // for(i ..) for(j ..)
232
233 if(useMod)
234 res.modulus.push_back(biginteger(common_modulus));
235
236 mpz_clear(R_ij);
237 mpz_clear(tt);
238 mpz_clear(common_modulus);
239
240 return( bigintegerR::create_SEXP (res));
241 } // matrix_crossp_z()
242 #undef K_LOOP
243
244 /* \brief matrix multiplication
245 *
246 * returns matrix multiplication T(a) %*% b or b %*% T(a)
247 * \param a matrix
248 * \param b matrix
249 * \param op operation code: 0: %*%, 1: crossprod, 2: tcrossprod
250 * (same codes as in R's do_matprod() in src/main/array.c )
251 */
matrix_mul_z(SEXP a,SEXP b,SEXP op)252 SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op)
253 {
254 if(!strcmp(class_P(b), "bigq")) { // b "bigq", --> use q arithm:
255 return(matrix_mul_q(a, b, op));
256 }
257 // FIXME: we may know that a is 'bigz' - but we don't know at all about b !!
258 // ----- create_bignum(.) should be much more careful (better: have a careful option!)
259 bool useMod = FALSE;// if(useMod) use a *common* modulus
260 int o_ = Rf_asInteger(op); // INTEGER(op)[0]
261 bigvec mat_a = bigintegerR::create_bignum(a),
262 mat_b = bigintegerR::create_bignum(b);
263
264 int sizemod_a = mat_a.modulus.size(),
265 sizemod_b = mat_b.modulus.size();
266
267 int
268 a_nrow = mat_a.nrow, a_len = mat_a.size(),
269 b_nrow = mat_b.nrow, b_len = mat_b.size(),
270 a_ncol = -1, b_ncol = -1;// -Wall
271
272 // distinguish cases of vectors / matrices ---------------------
273 if(a_nrow < 0) {
274 if(b_nrow < 0) { // *both* are vectors
275 if(o_ == 0) {
276 a_nrow = 1;
277 a_ncol = a_len;
278 } else {
279 a_nrow = a_len;
280 a_ncol = 1;
281 }
282 b_nrow = b_len;
283 b_ncol = 1;
284
285 } else { // a : vector, b : matrix
286 b_ncol = b_len / b_nrow;
287 if(o_ == 0) {
288 if (a_len == b_nrow) { /* x as row vector */
289 a_nrow = 1;
290 a_ncol = b_nrow; /* == a_len */
291 }
292 else if (b_nrow == 1) { /* x as col vector */
293 a_nrow = a_len;
294 a_ncol = 1;
295 }
296 } else if(o_ == 1) { /* crossprod() */
297 if (a_len == b_nrow) { /* x is a col vector */
298 a_nrow = b_nrow; /* == a_len */
299 a_ncol = 1;
300 }
301 /* else if (b_nrow == 1) ... not being too tolerant
302 to treat x as row vector, as t(x) *is* row vector */
303 } else { // o_ == 2 -- tcrossprod()
304 if (a_len == b_ncol) { /* x as row vector */
305 a_nrow = 1;
306 a_ncol = b_ncol; /* == a_len */
307 }
308 else if (b_ncol == 1) { /* x as col vector */
309 a_nrow = a_len;
310 a_ncol = 1;
311 }
312 }
313 }
314 }
315 else if (b_nrow < 0) { // a : matrix, b : vector
316 a_ncol = a_len / a_nrow;
317 if (o_ == 0) {
318 if (b_len == a_ncol) { /* y as col vector */
319 b_nrow = a_ncol;
320 b_ncol = 1;
321 }
322 else if (a_ncol == 1) { /* y as row vector */
323 b_nrow = 1;
324 b_ncol = b_len;
325 }
326 }
327 else if (o_ == 1) { /* crossprod() */
328 if (b_len == a_nrow) { /* y is a col vector */
329 b_nrow = a_nrow;
330 b_ncol = 1;
331 }
332 }
333 else { /* tcrossprod -- y is a col vector */
334 b_nrow = b_len;
335 b_ncol = 1;
336 }
337
338 } else { // a, b *both* matrices
339 a_ncol = a_len / a_nrow;
340 b_ncol = b_len / b_nrow;
341 }
342
343 if(((o_ == 0) && (a_ncol != b_nrow)) ||
344 ((o_ == 1) && (a_nrow != b_nrow)) || // crossprod()
345 ((o_ == 2) && (a_ncol != b_ncol)) // tcrossprod()
346 )
347 error(_("Matrix dimensions do not match"));
348
349 // Result R is R[1..n, 1..m] -- and R_{ij} = sum_{k=1} ^ p A.. B..
350 int n,m, p;
351 if(o_ == 0) {
352 n= a_nrow; m= b_ncol; p= a_ncol;// = b_nrow
353 }else if (o_ == 1) {
354 n= a_ncol; m= b_ncol; p= a_nrow;// = b_nrow
355 }else if (o_ == 2) {
356 n= a_nrow; m= b_nrow; p= a_ncol;// = b_ncol
357 }else {
358 error(_("invalid 'op' code in matrix_mul_z()"));
359 n = m = p = -1;// -Wall
360 }
361
362 bigvec res(n*m);
363 res.nrow=n;
364
365 mpz_t common_modulus, tt;
366 mpz_init(tt);
367 mpz_init(common_modulus);
368
369 /* modulus when modulus are "global" (i.e. of size 1) and
370 * either are the same, or only one of a or b is specified
371 */
372 if( !(sizemod_a > 1 || sizemod_b > 1)) {
373 if((sizemod_a == 1) && (sizemod_b == 1)) {
374 mpz_set(common_modulus, mat_a.modulus[0].getValueTemp());
375 if(mpz_cmp(common_modulus, mat_b.modulus[0].getValueTemp()) == 0
376 && !mat_a.modulus[0].isNA())
377 useMod = TRUE;
378 }
379 else { // at least one of the sizemod_* is > 1 :
380 if ((sizemod_a == 1) && !mat_a.modulus[0].isNA()) {
381 mpz_set(common_modulus, mat_a[0].getModulus().getValueTemp());
382 useMod = TRUE;
383 } else if ((sizemod_b == 1) && !mat_b.modulus[0].isNA()) {
384 mpz_set(common_modulus, mat_b[0].getModulus().getValueTemp());
385 useMod = TRUE;
386 }
387 }
388 }
389 // bigmod tmp;
390
391 // here the computation:
392 for(int i=0; i < n; i++)
393 for(int j=0; j < m; j++)
394 {
395 #define R_IJ res.value[ i + j*n]
396 #define K_LOOP \
397 for(int k=0; k < p; k++) \
398 { \
399 if(A_I_K.isNA() || B_K_J.isNA()) { \
400 R_IJ.setValue(0); R_IJ.NA(true); \
401 break; \
402 } \
403 /* Z = A_I_K * B_K_J */ \
404 mpz_mul(tt, A_I_K.getValueTemp(), B_K_J.getValueTemp()); \
405 /* R_IJ = R_IJ + A_I_K * B_K_J */ \
406 mpz_add(tt, tt, R_IJ.getValueTemp()); \
407 if(useMod) \
408 mpz_mod(tt,tt,common_modulus); \
409 R_IJ.setValue(tt); \
410 }
411
412 R_IJ.setValue(0);
413
414 if(o_ == 0) { //------------- %*% ---------------------------
415
416 #define A_I_K mat_a.value [ i + k *a_nrow]
417 #define B_K_J mat_b.value [ k + j *b_nrow]
418 K_LOOP
419 #undef A_I_K
420 #undef B_K_J
421
422 } else if(o_ == 1){//------------- crossprod ---------------------------
423
424 #define A_I_K mat_a.value [ k + i *a_nrow]
425 #define B_K_J mat_b.value [ k + j *b_nrow]
426 K_LOOP
427 #undef A_I_K
428 #undef B_K_J
429
430 } else {//(o_ == 2) ------------- tcrossprod ---------------------------
431
432 #define A_I_K mat_a.value [ i + k *a_nrow]
433 #define B_K_J mat_b.value [ j + k *b_nrow]
434 K_LOOP
435 #undef A_I_K
436 #undef B_K_J
437
438 }
439 }
440 if(useMod)
441 res.modulus.push_back(biginteger(common_modulus));
442
443 mpz_clear(tt);
444 mpz_clear(common_modulus);
445
446 return( bigintegerR::create_SEXP (res));
447 } // matrix_mul_z()
448 #undef R_IJ
449 #undef K_LOOP
450
451
biginteger_rbind(SEXP args)452 SEXP biginteger_rbind(SEXP args)
453 {
454 int i=0,j=0;
455 bigvec result;
456 bigvec v;
457
458 result = bigintegerR::create_bignum(VECTOR_ELT(args,0));
459 if(result.nrow==0)
460 result.nrow = result.size();
461
462 result = matrixz::bigint_transpose(result, result.nrow,
463 result.size() / result.nrow);
464 for(i=1; i < LENGTH(args); i++)
465 {
466 v = bigintegerR::create_bignum(VECTOR_ELT(args,i));
467 if(v.nrow == 0 )
468 v.nrow = v.size();
469 v = matrixz::bigint_transpose(v,v.nrow,v.size() / v.nrow);
470
471 for(j=0; j< (int)v.size(); j++)
472 result.push_back(v[j]);
473 v.clear();
474 }
475
476 result = matrixz::bigint_transpose(result, result.nrow,
477 result.size() / result.nrow);
478 return bigintegerR::create_SEXP(result);
479 }
480
481
482
483 namespace matrixz
484 {
bigint_transpose(bigvec & mat,int nr,int nc)485 bigvec bigint_transpose ( bigvec & mat,int nr,int nc)
486 {
487 int i,j;
488
489 /* cas: square matrix */
490 bigvec matbis (nr * nc);
491 matbis.nrow = nc;
492
493 /* we compute transpose */
494 for(i =0; i<nr; i++)
495 for(j =0; j<nc; j++)
496 matbis.set(j+i*nc,mat[i+j*nr]);
497
498 return matbis;
499 }
500
501 /* return dimension in dima */
checkDims(int dima,int dimb)502 int checkDims(int dima, int dimb)
503 {
504 if(dima > 0 && dimb > 0) {
505 if (dimb != dima)
506 error(_("Matrix dimensions do not match"));
507 }
508 else { /* either a or b is a matrix */
509 if(dima == -1)
510 return(dimb);
511 }
512 return(dima);
513 }
514 }
515