The R Project SVN R

Rev

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

Rev 60327 Rev 60335
Line 1565... Line 1565...
1565 
   UNPROTECT(nprotect);
 1565 
   UNPROTECT(nprotect);
1566 
   return ans;
 1566 
   return ans;
1567 
}
 1567 
}
1568 
 
 1568 
 
1569 
 
 1569 
 
1570 
/* backsolve(r, x, k, upper.tri, transpose) */
 1570 
/* backsolve(r, b, k, upper.tri, transpose) */
1571 
SEXP attribute_hidden do_backsolve(SEXP call, SEXP op, SEXP args, SEXP rho)
 1571 
SEXP attribute_hidden do_backsolve(SEXP call, SEXP op, SEXP args, SEXP rho)
1572 
{
 1572 
{
1573 
    int nprot = 1;
 1573 
    int nprot = 1;
1574 
    checkArity(op, args);
 1574 
    checkArity(op, args);
1575 
 
 1575 
 
1576 
    SEXP r = CAR(args); args = CDR(args);
 1576 
    SEXP r = CAR(args); args = CDR(args);
1577 
    SEXP x = CAR(args); args = CDR(args);
 1577 
    SEXP b = CAR(args); args = CDR(args);
1578 
    int nb = ncols(x), nrx = nrows(r);
 1578 
    int nrr = nrows(r), nrb = nrows(b), ncb = ncols(b);
1579 
    int k = asInteger(CAR(args)); args = CDR(args);
 1579 
    int k = asInteger(CAR(args)); args = CDR(args);
- 
 
 1580 
    /* k is the number of rows to be used: there must be at least that
- 
 
 1581 
       many rows and cols in the rhs and at least that many rows on
- 
 
 1582 
       the rhs.
- 
 
 1583 
    */
1580 
    if (k == NA_INTEGER || k <= 0 || k > nrx)
1584 
    if (k == NA_INTEGER || k <= 0 || k > nrr || k > ncols(r) || k > nrb)
1581 
	error(_("invalid '%s' argument"), "k");
 1585 
	error(_("invalid '%s' argument"), "k");
1582 
    int upper = asLogical(CAR(args)); args = CDR(args);
 1586 
    int upper = asLogical(CAR(args)); args = CDR(args);
1583 
    if (upper == NA_INTEGER) error(_("invalid '%s' argument"), "upper.tri");
 1587 
    if (upper == NA_INTEGER) error(_("invalid '%s' argument"), "upper.tri");
1584 
    int trans = asLogical(CAR(args));
 1588 
    int trans = asLogical(CAR(args));
1585 
    if (trans == NA_INTEGER) error(_("invalid '%s' argument"), "transpose");
 1589 
    if (trans == NA_INTEGER) error(_("invalid '%s' argument"), "transpose");
1586 
    SEXP ans;
 1590 
    SEXP ans;
1587 
    if (TYPEOF(r) != REALSXP) {PROTECT(r = coerceVector(r, REALSXP)); nprot++;}
 1591 
    if (TYPEOF(r) != REALSXP) {PROTECT(r = coerceVector(r, REALSXP)); nprot++;}
1588 
    if (TYPEOF(x) != REALSXP) {PROTECT(x = coerceVector(x, REALSXP)); nprot++;}
 1592 
    if (TYPEOF(b) != REALSXP) {PROTECT(b = coerceVector(b, REALSXP)); nprot++;}
- 
 
 1593 
    double *rr = REAL(r);
- 
 
 1594 
 
- 
 
 1595 
    /* check for zeros on diagonal of r: only k row/cols are used. */
- 
 
 1596 
    size_t incr = nrr + 1;
- 
 
 1597 
    for(int i = 0; i < k; i++) { /* check for zeros on diagonal */
- 
 
 1598 
	if (rr[i * incr] == 0.0)
- 
 
 1599 
	    error(_("singular matrix in 'backsolve'. First zero in diagonal [%d]"),
- 
 
 1600 
		  i + 1);
- 
 
 1601 
    }
- 
 
 1602 
 
1589 
    PROTECT(ans = allocMatrix(REALSXP, k, nb));
 1603 
    PROTECT(ans = allocMatrix(REALSXP, k, ncb));
- 
 
 1604 
    if (k > 0 && ncb > 0) {
1590 
    int info = 0, job = upper + 10 * trans;
 1605 
       /* copy (part) cols of b to ans */
- 
 
 1606 
	for(R_xlen_t j = 0; j < ncb; j++)
1591 
    bakslv(REAL(r), &nrx, &k, REAL(x), &k, &nb, REAL(ans), &job, &info);
1607 
	    memcpy(REAL(ans) + j*k, REAL(b) + j*nrb, (size_t)k *sizeof(double));
1592 
    if (info)
 1608 
	double one = 1.0;
1593 
	error(_("singular matrix in 'backsolve'. First zero in diagonal [%d]"),
 1609 
	F77_CALL(dtrsm)("L", upper ? "U" : "L", trans ? "T" : "N", "N",
- 
 
 1610 
			&k, &ncb, &one, rr, &nrr, REAL(ans), &k);
1594 
	    info);
 1611 
    }
1595 
    UNPROTECT(nprot);
 1612 
    UNPROTECT(nprot);
1596 
    return ans;
 1613 
    return ans;
1597 
}
 1614 
}
1598 
 
 1615 
 
1599 
/* max.col(m, ties.method) */
 1616 
/* max.col(m, ties.method) */