The R Project SVN R-packages

Rev

Rev 1539 | Rev 1786 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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