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 1248, Thu Apr 13 22:05:22 2006 UTC pkg/Matrix/src/dtCMatrix.c revision 2889, Thu Aug 8 21:06:22 2013 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)  #define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
 {  
     return triangularMatrix_validate(x);  
     /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */  
 }  
6    
7  #if 0  /* This is used for *BOTH* triangular and symmetric Csparse: */
8  SEXP tsc_transpose(SEXP x)  SEXP tCMatrix_validate(SEXP x)
9  {  {
10      cholmod_sparse *cx = as_cholmod_sparse(x);      SEXP val = xCMatrix_validate(x);/* checks x slot */
11        if(isString(val))
12            return(val);
13        else {
14            SEXP
15                islot = GET_SLOT(x, Matrix_iSym),
16                pslot = GET_SLOT(x, Matrix_pSym);
17            int uploT = (*uplo_P(x) == 'U'),
18                k, nnz = length(islot),
19                *xi = INTEGER(islot),
20                *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
21    
22            expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23    
24      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
25          islot = GET_SLOT(x, Matrix_iSym);          if(uploT) {
26      int nnz = length(islot),              for (k = 0; k < nnz; k++)
27          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));                  if(xi[k] > xj[k]) {
28      int up = uplo_P(x)[0] == 'U';                      RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
29                    }
30      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));          }
31      adims[0] = xdims[1]; adims[1] = xdims[0];          else {
32                for (k = 0; k < nnz; k++)
33      if(*diag_P(x) == 'U')                  if(xi[k] < xj[k]) {
34          SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));                      RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
35      SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));                  }
   
     csc_compTr(xdims[0], xdims[1], nnz,  
                INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),  
                REAL(GET_SLOT(x, Matrix_xSym)),  
                INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),  
                INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),  
                REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));  
     UNPROTECT(1);  
     return ans;  
36  }  }
 #endif  
37    
38  SEXP tsc_to_dgTMatrix(SEXP x)          RETURN(ScalarLogical(1));
39        }
40    }
41    
42    /* This is used for *BOTH* triangular and symmetric Rsparse: */
43    SEXP tRMatrix_validate(SEXP x)
44  {  {
45      SEXP ans;      SEXP val = xRMatrix_validate(x);/* checks x slot */
46      if (*diag_P(x) != 'U')      if(isString(val))
47          ans = compressed_to_dgTMatrix(x, ScalarLogical(1));          return(val);
48      else {                      /* unit triangular matrix */      else {
49          SEXP islot = GET_SLOT(x, Matrix_iSym),          SEXP
50                jslot = GET_SLOT(x, Matrix_jSym),
51              pslot = GET_SLOT(x, Matrix_pSym);              pslot = GET_SLOT(x, Matrix_pSym);
52          int *ai, *aj, j,          int uploT = (*uplo_P(x) == 'U'),
53              n = length(pslot) - 1,              k, nnz = length(jslot),
54              nod = length(islot),              *xj = INTEGER(jslot),
55              nout = n + nod,              *xi = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
             *p = INTEGER(pslot);  
         double *ax;  
   
         ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));  
         SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));  
         SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));  
         ai = INTEGER(GET_SLOT(ans, Matrix_iSym));  
         Memcpy(ai, INTEGER(islot), nod);  
         SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));  
         aj = INTEGER(GET_SLOT(ans, Matrix_jSym));  
         SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));  
         ax = REAL(GET_SLOT(ans, Matrix_xSym));  
         Memcpy(ax, REAL(GET_SLOT(x, Matrix_xSym)), nod);  
         for (j = 0; j < n; j++) {  
             int jj, npj = nod + j,  p2 = p[j+1];  
             aj[npj] = ai[npj] = j;  
             ax[npj] = 1.;  
             for (jj = p[j]; jj < p2; jj++) aj[jj] = j;  
         }  
         UNPROTECT(1);  
     }  
     return ans;  
 }  
   
 /**  
  * 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 = Calloc(n, int), j;  
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")));
64      Free(sz);                  }
65      return ap[n];          }
66  }          else {
67                for (k = 0; k < nnz; k++)
68  /**                  if(xi[k] < xj[k]) {
69   * Derive the row index array for the inverse of L from the parent array                      RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
  *  
  * @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;  
70      }      }
71  }  }
72    
73  SEXP Parent_inverse(SEXP par, SEXP unitdiag)          RETURN(ScalarLogical(1));
74        }
75    }
76    
77    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
78  {  {
79      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));      int cl = asLogical(classed);
80      int *ap, *ai, *dims, *pr = INTEGER(par),      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
81          countDiag = 1 - asLogical(unitdiag),      CSP A = AS_CSP(a);
82          j, n = length(par), nnz;      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
83      double *ax;          *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
84                             getAttrib(b, R_DimSymbol));
85      if (!isInteger(par)) error(_("par argument must be an integer vector"));      int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
86      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));      double *bx;
87      ap = INTEGER(GET_SLOT(ans, Matrix_pSym));      R_CheckStack();
88      nnz = parent_inv_ap(n, countDiag, pr, ap);  
89      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
90      ai = INTEGER(GET_SLOT(ans, Matrix_iSym));          error(_("Dimensions of system to be solved are inconsistent"));
91      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
92      ax = REAL(GET_SLOT(ans, Matrix_xSym));      /* FIXME: copy dimnames or Dimnames as well */
93      for (j = 0; j < nnz; j++) ax[j] = 1.;      bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
94      SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));                  REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
95      dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));      for (j = 0; j < nrhs; j++)
96      dims[0] = dims[1] = n;          lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
97      SET_SLOT(ans, Matrix_uploSym, mkString("L"));      RETURN(ans);
     SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));  
     parent_inv_ai(n, countDiag, pr, ai);  
     UNPROTECT(1);  
     return ans;  
