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 741, Tue May 17 17:44:02 2005 UTC revision 758, Mon Jun 6 15:46:28 2005 UTC
# Line 373  Line 373 
373      return ans;      return ans;
374  }  }
375    
376  SEXP csc_matrix_mm(SEXP a, SEXP b)  SEXP csc_matrix_mm(SEXP a, SEXP b, SEXP classed, SEXP right)
377  {  {
378      int *adim = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int cl = asLogical(classed), rt = asLogical(right);
379        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
380        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
381          *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),          *ai = INTEGER(GET_SLOT(a, Matrix_iSym)),
382          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
383          *bdim = INTEGER(getAttrib(b, R_DimSymbol));          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
384      int j, k, m = adim[0], n = bdim[1], r = adim[1];                           getAttrib(b, R_DimSymbol)),
385      double *ax = REAL(GET_SLOT(a, Matrix_xSym));          *cdims = INTEGER(ALLOC_SLOT(val, Matrix_DimSym, INTSXP, 2)),
386      SEXP val;          chk, ione = 1, j, jj, k, m, n;
387        double *ax = REAL(GET_SLOT(a, Matrix_xSym)),
388            *bx = REAL(cl ? GET_SLOT(b, Matrix_xSym) : b), *cx;
389    
390      if (bdim[0] != r)      if (rt) {
391          error(_("Matrices of sizes (%d,%d) and (%d,%d) cannot be multiplied"),          m = bdims[0]; n = adims[1]; k = bdims[1]; chk = adims[0];
392                m, r, bdim[0], n);      } else {
393      val = PROTECT(allocMatrix(REALSXP, m, n));          m = adims[0]; n = bdims[1]; k = adims[1]; chk = bdims[0];
394      for (j = 0; j < n; j++) {   /* across columns of b */      }
395          double *ccol = REAL(val) + j * m,      if (chk != k)
396              *bcol = REAL(b) + j * r;          error(_("Matrices are not conformable for multiplication"));
397        if (m < 1 || n < 1 || k < 1)
398            error(_("Matrices with zero extents cannot be multiplied"));
399        cx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, m * n));
400        AZERO(cx, m * n); /* zero the accumulators */
401        for (j = 0; j < n; j++) { /* across columns of c */
402            if (rt) {
403                int kk, k2 = ap[j + 1];
404                for (kk = ap[j]; kk < k2; kk++) {
405                    F77_CALL(daxpy)(&m, &ax[kk], &bx[ai[kk]*m],
406                                    &ione, &cx[j*m], &ione);
407                }
408            } else {
409                double *ccol = cx + j * m,
410                    *bcol = bx + j * k;
411    
412          for (k = 0; k < m; k++) ccol[k] = 0.; /* zero the accumulators */              for (jj = 0; jj < k; jj++) { /* across columns of a */
413          for (k = 0; k < r; k++) { /* across columns of a */                  int kk, k2 = ap[jj + 1];
414              int kk, k2 = ap[k + 1];                  for (kk = ap[jj]; kk < k2; kk++) {
415              for (kk = ap[k]; kk < k2; kk++) {                      ccol[ai[kk]] += ax[kk] * bcol[jj];
416                  ccol[ai[kk]] += ax[kk] * bcol[k];                  }
417              }              }
418          }          }
419      }      }
420        cdims[0] = m; cdims[1] = n;
421      UNPROTECT(1);      UNPROTECT(1);
422      return val;      return val;
423  }  }

Legend:
Removed from v.741  
changed lines
  Added in v.758

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