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 2210, Wed Jul 16 22:20:22 2008 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_);
6    
7    /* This is used for *BOTH* triangular and symmetric Csparse: */
8    SEXP tCMatrix_validate(SEXP x)
9  {  {
10      return triangularMatrix_validate(x);      SEXP val = xCMatrix_validate(x);/* checks x slot */
11      /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */      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            /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
25            if(uploT) {
26                for (k = 0; k < nnz; k++)
27                    if(xi[k] > xj[k]) {
28                        RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
29                    }
30            }
31            else {
32                for (k = 0; k < nnz; k++)
33                    if(xi[k] < xj[k]) {
34                        RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
35                    }
36  }  }
37    
38  SEXP tsc_transpose(SEXP x)          RETURN(ScalarLogical(1));
39  {      }
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix"))),  
         islot = GET_SLOT(x, Matrix_iSym);  
     int nnz = length(islot),  
         *adims, *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));  
     int up = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U';  
   
     adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  
     adims[0] = xdims[1]; adims[1] = xdims[0];  
   
     if(diag_value(x) == 'U')  
         SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(x, Matrix_diagSym)));  
     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;  
40  }  }
41    
42  SEXP tsc_to_dgTMatrix(SEXP x)  /* 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 (CHAR(STRING_ELT(GET_SLOT(x, Matrix_diagSym), 0))[0] != '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)));
56              *p = INTEGER(pslot);  
57          double *ax;          expand_cmprPt(length(pslot) - 1, INTEGER(pslot), xi);
58    
59          ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix")));          /* Maybe FIXME: ">" should be ">="      for diag = 'U' (uplo = 'U') */
60          SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));          if(uploT) {
61          SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));              for (k = 0; k < nnz; k++)
62          ai = INTEGER(GET_SLOT(ans, Matrix_iSym));                  if(xi[k] > xj[k]) {
63          Memcpy(ai, INTEGER(islot), nod);                      RETURN(mkString(_("uplo='U' must not have sparse entries below the diagonal")));
         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;  
64          }          }
         UNPROTECT(1);  
