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 677, Mon Mar 28 14:03:32 2005 UTC revision 1710, Tue Dec 26 15:57:06 2006 UTC
# Line 1  Line 1 
1                                  /* Sparse triangular 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)
     SEXP val;  
   
     if (isString(val = check_scalar_string(GET_SLOT(x, Matrix_uploSym),  
                                            "LU", "uplo"))) return val;  
     if (isString(val = check_scalar_string(GET_SLOT(x, Matrix_diagSym),  
                                            "NU", "diag"))) return val;  
     return ScalarLogical(1);  
 }  
   
 SEXP tsc_transpose(SEXP x)  
7  {  {
8        SEXP val = xCMatrix_validate(x);/* checks x slot */
9        if(isString(val))
10            return(val);
11        else {
12      SEXP      SEXP
13          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),              islot = GET_SLOT(x, Matrix_iSym),
         islot = GET_SLOT(x, Matrix_iSym);  
     int nnz = length(islot),  
         *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));  
   
     adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  
     adims[0] = xdims[1]; adims[1] = xdims[0];  
     if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')  
         SET_SLOT(ans, Matrix_uploSym, mkString("L"));  
     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;  
 }  
   
 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),  
14              pslot = GET_SLOT(x, Matrix_pSym);              pslot = GET_SLOT(x, Matrix_pSym);
15          int *ai, *aj, j,          int uploT = (*uplo_P(x) == 'U'),
16              n = length(pslot) - 1,              k, nnz = length(islot),
17              nod = length(islot),              *xi = INTEGER(islot),
18              nout = n + nod,              *xj = INTEGER(PROTECT(allocVector(INTSXP, nnz)));
             *p = INTEGER(pslot);  
         double *ax;  
19    
20          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));  #define RETURN(_CH_)   UNPROTECT(1); return (_CH_);
21          SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));  
22          SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xj);
23          ai = INTEGER(GET_SLOT(ans, Matrix_iSym));  
24          Memcpy(ai, INTEGER(islot), nod);          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
25          SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));          if(uploT) {
26          aj = INTEGER(GET_SLOT(ans, Matrix_jSym));              for (k = 0; k < nnz; k++)
27          SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));                  if(xi[k] > xj[k]) {
28          ax = REAL(GET_SLOT(ans, Matrix_xSym));                      RETURN(mkString(_("uplo='U' must not have sparse entries in lower diagonal")));
         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;  
29          }          }
         UNPROTECT(1);  
30      }      }
31      return ans;          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        }
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
45   *   *
# Line 138  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        cs *A = Matrix_as_cs(a);
116        int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
117            lo = uplo_P(a)[0] == 'L',
118            bnz = 10 * A->n;        /* initial estimate of nnz in b */
119        int *ti = Calloc(bnz, int), i, j, nz, pos = 0;
120        double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
121    
122        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
123        SET_SLOT(ans, Matrix_DimNamesSym,
124                 duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
125        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
126        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
127        bp[0] = 0;
128        for (j = 0; j < A->n; j++) {
129            AZERO(wrk, A->n);
130            wrk[j] = 1;
131            lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);
132            for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;
133            bp[j + 1] = nz + bp[j];
134            if (bp[j + 1] > bnz) {
135                while (bp[j + 1] > bnz) bnz *= 2;
136                ti = Realloc(ti, bnz, int);
137                tx = Realloc(tx, bnz, double);
138            }
139            for (i = 0; i < A->n; i++)
140                if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}
141        }
142        nz = bp[A->n];
143        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
144        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
145    
146        Free(A); Free(ti); Free(tx);
147        UNPROTECT(1);
148        return ans;
149    }
150    
151    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
152    {
153        int cl = asLogical(classed);
154        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
155        cs *A = Matrix_as_cs(a);
156        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
157            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
158                             getAttrib(b, R_DimSymbol));
159        int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
160        double *bx;
161    
162        if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
163            error(_("Dimensions of system to be solved are inconsistent"));
164        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
165        /* copy dimnames or Dimnames as well */
166        bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
167                    REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
168        for (j = 0; j < nrhs; j++)
169            lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
170        Free(A);
171        UNPROTECT(1);
172        return ans;
173    }
174    
175    SEXP dtCMatrix_upper_solve(SEXP a)
176    {
177        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
178        int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
179            n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
180            *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
181            *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
182            *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
183        int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
184        int *ti = Calloc(bnz, int), j, nz;
185        double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),
186            *tmp = Calloc(n, double);
187    
188        if (lo || (!unit))
189            error(_("Code written for unit upper triangular unit matrices"));
190        bp[0] = 0;
191        for (j = 0; j < n; j++) {
192            int i, i1 = ap[j + 1];
193            AZERO(tmp, n);
194            for (i = ap[j]; i < i1; i++) {
195                int ii = ai[i], k;
196                if (unit) tmp[ii] -= ax[i];
197                for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
198            }
199            for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
200            bp[j + 1] = bp[j] + nz;
201            if (bp[j + 1] > bnz) {
202                while (bp[j + 1] > bnz) bnz *= 2;
203                ti = Realloc(ti, bnz, int);
204                tx = Realloc(tx, bnz, double);
205            }
206            i1 = bp[j];
207            for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
208        }
209        nz = bp[n];
210        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
211        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
212        Free(tmp); Free(tx); Free(ti);
213        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
214        SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
215        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
216        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
217        UNPROTECT(1);
218        return ans;
219    }

Legend:
Removed from v.677  
changed lines
  Added in v.1710

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