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 1248, Thu Apr 13 22:05:22 2006 UTC revision 1251, Sat Apr 15 13:08:26 2006 UTC
# Line 1  Line 1 
1                                  /* Sparse triangular numeric 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  {  {
# Line 141  Line 142 
142      return ans;      return ans;
143  }  }
144    
145    #if 0
146  SEXP dtCMatrix_solve(SEXP a)  SEXP dtCMatrix_solve(SEXP a)
147  {  {
148      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
# Line 188  Line 190 
190          bp[j] = bpnew;          bp[j] = bpnew;
191      }      }
192      /* insert code for calculating the actual values here */      /* insert code for calculating the actual values here */
193      for (j = 0; j < bnz; j++) bx[j] = NA_REAL;      for (j = 0; j < bnz; j++) bx[j] = 1;
194    
195      Free(ind); Free(ri);      Free(ind); Free(ri);
196      UNPROTECT(1);      UNPROTECT(1);
197      return ans;      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.1248  
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 University of Wisconsin - Madison Powered By FusionForge