SCM

SCM Repository

[matrix] Diff of /branches/Matrix-mer2/src/dgCMatrix.c
ViewVC logotype

Diff of /branches/Matrix-mer2/src/dgCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 492, Thu Feb 3 14:24:03 2005 UTC revision 945, Wed Sep 28 08:54:28 2005 UTC
# Line 1  Line 1 
1  #include "dgCMatrix.h"  #include "dgCMatrix.h"
2    
3    #ifdef USE_CHOLMOD
4    #include "chm_common.h"
5    #endif
6    
7  SEXP dgCMatrix_validate(SEXP x)  SEXP dgCMatrix_validate(SEXP x)
8  {  {
9      SEXP pslot = GET_SLOT(x, Matrix_pSym),      SEXP pslot = GET_SLOT(x, Matrix_pSym),
# Line 14  Line 18 
18    
19      nrow = dims[0];      nrow = dims[0];
20      if (length(islot) != length(xslot))      if (length(islot) != length(xslot))
21          return ScalarString(mkChar("lengths of slots i and x must match"));          return mkString(_("lengths of slots i and x must match"));
22      if (length(pslot) <= 0)      if (length(pslot) <= 0)
23          return ScalarString(mkChar("slot p must have length > 0"));          return mkString(_("slot p must have length > 0"));
24      if (xp[0] != 0)      if (xp[0] != 0)
25          return ScalarString(mkChar("first element of slot p must be zero"));          return mkString(_("first element of slot p must be zero"));
26      if (length(islot) != xp[ncol])      if (length(islot) != xp[ncol])
27          return ScalarString(          return mkString(_("last element of slot p must match length of slots i and x"));
             mkChar(  
                 "last element of slot p must match length of slots i and x"));  
28      for (j = 0; j < ncol; j++) {      for (j = 0; j < ncol; j++) {
29          if (xp[j] > xp[j+1])          if (xp[j] > xp[j+1])
30              return ScalarString(mkChar("slot p must be non-decreasing"));              return mkString(_("slot p must be non-decreasing"));
31      }      }
32      for (j = 0; j < length(islot); j++) {      for (j = 0; j < length(islot); j++) {
33          if (xi[j] < 0 || xi[j] >= nrow)          if (xi[j] < 0 || xi[j] >= nrow)
34              return ScalarString(              return mkString(_("all row indices must be between 0 and nrow-1"));
                 mkChar("all row indices must be between 0 and nrow-1"));  
35      }      }
36      if (csc_unsorted_columns(ncol, xp, xi)) {      if (csc_unsorted_columns(ncol, xp, xi))
37          csc_sort_columns(ncol, xp, xi, REAL(xslot));          csc_sort_columns(ncol, xp, xi, REAL(xslot));
38      }  
39      return ScalarLogical(1);      return ScalarLogical(1);
40  }  }
41    
# Line 105  Line 106 
106    
107  SEXP csc_tcrossprod(SEXP x)  SEXP csc_tcrossprod(SEXP x)
108  {  {
109    #ifdef USE_CHOLMOD
110        cholmod_sparse *cha = cholmod_aat(as_cholmod_sparse(x),
111            (int *) NULL, 0, 1, &c);
112    
113        cha->stype = -1;            /* set the symmetry */
114        cholmod_sort(cha, &c);      /* drop redundant entries */
115        return chm_sparse_to_SEXP(cha, -1);
116    #else
117    
118      SEXP pslot = GET_SLOT(x, Matrix_pSym),      SEXP pslot = GET_SLOT(x, Matrix_pSym),
119          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix")));
120      int *xp = INTEGER(pslot),      int *xp = INTEGER(pslot),
# Line 152  Line 162 
162      triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,      triplet_to_col(nrow, nrow, ntrip, iVal, jVal, xVal,
163                     ansp, itmp, xtmp);                     ansp, itmp, xtmp);
164      nnz = ansp[nrow];      nnz = ansp[nrow];
165      SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));      SET_SLOT(ans, Matrix_uploSym, mkString("L"));
166      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
167      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
168      Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);      Memcpy(INTEGER(GET_SLOT(ans, Matrix_iSym)), itmp, nnz);
# Line 162  Line 172 
172      Free(itmp); Free(xtmp); Free(iVal); Free(jVal); Free(xVal);      Free(itmp); Free(xtmp); Free(iVal); Free(jVal); Free(xVal);
173      UNPROTECT(1);      UNPROTECT(1);
174      return ans;      return ans;
175    #endif /* USE_CHOLMOD */
176  }  }
177    
178  SEXP csc_matrix_crossprod(SEXP x, SEXP y)  SEXP csc_matrix_crossprod(SEXP x, SEXP y, SEXP classed)
179  {  {
180      SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;      int cl = asLogical(classed);
181      int j,      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
182          *xp = INTEGER(pslot),      int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
183          *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),          *ydims = INTEGER(cl ? GET_SLOT(y, Matrix_DimSym) :
184          xncol = length(pslot) - 1,                           getAttrib(y, R_DimSymbol)),
185          xnrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],          *vdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
186          *ydims;      int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
187      double *xx = REAL(GET_SLOT(x, Matrix_xSym));          *xp = INTEGER(GET_SLOT(x, Matrix_pSym));
188        int j, k = xdims[0], m = xdims[1], n = ydims[1];
189      if (!(isMatrix(y) && isReal(y))) error("y must be a numeric matrix");      double *vx, *xx = REAL(GET_SLOT(x, Matrix_xSym)),
190      ydims = INTEGER(getAttrib(y, R_DimSymbol));          *yx = REAL(cl ? GET_SLOT(y, Matrix_xSym) : y);
191      if (xnrow != ydims[0]) error("x and y must have the same number of rows");  
192      ans = PROTECT(allocMatrix(REALSXP, xncol, ydims[1]));      if (!cl && !(isMatrix(y) && isReal(y)))
193      for (j = 0; j < ydims[1]; j++) {          error(_("y must be a numeric matrix"));
194          int i; double *ypt = REAL(y) + j * ydims[0];      if (ydims[0] != k)
195          for(i = 0; i < xncol; i++) {          error(_("x and y must have the same number of rows"));
196        if (m < 1 || n < 1 || k < 1)
197            error(_("Matrices with zero extents cannot be multiplied"));
198        vdims[0] = m; vdims[1] = n;
199        vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
200        for (j = 0; j < n; j++) {
201            int i; double *ypt = yx + j * k;
202            for(i = 0; i < m; i++) {
203              int ii; double accum = 0.;              int ii; double accum = 0.;
204              for (ii = xp[i]; ii < xp[i+1]; ii++) {              for (ii = xp[i]; ii < xp[i+1]; ii++) {
205                  accum += xx[ii] * ypt[xi[ii]];                  accum += xx[ii] * ypt[xi[ii]];
206              }              }
207              REAL(ans)[i + j * xncol] = accum;              vx[i + j * m] = accum;
208          }          }
209      }      }
210      UNPROTECT(1);      UNPROTECT(1);
211      return ans;      return val;
212  }  }
213    
214  SEXP csc_to_dgTMatrix(SEXP x)  SEXP compressed_to_dgTMatrix(SEXP x, SEXP colP)
215  {  {
216      SEXP      int col = asLogical(colP); /* 1 if "C"olumn compressed;  0 if "R"ow */
217          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),      SEXP indSym = col ? Matrix_iSym : Matrix_jSym;
218          dimslot = GET_SLOT(x, Matrix_DimSym),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
219          islot = GET_SLOT(x, Matrix_iSym),          indP = GET_SLOT(x, indSym),
220          pslot = GET_SLOT(x, Matrix_pSym);          pP = GET_SLOT(x, Matrix_pSym);
221      int *dims = INTEGER(dimslot), j, jj,      int npt = length(pP) - 1;
         *xp = INTEGER(pslot), *yj;  
222    
223      SET_SLOT(ans, Matrix_iSym, duplicate(islot));      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
     SET_SLOT(ans, Matrix_DimSym, duplicate(dimslot));  
224      SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));      SET_SLOT(ans, Matrix_xSym, duplicate(GET_SLOT(x, Matrix_xSym)));
225      SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, length(islot)));      SET_SLOT(ans, indSym, duplicate(indP));
226      yj = INTEGER(GET_SLOT(ans, Matrix_jSym));      expand_cmprPt(npt, INTEGER(pP),
227      jj = 0;                    INTEGER(ALLOC_SLOT(ans, col ? Matrix_jSym : Matrix_iSym,
228      for (j = 0; j < dims[1]; j++) {                                       INTSXP, length(indP))));
229          while (jj < xp[j + 1]) {      UNPROTECT(1);
230              yj[jj] = j;      return ans;
             jj++;  
         }  
231      }      }
232    
233    SEXP compressed_non_0_ij(SEXP x, SEXP colP)
234    {
235        int col = asLogical(colP); /* 1 if "C"olumn compressed;  0 if "R"ow */
236        SEXP ans, indSym = col ? Matrix_iSym : Matrix_jSym;
237        SEXP indP = GET_SLOT(x, indSym),
238            pP = GET_SLOT(x, Matrix_pSym);
239        int n_el = length(indP), i, *ij;
240    
241        ij = INTEGER(ans = PROTECT(allocMatrix(INTSXP, n_el, 2)));
242        /* expand the compressed margin to 'i' or 'j' : */
243        expand_cmprPt(length(pP) - 1, INTEGER(pP), &ij[col ? n_el : 0]);
244        /* and copy the other one: */
245        if (col)
246            for(i = 0; i < n_el; i++)
247                ij[i] = INTEGER(indP)[i];
248        else /* row compressed */
249            for(i = 0; i < n_el; i++)
250                ij[i + n_el] = INTEGER(indP)[i];
251    
252      UNPROTECT(1);      UNPROTECT(1);
253      return ans;      return ans;
254  }  }
# Line 275  Line 309 
309      double *vx;      double *vx;
310    
311      if (!(isMatrix(A) && isReal(A)))      if (!(isMatrix(A) && isReal(A)))
312          error("A must be a numeric matrix");          error(_("A must be a numeric matrix"));
313      nrow = adims[0]; ncol = adims[1];      nrow = adims[0]; ncol = adims[1];
314      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));      SET_SLOT(val, Matrix_factorSym, allocVector(VECSXP, 0));
315      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));      SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
# Line 352  Line 386 
386    
387  SEXP csc_transpose(SEXP x)  SEXP csc_transpose(SEXP x)
388  {  {
389      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix"))),  #ifdef USE_CHOLMOD
390          islot = GET_SLOT(x, Matrix_iSym);      cholmod_sparse *chx = as_cholmod_sparse(x);
391      int *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),      SEXP ans =
392          nnz = length(islot);          chm_sparse_to_SEXP(cholmod_transpose(chx, 1, &c), 1);
393        Free(chx);
394        return ans;
395    #else
396        SEXP xi = GET_SLOT(x, Matrix_iSym);
397        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
398        int *adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
399            *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
400            nz = length(xi);
401        int *xj = Calloc(nz, int);
402        SEXP adn = ALLOC_SLOT(ans, Matrix_DimNamesSym, VECSXP, 2),
403            xdn = GET_SLOT(x, Matrix_DimNamesSym);
404    
     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));  
     adims = INTEGER(GET_SLOT(ans, Matrix_DimSym));  
