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

pkg/src/dtCMatrix.c revision 2210, Wed Jul 16 22:20:22 2008 UTC pkg/Matrix/src/dtCMatrix.c revision 2889, Thu Aug 8 21:06:22 2013 UTC
# Line 74  Line 74 
74      }      }
75  }  }
76    
 #undef RETURN  
   
 SEXP dtCMatrix_solve(SEXP a)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));  
     CSP A = AS_CSP(Csparse_diagU2N(a));  
     int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),  
         bnz = 10 * A->n,        /* initial estimate of nnz in b */  
         lo = uplo_P(a)[0] == 'L', top;  
     /* These arrays must use Calloc because of possible Realloc */  
     int *ti = Calloc(bnz, int), p, j, nz, pos = 0;  
     double *tx = Calloc(bnz, double);  
     cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */  
     double  *wrk = Alloca(A->n, double);  
     int *xi = Alloca(2*A->n, int);      /* for cs_reach */  
     R_CheckStack();  
   
     slot_dup(ans, a, Matrix_DimSym);  
     SET_DimNames(ans, a);  
     slot_dup(ans, a, Matrix_uploSym);  
     slot_dup(ans, a, Matrix_diagSym);  
     /* initialize the "constant part" of the sparse unit vector */  
     u->x[0] = 1.;  
     u->p[0] = 0; u->p[1] = 1;  
     bp[0] = 0;  
     for (j = 0; j < A->n; j++) {  
         u->i[0] = j;                    /* u := j'th unit vector */  
         /* (wrk[top:n],xi[top:n]) :=  A^{-1} u  : */  
         top = cs_spsolve (A, u, 0, xi, wrk, 0, lo);  
         nz = A->n - top;  
         bp[j + 1] = nz + bp[j];  
         if (bp[j + 1] > bnz) {  
             while (bp[j + 1] > bnz) bnz *= 2;  
             ti = Realloc(ti, bnz, int);  
             tx = Realloc(tx, bnz, double);  
         }  
         if (lo)  
             for(p = top; p < A->n; p++, pos++) {  
                 ti[pos] = xi[p];  
                 tx[pos] = wrk[xi[p]];  
             }  
         else /* upper triagonal */  
             for(p = A->n - 1; p >= top; p--, pos++) {  
                 ti[pos] = xi[p];  
                 tx[pos] = wrk[xi[p]];  
             }  
     }  
     nz = bp[A->n];  
     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);  
     Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);  
   
     Free(ti); Free(tx); cs_spfree(u);  
     UNPROTECT(1);  
     return ans;  
 }  
   
77  SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)  SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
78  {  {
79      int cl = asLogical(classed);      int cl = asLogical(classed);
80      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
81      CSP A = AS_CSP(Csparse_diagU2N(a));      CSP A = AS_CSP(a);
82      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
84                           getAttrib(b, R_DimSymbol));                           getAttrib(b, R_DimSymbol));
# Line 150  Line 94 
94                  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);                  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
95      for (j = 0; j < nrhs; j++)      for (j = 0; j < nrhs; j++)
96          lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);          lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
97      UNPROTECT(1);      RETURN(ans);
     return ans;  
98  }  }
99    
100  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
101  {  {
102      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
103      CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));      CSP A = AS_CSP(a), B = AS_CSP(b);
     int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),  
         xnz = 10 * B->p[B->n],  /* initial estimate of nnz in x */  
         lo = uplo_P(a)[0] == 'L', top;  
     /* These arrays must use Calloc because of possible Realloc */  
     int *ti = Calloc(xnz, int), p, j, nz, pos = 0;  
     double *tx = Calloc(xnz, double);  
     double  *wrk = Alloca(A->n, double);  
     int *xi = Alloca(2*A->n, int);      /* for cs_reach */  
104      R_CheckStack();      R_CheckStack();
   
105      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
106          error(_("Dimensions of system to be solved are inconsistent"));          error(_("Dimensions of system to be solved are inconsistent"));
107        // *before* Calloc()ing below [memory leak]!
108    
109        int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
110            xnz = 10 * B->p[B->n];  /* initial estimate of nnz in x */
111        int k, lo = uplo_P(a)[0] == 'L', pos = 0;
112        int    *ti = Calloc(xnz, int),     *xi = Calloc(2*A->n, int); /* for cs_reach */
113        double *tx = Calloc(xnz, double), *wrk = Calloc(  A->n, double);
114    
115      slot_dup(ans, b, Matrix_DimSym);      slot_dup(ans, b, Matrix_DimSym);
116      SET_DimNames(ans, b);      SET_DimNames(ans, b);
117      xp[0] = 0;      xp[0] = 0;
118      for (j = 0; j < B->n; j++) {      for (k = 0; k < B->n; k++) {
119          /* (wrk[top:n],xi[top:n]) :=  A^{-1} B[,j] */          int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
120          top = cs_spsolve (A, B, j, xi, wrk, 0, lo);          int nz = A->n - top;
121          nz = A->n - top;  
122          xp[j + 1] = nz + xp[j];          xp[k + 1] = nz + xp[k];
123          if (xp[j + 1] > xnz) {          if (xp[k + 1] > xnz) {
124              while (xp[j + 1] > xnz) xnz *= 2;              while (xp[k + 1] > xnz) xnz *= 2;
125              ti = Realloc(ti, xnz, int);              ti = Realloc(ti, xnz, int);
126              tx = Realloc(tx, xnz, double);              tx = Realloc(tx, xnz, double);
127          }          }
128          if (lo)          if (lo)                 /* increasing row order */
129              for(p = top; p < A->n; p++, pos++) {              for(int p = top; p < A->n; p++, pos++) {
130                  ti[pos] = xi[p];                  ti[pos] = xi[p];
131                  tx[pos] = wrk[xi[p]];                  tx[pos] = wrk[xi[p]];
132              }              }
133          else /* upper triagonal */          else                    /* decreasing order, reverse copy */
134              for(p = A->n - 1; p >= top; p--, pos++) {              for(int p = A->n - 1; p >= top; p--, pos++) {
135                  ti[pos] = xi[p];                  ti[pos] = xi[p];
136                  tx[pos] = wrk[xi[p]];                  tx[pos] = wrk[xi[p]];
137              }              }
138      }      }
139      nz = xp[A->n];      xnz = xp[B->n];
140      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  xnz)), ti, xnz);
141      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
142    
143      Free(ti); Free(tx);      Free(ti); Free(tx);
144      UNPROTECT(1);      Free(wrk); Free(xi);
145      return ans;  
146        RETURN(ans);
147  }  }
148    #undef RETURN
149    

Legend:
Removed from v.2210  
changed lines
  Added in v.2889

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