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 1251, Sat Apr 15 13:08:26 2006 UTC
# Line 1  Line 1 
1                                  /* Sparse triangular matrices */                                  /* Sparse triangular numeric matrices */
2  #include "dtCMatrix.h"  #include "dtCMatrix.h"
3    #include "cs_utils.h"
4    
5  SEXP tsc_validate(SEXP x)  SEXP tsc_validate(SEXP x)
6  {  {
7      SEXP val;      return triangularMatrix_validate(x);
8        /* 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);  
9  }  }
10    
11    #if 0
12  SEXP tsc_transpose(SEXP x)  SEXP tsc_transpose(SEXP x)
13  {  {
14      SEXP      cholmod_sparse *cx = as_cholmod_sparse(x);
15          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),  
16        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),
17          islot = GET_SLOT(x, Matrix_iSym);          islot = GET_SLOT(x, Matrix_iSym);
18      int nnz = length(islot),      int nnz = length(islot),
19          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
20        int up = uplo_P(x)[0] == 'U';
21    
22      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
23      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
24      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')  
25          SET_SLOT(ans, Matrix_uploSym, mkString("L"));      if(*diag_P(x) == 'U')
26            SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));
27        SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));
28    
29      csc_compTr(xdims[0], xdims[1], nnz,      csc_compTr(xdims[0], xdims[1], nnz,
30                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
31                 REAL(GET_SLOT(x, Matrix_xSym)),                 REAL(GET_SLOT(x, Matrix_xSym)),
# Line 33  Line 35 
35      UNPROTECT(1);      UNPROTECT(1);
36      return ans;      return ans;
37  }  }
38    #endif
39    
40  SEXP tsc_to_dgTMatrix(SEXP x)  SEXP tsc_to_dgTMatrix(SEXP x)
41  {  {
42      SEXP ans;      SEXP ans;
43      if (CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] != 'U')      if (*diag_P(x) != 'U')
44          ans = compressed_to_dgTMatrix(x, ScalarLogical(1));          ans = compressed_to_dgTMatrix(x, ScalarLogical(1));
45      else {                      /* unit triangular matrix */      else {                      /* unit triangular matrix */
46          SEXP islot = GET_SLOT(x, Matrix_iSym),          SEXP islot = GET_SLOT(x, Matrix_iSym),
# Line 138  Line 141 
141      UNPROTECT(1);      UNPROTECT(1);
142      return ans;      return ans;
143  }  }
144    
145    #if 0
146    SEXP dtCMatrix_solve(SEXP a)
147    {
148        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
149        int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
150            n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,
151            *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
152            *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,
153            *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
154        int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
155        int *ri = Calloc(bnz, int), *ind = Calloc(n, int), j;
156        double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx;
157    
158        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
159        SET_SLOT(ans, Matrix_DimNamesSym,
160                 duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
161        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
162        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
163    
164        if (!(lo && unit))
165            error("code for non-unit or upper triangular not yet written");
166        /* Initially bp will contain increasing negative values ending at zero. */
167        /* Later we add the negative of bp[0] to all values. */
168        bp[n] = 0;
169        for (j = n - 1; j >= 0; j--) { /* columns in reverse order */
170            int i, i1 = ap[j], i2 = ap[j + 1], k, nr;
171            if (i1 < i2) AZERO(ind, n);
172            for (i = i1; i < i2; i++) {
173                ind[ai[i]] = 1;
174                for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1;
175            }
176            for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++;
177            if ((nr - bp[j + 1]) > bnz) {
178                while (nr > (bnz + bp[j + 1])) bnz *= 2;
179                ri = Realloc(ri, bnz, int);
180            }
181            bp[j] = bp[j + 1] - nr;
182            for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k;
183        }
184        bnz = -bp[0];
185        bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz));
186        bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz));
187        for (j = 0; j < n; j++) {
188            int bpnew = bp[j] + bnz;
189            Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]);
190            bp[j] = bpnew;
191        }
192        /* insert code for calculating the actual values here */
193        for (j = 0; j < bnz; j++) bx[j] = 1;
194    
195        Free(ind); Free(ri);
196        UNPROTECT(1);
197        return ans;
198    }
199    #else
200    SEXP dtCMatrix_solve(SEXP a)
201    {
202        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
203        cs *A = Matrix_as_cs(a);
204        int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
205            lo = uplo_P(a)[0] == 'L',
206            bnz = 10 * A->n;        /* initial estimate of nnz in b */
207        int *ti = Calloc(bnz, int), i, j, nz, pos = 0;
208        double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
209    
210        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
211        SET_SLOT(ans, Matrix_DimNamesSym,
212                 duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
213        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
214        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
215        bp[0] = 0;
216        for (j = 0; j < A->n; j++) {
217            AZERO(wrk, A->n);
218            wrk[j] = 1;
219            lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);
220            for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;
221            bp[j + 1] = nz + bp[j];
222            if (bp[j + 1] > bnz) {
223                while (bp[j + 1] > bnz) bnz *= 2;
224                ti = Realloc(ti, bnz, int);
225                tx = Realloc(tx, bnz, double);
226            }
227            for (i = 0; i < A->n; i++)
228                if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}
229        }
230        nz = bp[A->n];
231        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
232        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
233    
234        Free(A); Free(ti); Free(tx);
235        UNPROTECT(1);
236        return ans;
237    }
238    #endif
239    
240    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
241    {
242        int cl = asLogical(classed);
243        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
244        cs *A = Matrix_as_cs(a);
245        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
246            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
247                             getAttrib(b, R_DimSymbol));
248        int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
249        double *bx;
250    
251        if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
252            error(_("Dimensions of system to be solved are inconsistent"));
253        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
254        /* copy dimnames or Dimnames as well */
255        bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
256                    REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
257        for (j = 0; j < nrhs; j++)
258            lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
259        Free(A);
260        UNPROTECT(1);
261        return ans;
262    }

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

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