405      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
406      SET_SLOT(ans, Matrix_factorSym, allocVector(VECSXP, 0));      SET_VECTOR_ELT(adn, 0, VECTOR_ELT(xdn, 1));
407      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));      SET_VECTOR_ELT(adn, 1, VECTOR_ELT(xdn, 0));
408      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      triplet_to_col(adims[0], adims[1], nz,
409      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));                     expand_cmprPt(xdims[1], INTEGER(GET_SLOT(x, Matrix_pSym)), xj),
410      csc_components_transpose(xdims[0], xdims[1], nnz,                     INTEGER(xi),
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
411                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
412                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                     INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, adims[1] + 1)),
413                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                     INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
414                               REAL(GET_SLOT(ans, Matrix_xSym)));                     REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)));
415        Free(xj);
416      UNPROTECT(1);      UNPROTECT(1);
417      return ans;      return ans;
418    #endif /* USE_CHOLMOD */
419  }  }
420    
421  SEXP csc_matrix_mm(SEXP a, SEXP b)  SEXP csc_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
422  {  {
423      int *adim = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int cl = asLogical(classed), rt = asLogical(right);
424        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
425        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
426          *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),          *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),
427          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
428          *bdim = INTEGER(getAttrib(b, R_DimSymbol));          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
429      int j, k, m = adim[0], n = bdim[1], r = adim[1];                           getAttrib(b, R_DimSymbol)),
430      double *ax = REAL(GET_SLOT(a, Matrix_xSym));          *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)),
431      SEXP val;          chk, ione = 1, j, jj, k, m, n;
432        double *ax = REAL(GET_SLOT(a, Matrix_xSym)),
433      if (bdim[0] != r)          *bx = REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), *cx;
434          error("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied",  
435                m, r, bdim[0], n);      if (rt) {
436      val = PROTECT(allocMatrix(REALSXP, m, n));          m = bdims[0]; n = adims[1]; k = bdims[1]; chk = adims[0];
437      for (j = 0; j < n; j++) {   /* across columns of b */      } else {
438          double *ccol = REAL(val) + j * m,          m = adims[0]; n = bdims[1]; k = adims[1]; chk = bdims[0];
439              *bcol = REAL(b) + j * r;      }
440        if (chk != k)
441          for (k = 0; k < m; k++) ccol[k] = 0.; /* zero the accumulators */          error(_("Matrices are not conformable for multiplication"));
442          for (k = 0; k < r; k++) { /* across columns of a */      if (m < 1 || n < 1 || k < 1)
443              int kk, k2 = ap[k + 1];          error(_("Matrices with zero extents cannot be multiplied"));
444              for (kk = ap[k]; kk < k2; kk++) {      cx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
445                  ccol[ai[kk]] += ax[kk] * bcol[k];      AZERO(cx, m * n); /* zero the accumulators */
446        for (j = 0; j < n; j++) { /* across columns of c */
447            if (rt) {
448                int kk, k2 = ap[j + 1];
449                for (kk = ap[j]; kk < k2; kk++) {
450                    F77_CALL(daxpy)(&m, &ax[kk], &bx[ai[kk]*m],
451                                    &ione, &cx[j*m], &ione);
452                }
453            } else {
454                double *ccol = cx + j * m,
455                    *bcol = bx + j * k;
456    
457                for (jj = 0; jj < k; jj++) { /* across columns of a */
458                    int kk, k2 = ap[jj + 1];
459                    for (kk = ap[jj]; kk < k2; kk++) {
460                        ccol[ai[kk]] += ax[kk] * bcol[jj];
461                    }
462              }              }
463          }          }
464      }      }
465        cdims[0] = m; cdims[1] = n;
466      UNPROTECT(1);      UNPROTECT(1);
467      return val;      return val;
468  }  }
# Line 416  Line 478 
478      SET_SLOT(val, Matrix_DimSym, duplicate(tmp));      SET_SLOT(val, Matrix_DimSym, duplicate(tmp));
479      ncol = INTEGER(tmp)[1];      ncol = INTEGER(tmp)[1];
480      if (!(isInteger(perm) && length(perm) == ncol))      if (!(isInteger(perm) && length(perm) == ncol))
481          error("perm must be an integer vector of length %d",          error(_("perm must be an integer vector of length %d"),
482                ncol);                ncol);
483      prm = INTEGER(perm);      prm = INTEGER(perm);
484      if (!R_ldl_valid_perm(ncol, prm))      if (!R_ldl_valid_perm(ncol, prm))
485          error("perm is not a valid 0-based permutation");          error(_("perm is not a valid 0-based permutation"));
486      iperm = Calloc(ncol, int);      iperm = Calloc(ncol, int);
487      for (j = 0; j < ncol; j++) iperm[prm[j]] = j;      for (j = 0; j < ncol; j++) iperm[prm[j]] = j;
488      tmp = GET_SLOT(x, Matrix_pSym);      tmp = GET_SLOT(x, Matrix_pSym);
# Line 451  Line 513 
513      UNPROTECT(1);      UNPROTECT(1);
514      return val;      return val;
515  }  }
   
   
   

Legend:
Removed from v.492  
changed lines
  Added in v.945

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge