The R Project SVN R-packages

Rev

Rev 3952 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3952 Rev 4427
Line 4... Line 4...
4 
{
 4 
{
5 
    int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values ;
 5 
    int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values ;
6 
    double *x, *Bx, *Cx ;
 6 
    double *x, *Bx, *Cx ;
7 
    cs *C ;
 7 
    cs *C ;
8 
    if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ;	    /* check inputs */
 8 
    if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ;	    /* check inputs */
- 
 
 9 
    if (A->m != B->m || A->n != B->n) return (NULL) ;
9 
    m = A->m ; anz = A->p [A->n] ;
 10 
    m = A->m ; anz = A->p [A->n] ;
10 
    n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [n] ;
 11 
    n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [n] ;
11 
    w = cs_calloc (m, sizeof (int)) ;			    /* get workspace */
 12 
    w = cs_calloc (m, sizeof (int)) ;			    /* get workspace */
12 
    values = (A->x != NULL) && (Bx != NULL) ;
 13 
    values = (A->x != NULL) && (Bx != NULL) ;
13 
    x = values ? cs_malloc (m, sizeof (double)) : NULL ;    /* get workspace */
 14 
    x = values ? cs_malloc (m, sizeof (double)) : NULL ;    /* get workspace */
Line 914... Line 915...
914 
	x [Vi [p]] -= Vx [p] * tau ;
 915 
	x [Vi [p]] -= Vx [p] * tau ;
915 
    }
 916 
    }
916 
    return (1) ;
 917 
    return (1) ;
917 
}
 918 
}
918 
/* create a Householder reflection [v,beta,s]=house(x), overwrite x with v,
 919 
/* create a Householder reflection [v,beta,s]=house(x), overwrite x with v,
919 
 * where (I-beta*v*v')*x = s*x.  See Algo 5.1.1, Golub & Van Loan, 3rd ed. */
 920 
 * where (I-beta*v*v')*x = s*e1.  See Algo 5.1.1, Golub & Van Loan, 3rd ed. */
920 
double cs_house (double *x, double *beta, int n)
 921 
double cs_house (double *x, double *beta, int n)
921 
{
 922 
{
922 
    double s, sigma = 0 ;
 923 
    double s, sigma = 0 ;
923 
    int i ;
 924 
    int i ;
924 
    if (!x || !beta) return (-1) ;	    /* check inputs */
 925 
    if (!x || !beta) return (-1) ;	    /* check inputs */
Line 1254... Line 1255...
1254 
{
 1255 
{
1255 
    int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values, *Bi ;
 1256 
    int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values, *Bi ;
1256 
    double *x, *Bx, *Cx ;
 1257 
    double *x, *Bx, *Cx ;
1257 
    cs *C ;
 1258 
    cs *C ;
1258 
    if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ;	 /* check inputs */
 1259 
    if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ;	 /* check inputs */
- 
 
 1260 
    if (A->n != B->m) return (NULL) ;
1259 
    m = A->m ; anz = A->p [A->n] ;
 1261 
    m = A->m ; anz = A->p [A->n] ;
1260 
    n = B->n ; Bp = B->p ; Bi = B->i ; Bx = B->x ; bnz = Bp [n] ;
 1262 
    n = B->n ; Bp = B->p ; Bi = B->i ; Bx = B->x ; bnz = Bp [n] ;
1261 
    w = cs_calloc (m, sizeof (int)) ;			 /* get workspace */
 1263 
    w = cs_calloc (m, sizeof (int)) ;			 /* get workspace */
1262 
    values = (A->x != NULL) && (Bx != NULL) ;
 1264 
    values = (A->x != NULL) && (Bx != NULL) ;
1263 
    x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */
 1265 
    x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */
Line 1397... Line 1399...
1397 
    return (1) ;
 1399 
    return (1) ;
1398 
}
 1400 
}
1399 
/* sparse QR factorization [V,beta,pinv,R] = qr (A) */
 1401 
/* sparse QR factorization [V,beta,pinv,R] = qr (A) */
1400 
csn *cs_qr (const cs *A, const css *S)
 1402 
csn *cs_qr (const cs *A, const css *S)
1401 
{
 1403 
{
1402 
    double *Rx, *Vx, *Ax, *Beta, *x ;
 1404 
    double *Rx, *Vx, *Ax, *x,  *Beta ;
1403 
    int i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai,
 1405 
    int i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai,
1404 
	*parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ;
 1406 
	*parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ;
1405 
    cs *R, *V ;
 1407 
    cs *R, *V ;
1406 
    csn *N ;
 1408 
    csn *N ;
1407 
    if (!CS_CSC (A) || !S) return (NULL) ;
 1409 
    if (!CS_CSC (A) || !S) return (NULL) ;
Line 1845... Line 1847...
1845 
    return (cs_done (C, w, NULL, 1)) ;	/* success; free w and return C */
 1847 
    return (cs_done (C, w, NULL, 1)) ;	/* success; free w and return C */
1846 
}
 1848 
}
1847 
/* sparse Cholesky update/downdate, L*L' + sigma*w*w' (sigma = +1 or -1) */
 1849 
/* sparse Cholesky update/downdate, L*L' + sigma*w*w' (sigma = +1 or -1) */
1848 
int cs_updown (cs *L, int sigma, const cs *C, const int *parent)
 1850 
int cs_updown (cs *L, int sigma, const cs *C, const int *parent)
1849 
{
 1851 
{
1850 
    int p, f, j, *Lp, *Li, *Cp, *Ci ;
 1852 
    int n, p, f, j, *Lp, *Li, *Cp, *Ci ;
1851 
    double *Lx, *Cx, alpha, beta = 1, delta, gamma, w1, w2, *w, n,  beta2 = 1 ;
 1853 
    double *Lx, *Cx, alpha, beta = 1, delta, gamma, w1, w2, *w, beta2 = 1 ;
1852 
    if (!CS_CSC (L) || !CS_CSC (C) || !parent) return (0) ;  /* check inputs */
 1854 
    if (!CS_CSC (L) || !CS_CSC (C) || !parent) return (0) ;  /* check inputs */
1853 
    Lp = L->p ; Li = L->i ; Lx = L->x ; n = L->n ;
 1855 
    Lp = L->p ; Li = L->i ; Lx = L->x ; n = L->n ;
1854 
    Cp = C->p ; Ci = C->i ; Cx = C->x ;
 1856 
    Cp = C->p ; Ci = C->i ; Cx = C->x ;
1855 
    if ((p = Cp [0]) >= Cp [1]) return (1) ;	    /* return if C empty */
 1857 
    if ((p = Cp [0]) >= Cp [1]) return (1) ;	    /* return if C empty */
1856 
    w = cs_malloc (n, sizeof (double)) ;	    /* get workspace */
 1858 
    w = cs_malloc (n, sizeof (double)) ;	    /* get workspace */