/************************************************************/ /*! \file matrix.cc * \brief C++ function to add matrix support * * \version 1 * * \date Created: 19/02/06 * \date Last modified: Time-stamp: <2006-06-17 23:10:44 antoine> * * \author A. Lucas * * \note * as usually, matrix x[i,j] (n x p) is represented by a vector * x[i + j x n] (i=0..n-1 ; j=0..p-1) * * \note Licence: GPL (>= 2) */ #include using namespace std; #include "bigintegerR.h" #include "matrix.h" // need to call matrix_mul_q() #include "matrixq.h" // given that x is "bigz" or "bigq", // return TRUE if x is a bigz/q *matrix*: R's is.matrixZQ(.) SEXP is_matrix_zq(SEXP x) { SEXP nrowSexp = Rf_mkString("nrow"); PROTECT(nrowSexp); SEXP attributeRow = Rf_getAttrib(x,nrowSexp ); PROTECT(attributeRow); SEXP ans = Rf_ScalarLogical(attributeRow != R_NilValue); UNPROTECT(2); return ans; } // C++ side of R function matrix.bigz() SEXP as_matrixz (SEXP x, SEXP nrR, SEXP ncR, SEXP byrowR, SEXP mod) { int nc=INTEGER(ncR)[0], nr=INTEGER(nrR)[0], byrow=INTEGER(byrowR)[0]; // get "bigz" vector, this makes all conversion int to bigz etc... bigvec mat = bigintegerR::create_bignum(x); int lendat = mat.value.size(); // int sizemod = mat.modulus.size(); // when modulus specified bigvec modulus = bigintegerR::create_bignum(mod); if(modulus.value.size()>0) // should be allways the case if(!modulus.value[0].isNA()) { mat.modulus.resize(modulus.size()); for (unsigned int i = 0; i < modulus.size(); i++) mat.modulus[i] = modulus.value[i]; // sizemod = modulus.size(); } // A copy from R project to get correct dimension // all warnings... // if (nr == NA_INTEGER) // This is < 0 error(_("matrix: invalid 'nrow' value (too large or NA)")); if (nr < 0) error(_("matrix: invalid 'nrow' value (< 0)")); if (nc < 0) error(_("matrix: invalid 'ncol' value (< 0)")); if (nc == NA_INTEGER) error(_("matrix: invalid 'ncol' value (too large or NA)")); if(lendat > 0 ) { if (lendat > 1 && (nr * nc) % lendat != 0) { if (((lendat > nr) && (lendat / nr) * nr != lendat) || ((lendat < nr) && (nr / lendat) * lendat != nr)) warning(_("data length [%d] is not a sub-multiple or multiple of the number of rows [%d] in matrix"), lendat, nr); else if (((lendat > nc) && (lendat / nc) * nc != lendat) || ((lendat < nc) && (nc / lendat) * lendat != nc)) warning(_("data length [%d] is not a sub-multiple or multiple of the number of columns [%d] in matrix"), lendat, nc); } else if ((lendat > 1) && (nr * nc == 0)){ warning(_("data length exceeds size of matrix")); } } // update dimension parameters if(nr == 1) nr = (int)ceil(lendat / (double) nc); if(nc == 1) nc = (int)ceil(lendat / (double)nr); // when we extend "x" if(nc*nr > lendat) { mat.value.resize(nr*nc); for(int i = lendat; i < nr*nc; i++) mat.value[i] = mat.value[i % lendat]; } mat.nrow = nr; if(byrow) { bigvec mat2 = matrixz::bigint_transpose (mat, nc,nr); mat2.nrow = nr;// FIXME - needed ?? return( bigintegerR::create_SEXP (mat2)); } return( bigintegerR::create_SEXP (mat)); } /* * Transposition */ SEXP bigint_transposeR(SEXP x) { SEXP dimKey =Rf_mkString("nrow"); PROTECT(dimKey); SEXP dimAttr = Rf_getAttrib(x,dimKey ); PROTECT(dimAttr); bigvec mat = bigintegerR::create_bignum(x); int nr, n = mat.size(); if (dimAttr == R_NilValue) { // vector nr = n; } else if (TYPEOF(dimAttr) == INTSXP) { nr = INTEGER(dimAttr)[0]; } else { nr = -1;// -Wall error(_("argument must be a matrix of class \"bigz\"")); } UNPROTECT(2); int nc = (int) n / nr; // Rprintf(" o bigI_tr(<%d x %d>) ..\n", nr,nc); return( bigintegerR::create_SEXP(matrixz::bigint_transpose(mat, nr,nc))); } /* \brief matrix cross product * * \param a matrix (n x p) * \param trans if(trans), compute tcrossprod(), else crossprod() * \return crossprod(a) := t(a) %*% a [p x p] or * tcrossprod(a) := a %*% t(a) [n x n] */ SEXP matrix_crossp_z (SEXP a, SEXP trans) { bool useMod = FALSE, tr = (bool)Rf_asLogical(trans); bigvec mat_a = bigintegerR::create_bignum(a); int sizemod = mat_a.modulus.size(), a_nrow = mat_a.nrow, a_len = mat_a.size(); // in case of a vector; crossprod() returns scalar product, // whereas tcrossprod() gives n x n matrix. if(a_nrow < 0) a_nrow = a_len; int a_ncol = a_len / a_nrow; // Result R is R[1..m, 1..m] -- and R_{ij} = sum_{k=1}^p A.. B.. int m, p; if(tr) { // tcrossprod() m= a_nrow; p= a_ncol; } else { // crossprod() m= a_ncol; p= a_nrow; } bigvec res(m*m); res.nrow= m; mpz_t R_ij, tt; mpz_init(R_ij); mpz_init(tt); mpz_t common_modulus; mpz_init(common_modulus); if(sizemod <= 1) { // maybe 'useMod' i.e., can use common modulus: if(sizemod == 1) { mpz_set(common_modulus, mat_a.modulus[0].getValueTemp()); if(!mat_a.modulus[0].isNA()) useMod = TRUE; } } // here the computation: for(int i=0; i < m; i++) for(int j=0; j < m; j++) { mpz_set_ui(R_ij, 0); bool isna = false; #define K_LOOP \ for(int k=0; k < p; k++) { \ /* R_ij = \sum_{k=1}^p a_{ik} b_{kj} */ \ if( !(A_I_K.isNA() || B_K_J.isNA())) { \ mpz_mul(tt, A_I_K.getValueTemp(), B_K_J.getValueTemp()); \ mpz_add(R_ij, tt,R_ij); \ } \ else { \ isna = true; break; \ } \ } if(tr) {//------------- tcrossprod --------------------------- #define A_I_K mat_a.value [ i + k *a_nrow] #define B_K_J mat_a.value [ j + k *a_nrow] K_LOOP #undef A_I_K #undef B_K_J } else {//------------- crossprod --------------------------- #define A_I_K mat_a.value [ k + i *a_nrow] #define B_K_J mat_a.value [ k + j *a_nrow] K_LOOP #undef A_I_K #undef B_K_J } if(isna) { res.value[i + j*m].setValue(0); res.value[i + j*m].NA(true); } else res.value[i + j*m].setValue(R_ij); } // for(i ..) for(j ..) if(useMod) res.modulus.push_back(biginteger(common_modulus)); mpz_clear(R_ij); mpz_clear(tt); mpz_clear(common_modulus); return( bigintegerR::create_SEXP (res)); } // matrix_crossp_z() #undef K_LOOP /* \brief matrix multiplication * * returns matrix multiplication T(a) %*% b or b %*% T(a) * \param a matrix * \param b matrix * \param op operation code: 0: %*%, 1: crossprod, 2: tcrossprod * (same codes as in R's do_matprod() in src/main/array.c ) */ SEXP matrix_mul_z (SEXP a, SEXP b, SEXP op) { if(!strcmp(class_P(b), "bigq")) { // b "bigq", --> use q arithm: return(matrix_mul_q(a, b, op)); } // FIXME: we may know that a is 'bigz' - but we don't know at all about b !! // ----- create_bignum(.) should be much more careful (better: have a careful option!) bool useMod = FALSE;// if(useMod) use a *common* modulus int o_ = Rf_asInteger(op); // INTEGER(op)[0] bigvec mat_a = bigintegerR::create_bignum(a), mat_b = bigintegerR::create_bignum(b); int sizemod_a = mat_a.modulus.size(), sizemod_b = mat_b.modulus.size(); int a_nrow = mat_a.nrow, a_len = mat_a.size(), b_nrow = mat_b.nrow, b_len = mat_b.size(), a_ncol = -1, b_ncol = -1;// -Wall // distinguish cases of vectors / matrices --------------------- if(a_nrow < 0) { if(b_nrow < 0) { // *both* are vectors if(o_ == 0) { a_nrow = 1; a_ncol = a_len; } else { a_nrow = a_len; a_ncol = 1; } b_nrow = b_len; b_ncol = 1; } else { // a : vector, b : matrix b_ncol = b_len / b_nrow; if(o_ == 0) { if (a_len == b_nrow) { /* x as row vector */ a_nrow = 1; a_ncol = b_nrow; /* == a_len */ } else if (b_nrow == 1) { /* x as col vector */ a_nrow = a_len; a_ncol = 1; } } else if(o_ == 1) { /* crossprod() */ if (a_len == b_nrow) { /* x is a col vector */ a_nrow = b_nrow; /* == a_len */ a_ncol = 1; } /* else if (b_nrow == 1) ... not being too tolerant to treat x as row vector, as t(x) *is* row vector */ } else { // o_ == 2 -- tcrossprod() if (a_len == b_ncol) { /* x as row vector */ a_nrow = 1; a_ncol = b_ncol; /* == a_len */ } else if (b_ncol == 1) { /* x as col vector */ a_nrow = a_len; a_ncol = 1; } } } } else if (b_nrow < 0) { // a : matrix, b : vector a_ncol = a_len / a_nrow; if (o_ == 0) { if (b_len == a_ncol) { /* y as col vector */ b_nrow = a_ncol; b_ncol = 1; } else if (a_ncol == 1) { /* y as row vector */ b_nrow = 1; b_ncol = b_len; } } else if (o_ == 1) { /* crossprod() */ if (b_len == a_nrow) { /* y is a col vector */ b_nrow = a_nrow; b_ncol = 1; } } else { /* tcrossprod -- y is a col vector */ b_nrow = b_len; b_ncol = 1; } } else { // a, b *both* matrices a_ncol = a_len / a_nrow; b_ncol = b_len / b_nrow; } if(((o_ == 0) && (a_ncol != b_nrow)) || ((o_ == 1) && (a_nrow != b_nrow)) || // crossprod() ((o_ == 2) && (a_ncol != b_ncol)) // tcrossprod() ) error(_("Matrix dimensions do not match")); // Result R is R[1..n, 1..m] -- and R_{ij} = sum_{k=1} ^ p A.. B.. int n,m, p; if(o_ == 0) { n= a_nrow; m= b_ncol; p= a_ncol;// = b_nrow }else if (o_ == 1) { n= a_ncol; m= b_ncol; p= a_nrow;// = b_nrow }else if (o_ == 2) { n= a_nrow; m= b_nrow; p= a_ncol;// = b_ncol }else { error(_("invalid 'op' code in matrix_mul_z()")); n = m = p = -1;// -Wall } bigvec res(n*m); res.nrow=n; mpz_t common_modulus, tt; mpz_init(tt); mpz_init(common_modulus); /* modulus when modulus are "global" (i.e. of size 1) and * either are the same, or only one of a or b is specified */ if( !(sizemod_a > 1 || sizemod_b > 1)) { if((sizemod_a == 1) && (sizemod_b == 1)) { mpz_set(common_modulus, mat_a.modulus[0].getValueTemp()); if(mpz_cmp(common_modulus, mat_b.modulus[0].getValueTemp()) == 0 && !mat_a.modulus[0].isNA()) useMod = TRUE; } else { // at least one of the sizemod_* is > 1 : if ((sizemod_a == 1) && !mat_a.modulus[0].isNA()) { mpz_set(common_modulus, mat_a[0].getModulus().getValueTemp()); useMod = TRUE; } else if ((sizemod_b == 1) && !mat_b.modulus[0].isNA()) { mpz_set(common_modulus, mat_b[0].getModulus().getValueTemp()); useMod = TRUE; } } } // bigmod tmp; // here the computation: for(int i=0; i < n; i++) for(int j=0; j < m; j++) { #define R_IJ res.value[ i + j*n] #define K_LOOP \ for(int k=0; k < p; k++) \ { \ if(A_I_K.isNA() || B_K_J.isNA()) { \ R_IJ.setValue(0); R_IJ.NA(true); \ break; \ } \ /* Z = A_I_K * B_K_J */ \ mpz_mul(tt, A_I_K.getValueTemp(), B_K_J.getValueTemp()); \ /* R_IJ = R_IJ + A_I_K * B_K_J */ \ mpz_add(tt, tt, R_IJ.getValueTemp()); \ if(useMod) \ mpz_mod(tt,tt,common_modulus); \ R_IJ.setValue(tt); \ } R_IJ.setValue(0); if(o_ == 0) { //------------- %*% --------------------------- #define A_I_K mat_a.value [ i + k *a_nrow] #define B_K_J mat_b.value [ k + j *b_nrow] K_LOOP #undef A_I_K #undef B_K_J } else if(o_ == 1){//------------- crossprod --------------------------- #define A_I_K mat_a.value [ k + i *a_nrow] #define B_K_J mat_b.value [ k + j *b_nrow] K_LOOP #undef A_I_K #undef B_K_J } else {//(o_ == 2) ------------- tcrossprod --------------------------- #define A_I_K mat_a.value [ i + k *a_nrow] #define B_K_J mat_b.value [ j + k *b_nrow] K_LOOP #undef A_I_K #undef B_K_J } } if(useMod) res.modulus.push_back(biginteger(common_modulus)); mpz_clear(tt); mpz_clear(common_modulus); return( bigintegerR::create_SEXP (res)); } // matrix_mul_z() #undef R_IJ #undef K_LOOP SEXP biginteger_rbind(SEXP args) { int i=0,j=0; bigvec result; bigvec v; result = bigintegerR::create_bignum(VECTOR_ELT(args,0)); if(result.nrow==0) result.nrow = result.size(); result = matrixz::bigint_transpose(result, result.nrow, result.size() / result.nrow); for(i=1; i < LENGTH(args); i++) { v = bigintegerR::create_bignum(VECTOR_ELT(args,i)); if(v.nrow == 0 ) v.nrow = v.size(); v = matrixz::bigint_transpose(v,v.nrow,v.size() / v.nrow); for(j=0; j< (int)v.size(); j++) result.push_back(v[j]); v.clear(); } result = matrixz::bigint_transpose(result, result.nrow, result.size() / result.nrow); return bigintegerR::create_SEXP(result); } namespace matrixz { bigvec bigint_transpose ( bigvec & mat,int nr,int nc) { int i,j; /* cas: square matrix */ bigvec matbis (nr * nc); matbis.nrow = nc; /* we compute transpose */ for(i =0; i 0 && dimb > 0) { if (dimb != dima) error(_("Matrix dimensions do not match")); } else { /* either a or b is a matrix */ if(dima == -1) return(dimb); } return(dima); } }