SCM

SCM Repository

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

Diff of /pkg/src/dtCMatrix.c

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

revision 357, Sat Nov 27 17:56:13 2004 UTC revision 358, Sat Nov 27 17:57:19 2004 UTC
# Line 70  Line 70 
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)  SEXP Parent_inverse(SEXP par, SEXP unitdiag)
116  {  {
117      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));
118      int *ap, *ai, *dims, *pr = INTEGER(par),      int *ap, *ai, *dims, *pr = INTEGER(par),
119          countDiag = 1 - asLogical(unitdiag),          countDiag = 1 - asLogical(unitdiag),
120          j, k, n = length(par), nnz, pos;          j, k, n = length(par), nnz;
     int *sz = Calloc(n, int);  
121      double *ax;      double *ax;
122    
123      if (!isInteger(par)) error("par argument must be an integer vector");      if (!isInteger(par)) error("par argument must be an integer vector");
     for (j = n - 1; j >= 0; j--) {  
         int parent = pr[j];  
         sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);  
     }  
124      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
125      ap = INTEGER(GET_SLOT(ans, Matrix_pSym));      ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
126      ap[0] = 0;      nnz = parent_inv_ap(n, countDiag, pr, ap);
     for (j = 0; j < n; j++)  
         ap[j+1] = ap[j] + sz[j];  
     nnz = ap[n];  
127      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
128      ai = INTEGER(GET_SLOT(ans, Matrix_iSym));      ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
129      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
# Line 102  Line 136 
136      SET_SLOT(ans, Matrix_diagSym,      SET_SLOT(ans, Matrix_diagSym,
137               (countDiag ? ScalarString(mkChar("N")) :               (countDiag ? ScalarString(mkChar("N")) :
138                   ScalarString(mkChar("U"))));                   ScalarString(mkChar("U"))));
139      pos = 0;      parent_inv_ai(n, countDiag, pr, ai);
     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);  
140      UNPROTECT(1);      UNPROTECT(1);
141      return ans;      return ans;
142  }  }

Legend:
Removed from v.357  
changed lines
  Added in v.358

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