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 1710, Tue Dec 26 15:57:06 2006 UTC pkg/Matrix/src/dtCMatrix.c revision 2889, Thu Aug 8 21:06:22 2013 UTC
# Line 2  Line 2 
2  #include "dtCMatrix.h"  #include "dtCMatrix.h"
3  #include "cs_utils.h"  #include "cs_utils.h"
4    
5  /* This should be use for *BOTH* triangular and symmetric Csparse: */  #define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
6    
7    /* This is used for *BOTH* triangular and symmetric Csparse: */
8  SEXP tCMatrix_validate(SEXP x)  SEXP tCMatrix_validate(SEXP x)
9  {  {
10      SEXP val = xCMatrix_validate(x);/* checks x slot */      SEXP val = xCMatrix_validate(x);/* checks x slot */
# Line 17  Line 19 
19              *xi = INTEGER(islot),              *xi = INTEGER(islot),
20              *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));              *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
21    
 #define RETURN(_CH_)   UNPROTECT(1); return (_CH_);  
   
22          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23    
24          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
25          if(uploT) {          if(uploT) {
26              for (k = 0; k < nnz; k++)              for (k = 0; k < nnz; k++)
27                  if(xi[k] > xj[k]) {                  if(xi[k] > xj[k]) {
28                      RETURN(mkString(_("uplo='U' must not have sparse entries in lower diagonal")));                      RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
29                  }                  }
30          }          }
31          else {          else {
32              for (k = 0; k < nnz; k++)              for (k = 0; k < nnz; k++)
33                  if(xi[k] < xj[k]) {                  if(xi[k] < xj[k]) {
34                      RETURN(mkString(_("uplo='L' must not have sparse entries in upper diagonal")));                      RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
35                  }                  }
36          }          }
37    
38          RETURN(ScalarLogical(1));          RETURN(ScalarLogical(1));
39      }      }
40  }  }
 #undef RETURN  
41    
42  /**  /* This is used for *BOTH* triangular and symmetric Rsparse: */
43   * Derive the column pointer vector for the inverse of L from the parent array  SEXP tRMatrix_validate(SEXP x)
  *  
  * @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[])  
44  {  {
45      int *sz = Calloc(n, int), j;      SEXP val = xRMatrix_validate(x);/* checks x slot */
46        if(isString(val))
47            return(val);
48        else {
49            SEXP
50                jslot = GET_SLOT(x, Matrix_jSym),
51                pslot = GET_SLOT(x, Matrix_pSym);
52            int uploT = (*uplo_P(x) == 'U'),
53                k, nnz = length(jslot),
54                *xj = INTEGER(jslot),
55                *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
56    
57      for (j = n - 1; j >= 0; j--) {          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
58          int parent = pr[j];  
59          sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
60      }          if(uploT) {
61      ap[0] = 0;              for (k = 0; k < nnz; k++)
62      for (j = 0; j < n; j++)                  if(xi[k] > xj[k]) {
63          ap[j+1] = ap[j] + sz[j];                      RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
     Free(sz);  
     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;  
64      }      }
65  }  }
66            else {
67  SEXP Parent_inverse(SEXP par, SEXP unitdiag)              for (k = 0; k < nnz; k++)
68  {                  if(xi[k] < xj[k]) {
69      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));                      RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
70      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;  
71  }  }
72    
73  SEXP dtCMatrix_solve(SEXP a)          RETURN(ScalarLogical(1));
74  {      }
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));  
     cs *A = Matrix_as_cs(a);  
     int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),  
         lo = uplo_P(a)[0] == 'L',  
         bnz = 10 * A->n;        /* initial estimate of nnz in b */  
     int *ti = Calloc(bnz, int), i, j, nz, pos = 0;  
     double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);  
   
     SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));  
     SET_SLOT(ans, Matrix_DimNamesSym,  
              duplicate(GET_SLOT(a, Matrix_DimNamesSym)));  
     SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));  
     SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));  
     bp[0] = 0;  
     for (j = 0; j < A->n; j++) {  
         AZERO(wrk, A->n);  
         wrk[j] = 1;  
         lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);  
         for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;  
         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);  
         }  
         for (i = 0; i < A->n; i++)  
             if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}  
     }  
     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(A); Free(ti); Free(tx);  
     UNPROTECT(1);  
     return ans;  
75  }  }
76    
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      cs *A = Matrix_as_cs(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));
85      int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');      int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
86      double *bx;      double *bx;
87        R_CheckStack();
88    
89      if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])      if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
90          error(_("Dimensions of system to be solved are inconsistent"));          error(_("Dimensions of system to be solved are inconsistent"));
91      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
92      /* copy dimnames or Dimnames as well */      /* FIXME: copy dimnames or Dimnames as well */
93      bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),      bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
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      Free(A);      RETURN(ans);
     UNPROTECT(1);  
     return ans;  
98  }  }
99    
100  SEXP dtCMatrix_upper_solve(SEXP a)  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      int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',      CSP A = AS_CSP(a), B = AS_CSP(b);
104          n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],      R_CheckStack();
105          *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
106          *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),          error(_("Dimensions of system to be solved are inconsistent"));
107          *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));      // *before* Calloc()ing below [memory leak]!
108      int bnz = 10 * ap[n];         /* initial estimate of nnz in b */  
109      int *ti = Calloc(bnz, int), j, nz;      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
110      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),          xnz = 10 * B->p[B->n];  /* initial estimate of nnz in x */
111          *tmp = Calloc(n, double);      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      if (lo || (!unit))      double *tx = Calloc(xnz, double), *wrk = Calloc(  A->n, double);
114          error(_("Code written for unit upper triangular unit matrices"));  
115      bp[0] = 0;      slot_dup(ans, b, Matrix_DimSym);
116      for (j = 0; j < n; j++) {      SET_DimNames(ans, b);
117          int i, i1 = ap[j + 1];      xp[0] = 0;
118          AZERO(tmp, n);      for (k = 0; k < B->n; k++) {
119          for (i = ap[j]; i < i1; i++) {          int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
120              int ii = ai[i], k;          int nz = A->n - top;
121              if (unit) tmp[ii] -= ax[i];  
122              for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];          xp[k + 1] = nz + xp[k];
123          }          if (xp[k + 1] > xnz) {
124          for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;              while (xp[k + 1] > xnz) xnz *= 2;
125          bp[j + 1] = bp[j] + nz;              ti = Realloc(ti, xnz, int);
126          if (bp[j + 1] > bnz) {              tx = Realloc(tx, xnz, double);
127              while (bp[j + 1] > bnz) bnz *= 2;          }
128              ti = Realloc(ti, bnz, int);          if (lo)                 /* increasing row order */
129              tx = Realloc(tx, bnz, double);              for(int p = top; p < A->n; p++, pos++) {
130          }                  ti[pos] = xi[p];
131          i1 = bp[j];                  tx[pos] = wrk[xi[p]];
132          for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}              }
133      }          else                    /* decreasing order, reverse copy */
134      nz = bp[n];              for(int p = A->n - 1; p >= top; p--, pos++) {
135      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);                  ti[pos] = xi[p];
136      Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);                  tx[pos] = wrk[xi[p]];
137      Free(tmp); Free(tx); Free(ti);              }
138      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));      }
139      SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(a, Matrix_DimNamesSym)));      xnz = xp[B->n];
140      SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  xnz)), ti, xnz);
141      SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));      Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
142      UNPROTECT(1);  
143      return ans;      Free(ti);  Free(tx);
144        Free(wrk); Free(xi);
145    
146        RETURN(ans);
147  }  }
148    #undef RETURN
149    

Legend:
Removed from v.1710  
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