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 1661, Thu Nov 2 13:10:50 2006 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"
 SEXP tsc_validate(SEXP x)  
 {  
     return triangularMatrix_validate(x);  
     /* see ./dsCMatrix.c or ./dtpMatrix.c  on how to do more testing here */  
 }  
   
 SEXP tsc_transpose(SEXP x)  
 {  
     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;  
 }  
   
 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;  
         }  
         UNPROTECT(1);  
     }  
     return ans;  
 }  
4    
5  /**  /**
6   * 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 136  Line 70 
70      UNPROTECT(1);      UNPROTECT(1);
71      return ans;      return ans;
72  }  }
73    
74    SEXP dtCMatrix_solve(SEXP a)
75    {
76        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
77        cs *A = Matrix_as_cs(a);
78        int *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, (A->n) + 1)),
79            lo = uplo_P(a)[0] == 'L',
80            bnz = 10 * A->n;        /* initial estimate of nnz in b */
81        int *ti = Calloc(bnz, int), i, j, nz, pos = 0;
82        double *tx = Calloc(bnz, double), *wrk = Calloc(A->n, double);
83    
84        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
85        SET_SLOT(ans, Matrix_DimNamesSym,
86                 duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
87        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
88        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
89        bp[0] = 0;
90        for (j = 0; j < A->n; j++) {
91            AZERO(wrk, A->n);
92            wrk[j] = 1;
93            lo ? cs_lsolve(A, wrk) : cs_usolve(A, wrk);
94            for (i = 0, nz = 0; i < A->n; i++) if (wrk[i]) nz++;
95            bp[j + 1] = nz + bp[j];
96            if (bp[j + 1] > bnz) {
97                while (bp[j + 1] > bnz) bnz *= 2;
98                ti = Realloc(ti, bnz, int);
99                tx = Realloc(tx, bnz, double);
100            }
101            for (i = 0; i < A->n; i++)
102                if (wrk[i]) {ti[pos] = i; tx[pos] = wrk[i]; pos++;}
103        }
104        nz = bp[A->n];
105        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
106        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
107    
108        Free(A); Free(ti); Free(tx);
109        UNPROTECT(1);
110        return ans;
111    }
112    
113    SEXP dtCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
114    {
115        int cl = asLogical(classed);
116        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
117        cs *A = Matrix_as_cs(a);
118        int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
119            *bdims = INTEGER(cl ? GET_SLOT(b, Matrix_DimSym) :
120                             getAttrib(b, R_DimSymbol));
121        int j, n = bdims[0], nrhs = bdims[1], lo = (*uplo_P(a) == 'L');
122        double *bx;
123    
124        if (*adims != n || nrhs < 1 || *adims < 1 || *adims != adims[1])
125            error(_("Dimensions of system to be solved are inconsistent"));
126        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)), bdims, 2);
127        /* copy dimnames or Dimnames as well */
128        bx = Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * nrhs)),
129                    REAL(cl ? GET_SLOT(b, Matrix_xSym):b), n * nrhs);
130        for (j = 0; j < nrhs; j++)
131            lo ? cs_lsolve(A, bx + n * j) : cs_usolve(A, bx + n * j);
132        Free(A);
133        UNPROTECT(1);
134        return ans;
135    }
136    
137    SEXP dtCMatrix_upper_solve(SEXP a)
138    {
139        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
140        int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
141            n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0],
142            *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
143            *ap = INTEGER(GET_SLOT(a, Matrix_pSym)),
144            *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
145        int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
146        int *ti = Calloc(bnz, int), j, nz;
147        double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),
148            *tmp = Calloc(n, double);
149    
150        if (lo || (!unit))
151            error(_("Code written for unit upper triangular unit matrices"));
152        bp[0] = 0;
153        for (j = 0; j < n; j++) {
154            int i, i1 = ap[j + 1];
155            AZERO(tmp, n);
156            for (i = ap[j]; i < i1; i++) {
157                int ii = ai[i], k;
158                if (unit) tmp[ii] -= ax[i];
159                for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
160            }
161            for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
162            bp[j + 1] = bp[j] + nz;
163            if (bp[j + 1] > bnz) {
164                while (bp[j + 1] > bnz) bnz *= 2;
165                ti = Realloc(ti, bnz, int);
166                tx = Realloc(tx, bnz, double);
167            }
168            i1 = bp[j];
169            for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
170        }
171        nz = bp[n];
172        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
173        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
174        Free(tmp); Free(tx); Free(ti);
175        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
176        SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
177        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
178        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
179        UNPROTECT(1);
180        return ans;
181    }

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

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