SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/src/dtCMatrix.c
ViewVC logotype

Diff of /pkg/Matrix/src/dtCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

pkg/src/tscMatrix.c revision 70, Mon Apr 12 12:10:01 2004 UTC pkg/src/dtCMatrix.c revision 677, Mon Mar 28 14:03:32 2005 UTC
# Line 1  Line 1 
1                                  /* Sparse triangular matrices */                                  /* Sparse triangular matrices */
2  #include "tscMatrix.h"  #include "dtCMatrix.h"
3    
4  SEXP tsc_validate(SEXP x)  SEXP tsc_validate(SEXP x)
5  {  {
6        SEXP val;
7    
8        if (isString(val = check_scalar_string(GET_SLOT(x, Matrix_uploSym),
9                                               "LU", "uplo"))) return val;
10        if (isString(val = check_scalar_string(GET_SLOT(x, Matrix_diagSym),
11                                               "NU", "diag"))) return val;
12      return ScalarLogical(1);      return ScalarLogical(1);
13  }  }
14    
15  SEXP tsc_transpose(SEXP x)  SEXP tsc_transpose(SEXP x)
16  {  {
17      SEXP      SEXP
18          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix"))),          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),
19          islot = GET_SLOT(x, Matrix_iSym);          islot = GET_SLOT(x, Matrix_iSym);
20      int nnz = length(islot),      int nnz = length(islot),
21          *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
         *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));  
22    
23        adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
24      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
25      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')
26          SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));          SET_SLOT(ans, Matrix_uploSym, mkString("L"));
27      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));      csc_compTr(xdims[0], xdims[1], nnz,
28      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));  
     csc_components_transpose(xdims[0], xdims[1], nnz,  
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
29                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
30                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
31                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
32                               REAL(GET_SLOT(ans, Matrix_xSym)));                 REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
33      UNPROTECT(1);      UNPROTECT(1);
34      return ans;      return ans;
35  }  }
36    
37  SEXP tsc_to_triplet(SEXP x)  SEXP tsc_to_dgTMatrix(SEXP x)
38  {  {
39      SEXP ans;      SEXP ans;
40      if (toupper(CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0]) != 'U')      if (CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] != 'U')
41          ans = csc_to_triplet(x);          ans = compressed_to_dgTMatrix(x, ScalarLogical(1));
42      else {                      /* unit triangular matrix */      else {                      /* unit triangular matrix */
43          SEXP islot = GET_SLOT(x, Matrix_iSym),          SEXP islot = GET_SLOT(x, Matrix_iSym),
44              pslot = GET_SLOT(x, Matrix_pSym);              pslot = GET_SLOT(x, Matrix_pSym);
# Line 47  Line 49 
49              *p = INTEGER(pslot);              *p = INTEGER(pslot);
50          double *ax;          double *ax;
51    
52          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tripletMatrix")));          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));
53          SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));          SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
54          SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));          SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
55          ai = INTEGER(GET_SLOT(ans, Matrix_iSym));          ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
# Line 67  Line 69 
69      }      }
70      return ans;      return ans;
71  }  }
72    
73    /**
74     * Derive the column pointer vector for the inverse of L from the parent array
75     *
76     * @param n length of parent array
77     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
78     * @param pr parent vector describing the elimination tree
79     * @param ap array of length n+1 to be filled with the column pointers
80     *
81     * @return the number of non-zero entries (ap[n])
82     */
83    int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
84    {
85        int *sz = Calloc(n, int), j;
86    
87        for (j = n - 1; j >= 0; j--) {
88            int parent = pr[j];
89            sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);
90        }
91        ap[0] = 0;
92        for (j = 0; j < n; j++)
93            ap[j+1] = ap[j] + sz[j];
94        Free(sz);
95        return ap[n];
96    }
97    
98    /**
99     * Derive the row index array for the inverse of L from the parent array
100     *
101     * @param n length of parent array
102     * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1
103     * @param pr parent vector describing the elimination tree
104     * @param ai row index vector of length ap[n]
105     */
106    void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])
107    {
108        int j, k, pos = 0;
109        for (j = 0; j < n; j++) {
110            if (countDiag) ai[pos++] = j;
111            for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
112        }
113    }
114    
115    SEXP Parent_inverse(SEXP par, SEXP unitdiag)
116    {
117        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
118        int *ap, *ai, *dims, *pr = INTEGER(par),
119            countDiag = 1 - asLogical(unitdiag),
120            j, n = length(par), nnz;
121        double *ax;
122    
123        if (!isInteger(par)) error(_("par argument must be an integer vector"));
124        SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
125        ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
126        nnz = parent_inv_ap(n, countDiag, pr, ap);
127        SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
128        ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
129        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
130        ax = REAL(GET_SLOT(ans, Matrix_xSym));
131        for (j = 0; j < nnz; j++) ax[j] = 1.;
132        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
133        dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
134        dims[0] = dims[1] = n;
135        SET_SLOT(ans, Matrix_uploSym, mkString("L"));
136        SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));
137        parent_inv_ai(n, countDiag, pr, ai);
138        UNPROTECT(1);
139        return ans;
140    }

Legend:
Removed from v.70  
changed lines
  Added in v.677

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