98  }  }
99    
100  SEXP dtCMatrix_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], nnz,      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)), *bi,          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 *ri = Calloc(bnz, int), *ind = Calloc(n, int), j;      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
110      double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *bx;          xnz = 10 * B->p[B->n];  /* initial estimate of nnz in x */
111        int k, lo = uplo_P(a)[0] == 'L', pos = 0;
112      SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));      int    *ti = Calloc(xnz, int),     *xi = Calloc(2*A->n, int); /* for cs_reach */
113      SET_SLOT(ans, Matrix_DimNamesSym,      double *tx = Calloc(xnz, double), *wrk = Calloc(  A->n, double);
114               duplicate(GET_SLOT(a, Matrix_DimNamesSym)));  
115      SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));      slot_dup(ans, b, Matrix_DimSym);
116      SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));      SET_DimNames(ans, b);
117        xp[0] = 0;
118      if (!(lo && unit))      for (k = 0; k < B->n; k++) {
119          error("code for non-unit or upper triangular not yet written");          int top = cs_spsolve (A, B, k, xi, wrk, (int *)NULL, lo);
120      /* Initially bp will contain increasing negative values ending at zero. */          int nz = A->n - top;
121      /* Later we add the negative of bp[0] to all values. */  
122      bp[n] = 0;          xp[k + 1] = nz + xp[k];
123      for (j = n - 1; j >= 0; j--) { /* columns in reverse order */          if (xp[k + 1] > xnz) {
124          int i, i1 = ap[j], i2 = ap[j + 1], k, nr;              while (xp[k + 1] > xnz) xnz *= 2;
125          if (i1 < i2) AZERO(ind, n);              ti = Realloc(ti, xnz, int);
126          for (i = i1; i < i2; i++) {              tx = Realloc(tx, xnz, double);
             ind[ai[i]] = 1;  
             for (k = -bp[ai[i] + 1]; k < -bp[ai[i]]; k++) ind[ri[k]] = 1;  
         }  
         for (k = 0, nr = 0; k < n; k++) if (ind[k]) nr++;  
         if ((nr - bp[j + 1]) > bnz) {  
             while (nr > (bnz + bp[j + 1])) bnz *= 2;  
             ri = Realloc(ri, bnz, int);  
         }  
         bp[j] = bp[j + 1] - nr;  
         for (k = 0, i = -bp[j + 1]; k < n; k++) if (ind[k]) ri[i++] = k;  
     }  
     bnz = -bp[0];  
     bi = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, bnz));  
     bx = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, bnz));  
     for (j = 0; j < n; j++) {  
         int bpnew = bp[j] + bnz;  
         Memcpy(bi + bpnew, ri - bp[j], bp[j + 1] - bp[j]);  
         bp[j] = bpnew;  
     }  
     /* insert code for calculating the actual values here */  
     for (j = 0; j < bnz; j++) bx[j] = NA_REAL;  
   
     Free(ind); Free(ri);  
     UNPROTECT(1);  
     return ans;  
127  }  }
128            if (lo)                 /* increasing row order */
129                for(int p = top; p < A->n; p++, pos++) {
130                    ti[pos] = xi[p];
131                    tx[pos] = wrk[xi[p]];
132                }
133            else                    /* decreasing order, reverse copy */
134                for(int p = A->n - 1; p >= top; p--, pos++) {
135                    ti[pos] = xi[p];
136                    tx[pos] = wrk[xi[p]];
137                }
138        }
139        xnz = xp[B->n];
140        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  xnz)), ti, xnz);
141        Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, xnz)), tx, xnz);
142    
143        Free(ti);  Free(tx);
144        Free(wrk); Free(xi);
145    
146        RETURN(ans);
147    }
148    #undef RETURN
149    

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