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 945, Wed Sep 28 08:54:28 2005 UTC revision 1960, Fri Jul 6 16:54:43 2007 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)  /* This should be use for *BOTH* triangular and symmetric Csparse: */
6    SEXP tCMatrix_validate(SEXP x)
7  {  {
8      return triangularMatrix_validate(x);      SEXP val = xCMatrix_validate(x);/* checks x slot */
9      /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */      if(isString(val))
10  }          return(val);
11        else {
12  SEXP tsc_transpose(SEXP x)          SEXP
13  {              islot = GET_SLOT(x, Matrix_iSym),
14      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),              pslot = GET_SLOT(x, Matrix_pSym);
15          islot = GET_SLOT(x, Matrix_iSym);          int uploT = (*uplo_P(x) == 'U'),
16      int nnz = length(islot),              k, nnz = length(islot),
17          *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));              *xi = INTEGER(islot),
18      int up = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U';              *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
19    
20      adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  #define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
     adims[0] = xdims[1]; adims[1] = xdims[0];  
21    
22      if(diag_value(x) == 'U')          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
         SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));  
     SET_SLOT(ans, Matrix_uploSym, mkString(up ? "L" : "U"));  
23    
24      csc_compTr(xdims[0], xdims[1], nnz,          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
25                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),          if(uploT) {
26                 REAL(GET_SLOT(x, Matrix_xSym)),              for (k = 0; k < nnz; k++)
27                 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),                  if(xi[k] > xj[k]) {
28                 INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),                      RETURN(mkString(_("uplo='U' must not have sparse entries in lower diagonal")));
                REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));  
     UNPROTECT(1);  
     return ans;  
29  }  }
   
 SEXP tsc_to_dgTMatrix(SEXP x)  
 {  
     SEXP ans;  
     if (CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] != 'U')  
         ans = compressed_to_dgTMatrix(x, ScalarLogical(1));  
     else {                      /* unit triangular matrix */  
         SEXP islot = GET_SLOT(x, Matrix_iSym),  
             pslot = GET_SLOT(x, Matrix_pSym);  
         int *ai, *aj, j,  
             n = length(pslot) - 1,  
             nod = length(islot),  
             nout = n + nod,  
             *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;  
30          }          }
31          UNPROTECT(1);          else {
32                for (k = 0; k < nnz; k++)
33                    if(xi[k] < xj[k]) {
34                        RETURN(mkString(_("uplo='L' must not have sparse entries in upper diagonal")));
35                    }
36            }
37    
38            RETURN(ScalarLogical(1));
39      }      }
     return ans;  
40  }  }
41    #undef RETURN
42    
43  /**  /**
44   * Derive the column pointer vector for the inverse of L from the parent array   * Derive the column pointer vector for the inverse of L from the parent array
# Line 80  Line 52 
52   */   */
53  int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])  int parent_inv_ap(int n, int countDiag, const int pr[], int ap[])
54  {  {
55      int *sz = Calloc(n, int), j;      int *sz = Alloca(n, int), j;
56        R_CheckStack();
57    
58      for (j = n - 1; j >= 0; j--) {      for (j = n - 1; j >= 0; j--) {
59          int parent = pr[j];          int parent = pr[j];
# Line 89  Line 62 
62      ap[0] = 0;      ap[0] = 0;
63      for (j = 0; j < n; j++)      for (j = 0; j < n; j++)
64          ap[j+1] = ap[j] + sz[j];          ap[j+1] = ap[j] + sz[j];
     Free(sz);  
65      return ap[n];      return ap[n];
66  }  }
67    
# Line 136  Line 108 
108      UNPROTECT(1);      UNPROTECT(1);
109      return ans;      return ans;
110  }  }
111    
112    SEXP dtCMatrix_solve(SEXP a)
113    {
114        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
115        CSP A = AS_CSP(a);
116        int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
117            bnz = 10 * A->n,        /* initial estimate of nnz in b */
118            lo = uplo_P(a)[0] == 'L', top;
119        /* These arrays must use Calloc because of possible Realloc */
120        int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
121        double *tx = Calloc(bnz, double);
122        cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */
123        double  *wrk = Alloca(A->n, double);
124        int *xi = Alloca(2*A->n, int);      /* for cs_reach */
125        R_CheckStack();
126    
127        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
128        SET_DimNames(ans, a);
129        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
130        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
131        /* initialize the "constant part" of the sparse unit vector */
132        u->x[0] = 1.;
133        u->p[0] = 0; u->p[1] = 1;
134        bp[0] = 0;
135        for (j = 0; j < A->n; j++) {
136            u->i[0] = j;                    /* u := j'th unit vector */
137            /* (wrk[top:n],xi[top:n]) :=  A^{-1} u  : */
138            top = cs_spsolve (A, u, 0, xi, wrk, 0, lo);
139            nz = A->n - top;
140            bp[j + 1] = nz + bp[j];
141            if (bp[j + 1] > bnz) {
142                while (bp[j + 1] > bnz) bnz *= 2;
143                ti = Realloc(ti, bnz, int);
144                tx = Realloc(tx, bnz, double);
145            }
146            if (lo)
147                for(p = top; p < A->n; p++, pos++) {
148                    ti[pos] = xi[p];
149                    tx[pos] = wrk[xi[p]];
150                }
151            else /* upper triagonal */
152                for(p = A->n - 1; p >= top; p--, pos++) {
153                    ti[pos] = xi[p];
154                    tx[pos] = wrk[xi[p]];
155                }
156        }
157        nz = bp[A->n];
158        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
159        Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
160    
161        Free(ti); Free(tx); cs_spfree(u);
162        UNPROTECT(1);
163        return ans;
164    }
165    
166    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
167    {
168        int cl = asLogical(classed);
169        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
170        CSP A = AS_CSP(a);
171        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
172            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
173                             getAttrib(b, R_DimSymbol));
174        int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
175        double *bx;
176        R_CheckStack();
177    
178        if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
179            error(_("Dimensions of system to be solved are inconsistent"));
180        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
181        /* FIXME: copy dimnames or Dimnames as well */
182        bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
183                    REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
184        for (j = 0; j < nrhs; j++)
185            lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
186        UNPROTECT(1);
187        return ans;
188    }
189    
190    SEXP dtCMatrix_upper_solve(SEXP a)
191    {
192        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
193        int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
194            n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
195            *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
196            *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
197            *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
198        int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
199        int *ti = Alloca(bnz, int), j, nz;
200        double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Alloca(bnz, double),
201            *tmp = Alloca(n, double);
202        R_CheckStack();
203    
204        if (lo || (!unit))
205            error(_("Code written for unit upper triangular unit matrices"));
206        bp[0] = 0;
207        for (j = 0; j < n; j++) {
208            int i, i1 = ap[j + 1];
209            AZERO(tmp, n);
210            for (i = ap[j]; i < i1; i++) {
211                int ii = ai[i], k;
212                if (unit) tmp[ii] -= ax[i];
213                for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
214            }
215            for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
216            bp[j + 1] = bp[j] + nz;
217            if (bp[j + 1] > bnz) {
218                while (bp[j + 1] > bnz) bnz *= 2;
219                ti = Realloc(ti, bnz, int);
220                tx = Realloc(tx, bnz, double);
221            }
222            i1 = bp[j];
223            for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
224        }
225        nz = bp[n];
226        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
227        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
228        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
229        SET_DimNames(ans, a);
230        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
231        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
232        UNPROTECT(1);
233        return ans;
234    }

Legend:
Removed from v.945  
changed lines
  Added in v.1960

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