The R Project SVN R-packages

Rev

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

Rev 4979 Rev 4994
Line 473... Line 473...
473 
			      (CHM_SP)NULL, (char*) NULL, &c))
 473 
			      (CHM_SP)NULL, (char*) NULL, &c))
474 
	error(_("cholmod_write_sparse returned error code"));
 474 
	error(_("cholmod_write_sparse returned error code"));
475 
    fclose(f);
 475 
    fclose(f);
476 
    return R_NilValue;
 476 
    return R_NilValue;
477 
}
 477 
}
- 
 
 478 
 
- 
 
 479 
 
- 
 
 480 
/**
- 
 
 481 
 * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
- 
 
 482 
 * cholmod_sparse factor (LDL = TRUE).
- 
 
 483 
 *
- 
 
 484 
 * @param n  dimension of the matrix.
- 
 
 485 
 * @param x_p  'p' (column pointer) slot contents
- 
 
 486 
 * @param x_x  'x' (non-zero entries) slot contents
- 
 
 487 
 * @param perm 'perm' (= permutation vector) slot contents
- 
 
 488 
 * @param resultKind a (SEXP) string indicating which kind of result is desired.
- 
 
 489 
 *
- 
 
 490 
 * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
- 
 
 491 
 */
- 
 
 492 
SEXP diag_tC_ptr(int n, int *x_p, double *x_x, int *perm, SEXP resultKind)
- 
 
 493 
/*                                ^^^^^^ FIXME[Generalize] to int / ... */
- 
 
 494 
{
- 
 
 495 
    const char* res_ch = CHAR(STRING_ELT(resultKind,0));
- 
 
 496 
    enum diag_kind { diag, diag_backpermuted, trace, prod, sum_log
- 
 
 497 
    } res_kind = ((!strcmp(res_ch, "trace")) ? trace :
- 
 
 498 
		  ((!strcmp(res_ch, "sumLog")) ? sum_log :
- 
 
 499 
		   ((!strcmp(res_ch, "prod")) ? prod :
- 
 
 500 
		    ((!strcmp(res_ch, "diag")) ? diag :
- 
 
 501 
		     ((!strcmp(res_ch, "diagBack")) ? diag_backpermuted :
- 
 
 502 
		      -1)))));
- 
 
 503 
    int i, n_x, i_from = 0;
- 
 
 504 
    SEXP ans = PROTECT(allocVector(REALSXP,
- 
 
 505 
/*                                 ^^^^  FIXME[Generalize] */
- 
 
 506 
				   (res_kind == diag ||
- 
 
 507 
				    res_kind == diag_backpermuted) ? n : 1));
- 
 
 508 
    double *v = REAL(ans);
- 
 
 509 
/*  ^^^^^^      ^^^^  FIXME[Generalize] */
- 
 
 510 
 
- 
 
 511 
#define for_DIAG(v_ASSIGN)						\
- 
 
 512 
    for(i = 0; i < n; i++, i_from += n_x) {				\
- 
 
 513 
	/* looking at i-th column */					\
- 
 
 514 
	n_x = x_p[i+1] - x_p[i];/* #{entries} in this column */	\
- 
 
 515 
	v_ASSIGN;							\
- 
 
 516 
    }
- 
 
 517 
 
- 
 
 518 
    /* NOTA BENE: we assume  -- uplo = "L" i.e. lower triangular matrix
- 
 
 519 
     *            for uplo = "U" (makes sense with a "dtCMatrix" !),
- 
 
 520 
     *            should use  x_x[i_from + (nx - 1)] instead of x_x[i_from],
- 
 
 521 
     *            where nx = (x_p[i+1] - x_p[i])
- 
 
 522 
     */
- 
 
 523 
 
- 
 
 524 
    switch(res_kind) {
- 
 
 525 
    case trace:
- 
 
 526 
	v[0] = 0.;
- 
 
 527 
	for_DIAG(v[0] += x_x[i_from]);
- 
 
 528 
	break;
- 
 
 529 
 
- 
 
 530 
    case sum_log:
- 
 
 531 
	v[0] = 0.;
- 
 
 532 
	for_DIAG(v[0] += log(x_x[i_from]));
- 
 
 533 
	break;
- 
 
 534 
 
- 
 
 535 
    case prod:
- 
 
 536 
	v[0] = 1.;
- 
 
 537 
	for_DIAG(v[0] *= x_x[i_from]);
- 
 
 538 
	break;
- 
 
 539 
 
- 
 
 540 
    case diag:
- 
 
 541 
	for_DIAG(v[i] = x_x[i_from]);
- 
 
 542 
	break;
- 
 
 543 
 
- 
 
 544 
    case diag_backpermuted:
- 
 
 545 
	for_DIAG(v[i] = x_x[i_from]);
- 
 
 546 
	/* now back_permute : */
- 
 
 547 
 
- 
 
 548 
	error(_("resultKind = 'diagBack' (back-permuted) is not yet implemented"));
- 
 
 549 
	break;
- 
 
 550 
 
- 
 
 551 
    default: /* -1 from above */
- 
 
 552 
	error("diag_tC(): invalid 'resultKind'");
- 
 
 553 
	/* Wall: */ ans = R_NilValue; v = REAL(ans);
- 
 
 554 
    }
- 
 
 555 
 
- 
 
 556 
    UNPROTECT(1);
- 
 
 557 
    return ans;
- 
 
 558 
}
- 
 
 559 
 
- 
 
 560 
/**
- 
 
 561 
 * Extract the diagonal entries from *triangular* Csparse matrix  __or__ a
- 
 
 562 
 * cholmod_sparse factor (LDL = TRUE).
- 
 
 563 
 *
- 
 
 564 
 * @param pslot  'p' (column pointer)   slot of Csparse matrix/factor
- 
 
 565 
 * @param xslot  'x' (non-zero entries) slot of Csparse matrix/factor
- 
 
 566 
 * @param perm_slot  'perm' (= permutation vector) slot of corresponding CHMfactor
- 
 
 567 
 * @param resultKind a (SEXP) string indicating which kind of result is desired.
- 
 
 568 
 *
- 
 
 569 
 * @return  a SEXP, either a (double) number or a length n-vector of diagonal entries
- 
 
 570 
 */
- 
 
 571 
SEXP diag_tC(SEXP pslot, SEXP xslot, SEXP perm_slot, SEXP resultKind)
- 
 
 572 
{
- 
 
 573 
    int n = length(pslot) - 1, /* n = ncol(.) = nrow(.) */
- 
 
 574 
	*x_p  = INTEGER(pslot),
- 
 
 575 
	*perm = INTEGER(perm_slot);
- 
 
 576 
    double *x_x = REAL(xslot);
- 
 
 577 
/*  ^^^^^^        ^^^^ FIXME[Generalize] to INTEGER(.) / LOGICAL(.) / ... xslot !*/
- 
 
 578 
 
- 
 
 579 
    return diag_tC_ptr(n, x_p, x_x, perm, resultKind);
- 
 
 580 
}