65      }      }
66      return ans;          else {
67                for (k = 0; k < nnz; k++)
68                    if(xi[k] < xj[k]) {
69                        RETURN(mkString(_("uplo='L' must not have sparse entries above the diagonal")));
70                    }
71  }  }
72    
73  /**          RETURN(ScalarLogical(1));
74   * Derive the column pointer vector for the inverse of L from the parent array      }
75   *  }
  * @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;  
76    
77      for (j = n - 1; j >= 0; j--) {  #undef RETURN
78          int parent = pr[j];  
79          sz[j] = (parent < 0) ?  countDiag : (1 + sz[parent]);  SEXP dtCMatrix_solve(SEXP a)
     }  
     ap[0] = 0;  
     for (j = 0; j < n; j++)  
         ap[j+1] = ap[j] + sz[j];  
     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[])  
80  {  {
81      int j, k, pos = 0;      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
82      for (j = 0; j < n; j++) {      CSP A = AS_CSP(Csparse_diagU2N(a));
83          if (countDiag) ai[pos++] = j;      int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
84          for (k = pr[j]; k >= 0; k = pr[k]) ai[pos++] = k;          bnz = 10 * A->n,        /* initial estimate of nnz in b */
85            lo = uplo_P(a)[0] == 'L', top;
86        /* These arrays must use Calloc because of possible Realloc */
87        int *ti = Calloc(bnz, int), p, j, nz, pos = 0;
88        double *tx = Calloc(bnz, double);
89        cs *u = cs_spalloc(A->n, 1,1,1,0);  /* Sparse unit vector */
90        double  *wrk = Alloca(A->n, double);
91        int *xi = Alloca(2*A->n, int);      /* for cs_reach */
92        R_CheckStack();
93    
94        slot_dup(ans, a, Matrix_DimSym);
95        SET_DimNames(ans, a);
96        slot_dup(ans, a, Matrix_uploSym);
97        slot_dup(ans, a, Matrix_diagSym);
98        /* initialize the "constant part" of the sparse unit vector */
99        u->x[0] = 1.;
100        u->p[0] = 0; u->p[1] = 1;
101        bp[0] = 0;
102        for (j = 0; j < A->n; j++) {
103            u->i[0] = j;                    /* u := j'th unit vector */
104            /* (wrk[top:n],xi[top:n]) :=  A^{-1} u  : */
105            top = cs_spsolve (A, u, 0, xi, wrk, 0, lo);
106            nz = A->n - top;
107            bp[j + 1] = nz + bp[j];
108            if (bp[j + 1] > bnz) {
109                while (bp[j + 1] > bnz) bnz *= 2;
110                ti = Realloc(ti, bnz, int);
111                tx = Realloc(tx, bnz, double);
112            }
113            if (lo)
114                for(p = top; p < A->n; p++, pos++) {
115                    ti[pos] = xi[p];
116                    tx[pos] = wrk[xi[p]];
117                }
118            else /* upper triagonal */
119                for(p = A->n - 1; p >= top; p--, pos++) {
120                    ti[pos] = xi[p];
121                    tx[pos] = wrk[xi[p]];
122                }
123        }
124        nz = bp[A->n];
125        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
126        Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
127    
128        Free(ti); Free(tx); cs_spfree(u);
129        UNPROTECT(1);
130        return ans;
131      }      }
132    
133    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
134    {
135        int cl = asLogical(classed);
136        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
137        CSP A = AS_CSP(Csparse_diagU2N(a));
138        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
139            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
140                             getAttrib(b, R_DimSymbol));
141        int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
142        double *bx;
143        R_CheckStack();
144    
145        if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
146            error(_("Dimensions of system to be solved are inconsistent"));
147        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
148        /* FIXME: copy dimnames or Dimnames as well */
149        bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
150                    REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
151        for (j = 0; j < nrhs; j++)
152            lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
153        UNPROTECT(1);
154        return ans;
155  }  }
156    
157  SEXP Parent_inverse(SEXP par, SEXP unitdiag)  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 *ap, *ai, *dims, *pr = INTEGER(par),      CSP A = AS_CSP(Csparse_diagU2N(a)), B = AS_CSP(Csparse_diagU2N(b));
161          countDiag = 1 - asLogical(unitdiag),      int *xp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (B->n) + 1)),
162          j, n = length(par), nnz;          xnz = 10 * B->p[B->n],  /* initial estimate of nnz in x */
163      double *ax;          lo = uplo_P(a)[0] == 'L', top;
164        /* These arrays must use Calloc because of possible Realloc */
165      if (!isInteger(par)) error(_("par argument must be an integer vector"));      int *ti = Calloc(xnz, int), p, j, nz, pos = 0;
166      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, n + 1));      double *tx = Calloc(xnz, double);
167      ap = INTEGER(GET_SLOT(ans, Matrix_pSym));      double  *wrk = Alloca(A->n, double);
168      nnz = parent_inv_ap(n, countDiag, pr, ap);      int *xi = Alloca(2*A->n, int);      /* for cs_reach */
169      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));      R_CheckStack();
170      ai = INTEGER(GET_SLOT(ans, Matrix_iSym));  
171      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      if (A->m != A->n || B->n < 1 || A->n < 1 || A->n != B->m)
172      ax = REAL(GET_SLOT(ans, Matrix_xSym));          error(_("Dimensions of system to be solved are inconsistent"));
173      for (j = 0; j < nnz; j++) ax[j] = 1.;      slot_dup(ans, b, Matrix_DimSym);
174      SET_SLOT(ans, Matrix_DimSym, allocVector(INTSXP, 2));      SET_DimNames(ans, b);
175      dims = INTEGER(GET_SLOT(ans, Matrix_DimSym));      xp[0] = 0;
176      dims[0] = dims[1] = n;      for (j = 0; j < B->n; j++) {
177      SET_SLOT(ans, Matrix_uploSym, mkString("L"));          /* (wrk[top:n],xi[top:n]) :=  A^{-1} B[,j] */
178      SET_SLOT(ans, Matrix_diagSym, (countDiag ? mkString("N") : mkString("U")));          top = cs_spsolve (A, B, j, xi, wrk, 0, lo);
179      parent_inv_ai(n, countDiag, pr, ai);          nz = A->n - top;
180            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            if (lo)
187                for(p = top; p < A->n; p++, pos++) {
188                    ti[pos] = xi[p];
189                    tx[pos] = wrk[xi[p]];
190                }
191            else /* upper triagonal */
192                for(p = A->n - 1; p >= top; p--, pos++) {
193                    ti[pos] = xi[p];
194                    tx[pos] = wrk[xi[p]];
195                }
196        }
197        nz = xp[A->n];
198        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP,  nz)), ti, nz);
199        Memcpy(   REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
200    
201        Free(ti); Free(tx);
202      UNPROTECT(1);      UNPROTECT(1);
203      return ans;      return ans;
204  }  }

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

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