SCM

SCM Repository

[matrix] Diff of /pkg/src/dgCMatrix.c
ViewVC logotype

Diff of /pkg/src/dgCMatrix.c

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

revision 682, Thu Mar 31 08:07:02 2005 UTC revision 683, Thu Mar 31 17:10:16 2005 UTC
# Line 161  Line 161 
161      return ans;      return ans;
162  }  }
163    
164  SEXP csc_matrix_crossprod(SEXP x, SEXP y)  SEXP csc_matrix_crossprod(SEXP x, SEXP y, SEXP classed)
165  {  {
166      SEXP pslot = GET_SLOT(x, Matrix_pSym), ans;      int cl = asLogical(classed);
167      int j,      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
168          *xp = INTEGER(pslot),      int *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym)),
169          *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),          *ydims = INTEGER(cl ? GET_SLOT(y, Matrix_DimSym) :
170          xncol = length(pslot) - 1,                           getAttrib(y, R_DimSymbol)),
171          xnrow = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],          *vdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2));
172          *ydims;      int *xi = INTEGER(GET_SLOT(x, Matrix_iSym)),
173      double *xx = REAL(GET_SLOT(x, Matrix_xSym));          *xp = INTEGER(GET_SLOT(x, Matrix_pSym));
174        int j, k = xdims[0], m = xdims[1], n = ydims[1];
175        double *vx, *xx = REAL(GET_SLOT(x, Matrix_xSym)),
176            *yx = REAL(cl ? GET_SLOT(y, Matrix_xSym) : y);
177    
178      if (!(isMatrix(y) && isReal(y))) error(_("y must be a numeric matrix"));      if (!cl && !(isMatrix(y) && isReal(y)))
179      ydims = INTEGER(getAttrib(y, R_DimSymbol));          error(_("y must be a numeric matrix"));
180      if (xnrow != ydims[0]) error(_("x and y must have the same number of rows"));      if (ydims[0] != k)
181      ans = PROTECT(allocMatrix(REALSXP, xncol, ydims[1]));          error(_("x and y must have the same number of rows"));
182      for (j = 0; j < ydims[1]; j++) {      if (m < 1 || n < 1 || k < 1)
183          int i; double *ypt = REAL(y) + j * ydims[0];          error(_("Matrices with zero extents cannot be multiplied"));
184          for(i = 0; i < xncol; i++) {      vdims[0] = m; vdims[1] = n;
185        vx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
186        for (j = 0; j < n; j++) {
187            int i; double *ypt = yx + j * k;
188            for(i = 0; i < m; i++) {
189              int ii; double accum = 0.;              int ii; double accum = 0.;
190              for (ii = xp[i]; ii < xp[i+1]; ii++) {              for (ii = xp[i]; ii < xp[i+1]; ii++) {
191                  accum += xx[ii] * ypt[xi[ii]];                  accum += xx[ii] * ypt[xi[ii]];
192              }              }
193              REAL(ans)[i + j * xncol] = accum;              vx[i + j * m] = accum;
194          }          }
195      }      }
196      UNPROTECT(1);      UNPROTECT(1);
197      return ans;      return val;
198  }  }
199    
200  SEXP compressed_to_dgTMatrix(SEXP x, SEXP colP)  SEXP compressed_to_dgTMatrix(SEXP x, SEXP colP)

Legend:
Removed from v.682  
changed lines
  Added in v.683

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