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

Legend:
Removed from v.338  
changed lines
  Added in v.945

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