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 2216, Thu Jul 17 20:19:29 2008 UTC revision 2217, Thu Jul 17 20:20:26 2008 UTC
# Line 80  Line 80 
80  {  {
81      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
82      CSP A = AS_CSP(Csparse_diagU2N(a));      CSP A = AS_CSP(Csparse_diagU2N(a));
83        CSP eye = csp_eye(A->n);
84      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
85          bnz = 10 * A->n,        /* initial estimate of nnz in b */          bnz = 10 * A->n,        /* initial estimate of nnz in b */
86          lo = uplo_P(a)[0] == 'L', top;          lo = uplo_P(a)[0] == 'L', top;
87      /* These arrays must use Calloc because of possible Realloc */      /* These arrays must use Calloc because of possible Realloc */
88      int *ti = Calloc(bnz, int), p, j, nz, pos = 0;      int *ti = Calloc(bnz, int), pos = 0;
89      double *tx = Calloc(bnz, double);      double *tx = Calloc(bnz, double);
     cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */  
90      double  *wrk = Alloca(A->n, double);      double  *wrk = Alloca(A->n, double);
91      int *xi = Alloca(2*A->n, int);      /* for cs_reach */      int *xi = Alloca(2*A->n, int);      /* for cs_reach */
92      R_CheckStack();      R_CheckStack();
# Line 95  Line 95 
95      SET_DimNames(ans, a);      SET_DimNames(ans, a);
96      slot_dup(ans, a, Matrix_uploSym);      slot_dup(ans, a, Matrix_uploSym);
97      slot_dup(ans, a, Matrix_diagSym);      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;  
98      bp[0] = 0;      bp[0] = 0;
99      for (j = 0; j < A->n; j++) {      for (int k = 0; k < A->n; k++)
100          u->i[0] = j;                    /* u := j'th unit vector */          col_spsolve(A, eye, k, xi, wrk, (int*)NULL, lo, bnz, pos, bp, ti, tx);
         /* (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]];  
             }  
     }  
101      nz = bp[A->n];      nz = bp[A->n];
102      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
103      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
# Line 130  Line 107 
107      return ans;      return ans;
108  }  }
109    
110    #endif
111    
112  SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)  SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
113  {  {
114      int cl = asLogical(classed);      int cl = asLogical(classed);
# Line 156  Line 135 
135    
136  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
137  {  {
138      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix")));
139      CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));      CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));
140      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
141          xnz = 10 * B->p[B->n],  /* initial estimate of nnz in x */          xnz = 10 * B->p[B->n];  /* initial estimate of nnz in x */
142          lo = uplo_P(a)[0] == 'L', top;      int *ti = Calloc(xnz, int), k, lo = uplo_P(a)[0] == 'L', pos = 0;
     /* These arrays must use Calloc because of possible Realloc */  
     int *ti = Calloc(xnz, int), p, j, nz, pos = 0;  
143      double *tx = Calloc(xnz, double);      double *tx = Calloc(xnz, double);
144      double  *wrk = Alloca(A->n, double);      double  *wrk = Alloca(A->n, double);
145      int *xi = Alloca(2*A->n, int);      /* for cs_reach */      int *xi = Alloca(2*A->n, int);      /* for cs_reach */
# Line 173  Line 150 
150      slot_dup(ans, b, Matrix_DimSym);      slot_dup(ans, b, Matrix_DimSym);
151      SET_DimNames(ans, b);      SET_DimNames(ans, b);
152      xp[0] = 0;      xp[0] = 0;
153      for (j = 0; j < B->n; j++) {      for (k = 0; k < B->n; k++) {
154          /* (wrk[top:n],xi[top:n]) :=  A^{-1} B[,j] */          int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
155          top = cs_spsolve (A, B, j, xi, wrk, 0, lo);          int nz = A->n - top, p;
156          nz = A->n - top;  
157          xp[j + 1] = nz + xp[j];          xp[k + 1] = nz + xp[k];
158          if (xp[j + 1] > xnz) {          if (xp[k + 1] > xnz) {
159              while (xp[j + 1] > xnz) xnz *= 2;              while (xp[k + 1] > xnz) xnz *= 2;
160              ti = Realloc(ti, xnz, int);              ti = Realloc(ti, xnz, int);
161              tx = Realloc(tx, xnz, double);              tx = Realloc(tx, xnz, double);
162          }          }
163          if (lo)          if (lo)                 /* increasing row order */
164              for(p = top; p < A->n; p++, pos++) {              for(p = top; p < A->n; p++, pos++) {
165                  ti[pos] = xi[p];                  ti[pos] = xi[p];
166                  tx[pos] = wrk[xi[p]];                  tx[pos] = wrk[xi[p]];
167              }              }
168          else /* upper triagonal */          else                    /* decreasing order, reverse copy */
169              for(p = A->n - 1; p >= top; p--, pos++) {              for(p = A->n - 1; p >= top; p--, pos++) {
170                  ti[pos] = xi[p];                  ti[pos] = xi[p];
171                  tx[pos] = wrk[xi[p]];                  tx[pos] = wrk[xi[p]];
172              }              }
173      }      }
174      nz = xp[A->n];      xnz = xp[A->n];
175      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  xnz)), ti, xnz);
176      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
177    
178      Free(ti); Free(tx);      Free(ti); Free(tx);
179      UNPROTECT(1);      UNPROTECT(1);

Legend:
Removed from v.2216  
changed lines
  Added in v.2217

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