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

revision 70, Mon Apr 12 12:10:01 2004 UTC revision 338, Fri Nov 12 21:17:36 2004 UTC
# Line 12  Line 12 
12          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix"))),          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 = INTEGER(GET_SLOT(ans, Matrix_DimSym)),          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
         *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));  
16    
17    
18        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
19        adims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
20      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
21      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')
22          SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));          SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
# Line 67  Line 69 
69      }      }
70      return ans;      return ans;
71  }  }
72    
73    SEXP Parent_inverse(SEXP par, SEXP unitdiag)
74    {
75        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("tscMatrix")));
76        int *ap, *ai, *dims, *pr = INTEGER(par),
77            countDiag = 1 - asLogical(unitdiag),
78            j, k, n = length(par), nnz, pos;
79        int *sz = Calloc(n, int);
80        double *ax;
81    
82        if (!isInteger(par)) error("par argument must be an integer vector");
83        for (j = n - 1; j >= 0; j--) {
84            int parent = pr[j];
85            sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);
86        }
87        SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));
88        ap = INTEGER(GET_SLOT(ans, Matrix_pSym));
89        ap[0] = 0;
90        for (j = 0; j < n; j++)
91            ap[j+1] = ap[j] + sz[j];
92        nnz = ap[n];
93        SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));
94        ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
95        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));
96        ax = REAL(GET_SLOT(ans, Matrix_xSym));
97        for (j = 0; j < nnz; j++) ax[j] = 1.;
98        SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));
99        dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));
100        dims[0] = dims[1] = n;
101        SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));
102        SET_SLOT(ans, Matrix_diagSym,
103                 (countDiag ? ScalarString(mkChar("N")) :
104                     ScalarString(mkChar("U"))));
105        pos = 0;
106        for (j = 0; j < n; j++) {
107            if (countDiag) ai[pos++] = j;
108            for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;
109        }
110        Free(sz);
111        UNPROTECT(1);
112        return ans;
113    }

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

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge