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 2208, Sun Jul 13 15:51:45 2008 UTC revision 2210, Wed Jul 16 22:20:22 2008 UTC
# Line 74  Line 74 
74      }      }
75  }  }
76    
   
77  #undef RETURN  #undef RETURN
78    
 /**  
  * Derive the column pointer vector for the inverse of L from the parent array  
  *  
  * @param n length of parent array  
  * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1  
  * @param pr parent vector describing the elimination tree  
  * @param ap array of length n+1 to be filled with the column pointers  
  *  
  * @return the number of non-zero entries (ap[n])  
  */  
 int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])  
 {  
     int *sz = Alloca(n, int), j;  
     R_CheckStack();  
   
     for (j = n - 1; j >= 0; j--) {  
         int parent = pr[j];  
         sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);  
     }  
     ap[0] = 0;  
     for (j = 0; j < n; j++)  
         ap[j+1] = ap[j] + sz[j];  
     return ap[n];  
 }  
   
 /**  
  * Derive the row index array for the inverse of L from the parent array  
  *  
  * @param n length of parent array  
  * @param countDiag 0 for a unit triangular matrix with implicit diagonal, otherwise 1  
  * @param pr parent vector describing the elimination tree  
  * @param ai row index vector of length ap[n]  
  */  
 void parent_inv_ai(int n, int countDiag, const int pr[], int ai[])  
 {  
     int j, k, pos = 0;  
     for (j = 0; j < n; j++) {  
         if (countDiag) ai[pos++] = j;  
         for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;  
     }  
 }  
   
 SEXP Parent_inverse(SEXP par, SEXP unitdiag)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));  
     int *ap, *ai, *dims, *pr = INTEGER(par),  
         countDiag = 1 - asLogical(unitdiag),  
         j, n = length(par), nnz;  
     double *ax;  
   
     if (!isInteger(par)) error(_("par argument must be an integer vector"));  
     SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));  
     ap = INTEGER(GET_SLOT(ans, Matrix_pSym));  
     nnz = parent_inv_ap(n, countDiag, pr, ap);  
     SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));  
     ai = INTEGER(GET_SLOT(ans, Matrix_iSym));  
     SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));  
     ax = REAL(GET_SLOT(ans, Matrix_xSym));  
     for (j = 0; j < nnz; j++) ax[j] = 1.;  
     SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));  
     dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));  
     dims[0] = dims[1] = n;  
     SET_SLOT(ans, Matrix_uploSym, mkString("L"));  
     SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));  
     parent_inv_ai(n, countDiag, pr, ai);  
     UNPROTECT(1);  
     return ans;  
 }  
   
79  SEXP dtCMatrix_solve(SEXP a)  SEXP dtCMatrix_solve(SEXP a)
80  {  {
81      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
82      CSP A = AS_CSP(a);      CSP A = AS_CSP(Csparse_diagU2N(a));
83      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
84          bnz = 10 * A->n,        /* initial estimate of nnz in b */          bnz = 10 * A->n,        /* initial estimate of nnz in b */
85          lo = uplo_P(a)[0] == 'L', top;          lo = uplo_P(a)[0] == 'L', top;
# Line 204  Line 134 
134  {  {
135      int cl = asLogical(classed);      int cl = asLogical(classed);
136      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
137      CSP A = AS_CSP(a);      CSP A = AS_CSP(Csparse_diagU2N(a));
138      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
139          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
140                           getAttrib(b, R_DimSymbol));                           getAttrib(b, R_DimSymbol));
# Line 224  Line 154 
154      return ans;      return ans;
155  }  }
156    
157  SEXP dtCMatrix_upper_solve(SEXP a)  SEXP dtCMatrix_sparse_solve(SEXP a, SEXP b)
158  {  {
159      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
160      int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',      CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));
161          n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
162          *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),          xnz = 10 * B->p[B->n],  /* initial estimate of nnz in x */
163          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          lo = uplo_P(a)[0] == 'L', top;
164          *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));      /* These arrays must use Calloc because of possible Realloc */
165      int bnz = 10 * ap[n];         /* initial estimate of nnz in b */      int *ti = Calloc(xnz, int), p, j, nz, pos = 0;
166      int *ti = Alloca(bnz, int), j, nz;      double *tx = Calloc(xnz, double);
167      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Alloca(bnz, double),      double  *wrk = Alloca(A->n, double);
168          *tmp = Alloca(n, double);      int *xi = Alloca(2*A->n, int);      /* for cs_reach */
169      R_CheckStack();      R_CheckStack();
170    
171      if (lo || (!unit))      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
172          error(_("Code written for unit upper triangular unit matrices"));          error(_("Dimensions of system to be solved are inconsistent"));
173      bp[0] = 0;      slot_dup(ans, b, Matrix_DimSym);
174      for (j = 0; j < n; j++) {      SET_DimNames(ans, b);
175          int i, i1 = ap[j + 1];      xp[0] = 0;
176          AZERO(tmp, n);      for (j = 0; j < B->n; j++) {
177          for (i = ap[j]; i < i1; i++) {          /* (wrk[top:n],xi[top:n]) :=  A^{-1} B[,j] */
178              int ii = ai[i], k;          top = cs_spsolve (A, B, j, xi, wrk, 0, lo);
179              if (unit) tmp[ii] -= ax[i];          nz = A->n - top;
180              for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];          xp[j + 1] = nz + xp[j];
181            if (xp[j + 1] > xnz) {
182                while (xp[j + 1] > xnz) xnz *= 2;
183                ti = Realloc(ti, xnz, int);
184                tx = Realloc(tx, xnz, double);
185          }          }
186          for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;          if (lo)
187          bp[j + 1] = bp[j] + nz;              for(p = top; p < A->n; p++, pos++) {
188          if (bp[j + 1] > bnz) {                  ti[pos] = xi[p];
189              while (bp[j + 1] > bnz) bnz *= 2;                  tx[pos] = wrk[xi[p]];
             ti = Realloc(ti, bnz, int);  
             tx = Realloc(tx, bnz, double);  
190          }          }
191          i1 = bp[j];          else /* upper triagonal */
192          for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}              for(p = A->n - 1; p >= top; p--, pos++) {
193                    ti[pos] = xi[p];
194                    tx[pos] = wrk[xi[p]];
195      }      }
196      nz = bp[n];      }
197        nz = xp[A->n];
198      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
199      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
200      slot_dup(ans, a, Matrix_DimSym);  
201      SET_DimNames(ans, a);      Free(ti); Free(tx);
     slot_dup(ans, a, Matrix_uploSym);  
     slot_dup(ans, a, Matrix_diagSym);  
202      UNPROTECT(1);      UNPROTECT(1);
203      return ans;      return ans;
204  }  }

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

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