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 677, Mon Mar 28 14:03:32 2005 UTC revision 1248, Thu Apr 13 22:05:22 2006 UTC
# Line 1  Line 1 
1                                  /* Sparse triangular matrices */                                  /* Sparse triangular numeric matrices */
2  #include "dtCMatrix.h"  #include "dtCMatrix.h"
3    
4  SEXP tsc_validate(SEXP x)  SEXP tsc_validate(SEXP x)
5  {  {
6      SEXP val;      return triangularMatrix_validate(x);
7        /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */
     if (isString(val = check_scalar_string(GET_SLOT(x, Matrix_uploSym),  
                                            "LU", "uplo"))) return val;  
     if (isString(val = check_scalar_string(GET_SLOT(x, Matrix_diagSym),  
                                            "NU", "diag"))) return val;  
     return ScalarLogical(1);  
8  }  }
9    
10    #if 0
11  SEXP tsc_transpose(SEXP x)  SEXP tsc_transpose(SEXP x)
12  {  {
13      SEXP      cholmod_sparse *cx = as_cholmod_sparse(x);
14          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),  
15        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),
16          islot = GET_SLOT(x, Matrix_iSym);          islot = GET_SLOT(x, Matrix_iSym);
17      int nnz = length(islot),      int nnz = length(islot),
18          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
19        int up = uplo_P(x)[0] == 'U';
20    
21      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
22      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
23      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')  
24          SET_SLOT(ans, Matrix_uploSym, mkString("L"));      if(*diag_P(x) == 'U')
25            SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));
26        SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
27    
28      csc_compTr(xdims[0], xdims[1], nnz,      csc_compTr(xdims[0], xdims[1], nnz,
29                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
30                 REAL(GET_SLOT(x, Matrix_xSym)),                 REAL(GET_SLOT(x, Matrix_xSym)),
# Line 33  Line 34 
34      UNPROTECT(1);      UNPROTECT(1);
35      return ans;      return ans;
36  }  }
37    #endif
38    
39  SEXP tsc_to_dgTMatrix(SEXP x)  SEXP tsc_to_dgTMatrix(SEXP x)
40  {  {
41      SEXP ans;      SEXP ans;
42      if (CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] != 'U')      if (*diag_P(x) != 'U')
43          ans = compressed_to_dgTMatrix(x, ScalarLogical(1));          ans = compressed_to_dgTMatrix(x, ScalarLogical(1));
44      else {                      /* unit triangular matrix */      else {                      /* unit triangular matrix */
45          SEXP islot = GET_SLOT(x, Matrix_iSym),          SEXP islot = GET_SLOT(x, Matrix_iSym),
# Line 138  Line 140 
140      UNPROTECT(1);      UNPROTECT(1);
141      return ans;      return ans;
142  }  }
143    
144    SEXP dtCMatrix_solve(SEXP a)
145    {
146        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
147        int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
148            n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,
149            *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
150            *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,
151            *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
152        int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
153        int *ri = Calloc(bnz, int), *ind = Calloc(n, int), j;
154        double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx;
155    
156        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
157        SET_SLOT(ans, Matrix_DimNamesSym,
158                 duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
159        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
160        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
161    
162        if (!(lo && unit))
163            error("code for non-unit or upper triangular not yet written");
164        /* Initially bp will contain increasing negative values ending at zero. */
165        /* Later we add the negative of bp[0] to all values. */
166        bp[n] = 0;
167        for (j = n - 1; j >= 0; j--) { /* columns in reverse order */
168            int i, i1 = ap[j], i2 = ap[j + 1], k, nr;
169            if (i1 < i2) AZERO(ind, n);
170            for (i = i1; i < i2; i++) {
171                ind[ai[i]] = 1;
172                for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1;
173            }
174            for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++;
175            if ((nr - bp[j + 1]) > bnz) {
176                while (nr > (bnz + bp[j + 1])) bnz *= 2;
177                ri = Realloc(ri, bnz, int);
178            }
179            bp[j] = bp[j + 1] - nr;
180            for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k;
181        }
182        bnz = -bp[0];
183        bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz));
184        bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz));
185        for (j = 0; j < n; j++) {
186            int bpnew = bp[j] + bnz;
187            Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]);
188            bp[j] = bpnew;
189        }
190        /* insert code for calculating the actual values here */
191        for (j = 0; j < bnz; j++) bx[j] = NA_REAL;
192    
193        Free(ind); Free(ri);
194        UNPROTECT(1);
195        return ans;
196    }

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

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