SCM

SCM Repository

[matrix] Diff of /pkg/src/dsCMatrix.c
ViewVC logotype

Diff of /pkg/src/dsCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

pkg/src/sscMatrix.c revision 10, Mon Mar 22 20:20:05 2004 UTC pkg/src/dsCMatrix.c revision 683, Thu Mar 31 17:10:16 2005 UTC
# Line 1  Line 1 
1  #include "sscMatrix.h"  #include "dsCMatrix.h"
2    
3  SEXP sscMatrix_validate(SEXP obj)  SEXP dsCMatrix_validate(SEXP obj)
4  {  {
5      SEXP uplo = GET_SLOT(obj, Matrix_uploSym);      SEXP val = check_scalar_string(GET_SLOT(obj, Matrix_uploSym),
6                                       "LU", "uplo");
7      int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));      int *Dim = INTEGER(GET_SLOT(obj, Matrix_DimSym));
     char *val;  
8    
9      if (length(uplo) != 1)      if (isString(val)) return val;
         return ScalarString(mkChar("uplo slot must have length 1"));  
     val = CHAR(STRING_ELT(uplo, 0));  
     if (strlen(val) != 1)  
         return ScalarString(mkChar("uplo[1] must have string length 1"));  
     if (toupper(*val) != 'U' && toupper(*val) != 'L')  
         return ScalarString(mkChar("uplo[1] must be \"U\" or \"L\""));  
10      if (Dim[0] != Dim[1])      if (Dim[0] != Dim[1])
11          return ScalarString(mkChar("Symmetric matrix must be square"));          return mkString(_("Symmetric matrix must be square"));
12      csc_check_column_sorting(obj);      csc_check_column_sorting(obj);
13      return ScalarLogical(1);      return ScalarLogical(1);
14  }  }
15    
16  SEXP sscMatrix_chol(SEXP x, SEXP pivot)  SEXP dsCMatrix_chol(SEXP x, SEXP pivot)
17  {  {
18      SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("sscChol")));      SEXP pSlot = GET_SLOT(x, Matrix_pSym), xorig = x;
19      taucs_ccs_matrix *tm =      int *Ai = INTEGER(GET_SLOT(x, Matrix_iSym)),
20          csc_taucs_ptr(x, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC);          *Ap = INTEGER(pSlot),
21      int nnz, piv = asLogical(pivot);          *Lp, *Parent, info,
22      int *dims;          lo = CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L',
23      void *L;          n = length(pSlot)-1,
24            nnz, piv = asLogical(pivot);
25        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dCholCMatrix")));
26        int *P, *Pinv;
27        double *Ax;
28    
29        /* FIXME: Check if there is a Cholesky factorization.  If yes,
30           check if the permutation status matches that of the call.  If
31           so, return it. */
32    
33        if (lo) {
34            x = PROTECT(ssc_transpose(x));
35            Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
36            Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
37        }
38        SET_SLOT(val, Matrix_uploSym, mkString("L"));
39        SET_SLOT(val, Matrix_diagSym, mkString("U"));
40        SET_SLOT(val, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
41        SET_SLOT(val, Matrix_ParentSym, allocVector(INTSXP, n));
42        Parent = INTEGER(GET_SLOT(val, Matrix_ParentSym));
43        SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, n + 1));
44        Lp = INTEGER(GET_SLOT(val, Matrix_pSym));
45        SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, n));
46        P = INTEGER(GET_SLOT(val, Matrix_permSym));
47      if (piv) {      if (piv) {
48          int *iperm, *perm;          SEXP trip = PROTECT(dsCMatrix_to_dgTMatrix(x));
49            SEXP Ti = GET_SLOT(trip, Matrix_iSym);
50    
51            /* determine the permutation with Metis */
52            Pinv = Calloc(n, int);
53            ssc_metis_order(n, Ap, Ai, P, Pinv);
54            /* create a symmetrized form of x */
55            nnz = length(Ti);
56            Ai = Calloc(nnz, int);
57            Ax = Calloc(nnz, double);
58            Ap = Calloc(n + 1, int);
59            triplet_to_col(n, n, nnz, INTEGER(Ti),
60                           INTEGER(GET_SLOT(trip, Matrix_jSym)),
61                           REAL(GET_SLOT(trip, Matrix_xSym)),
62                           Ap, Ai, Ax);
63            UNPROTECT(1);
64        } else {
65            int i;
66            for (i = 0; i < n; i++) P[i] = i;
67            Ax = REAL(GET_SLOT(x, Matrix_xSym));
68    
69          SET_SLOT(val, Matrix_permSym, allocVector(INTSXP, tm->n));      }
70          SET_SLOT(val, Matrix_ipermSym, allocVector(INTSXP, tm->n));      R_ldl_symbolic(n, Ap, Ai, Lp, Parent, (piv) ? P : (int *)NULL,
71          perm = INTEGER(GET_SLOT(val, Matrix_permSym));                     (piv) ? Pinv : (int *)NULL);
72          iperm = INTEGER(GET_SLOT(val, Matrix_ipermSym));      nnz = Lp[n];
         ssc_metis_order(tm->n, (tm->colptr)[tm->n], tm->colptr,  
                         tm->rowind, perm, iperm);  
         tm = taucs_dccs_permute_symmetrically(tm, perm, iperm);  
     }  
     if (!(L = taucs_ccs_factor_llt_mf(tm)))  
         error("Matrix is not positive definite");  
     if (piv) taucs_dccs_free(tm);  
     tm = taucs_supernodal_factor_to_ccs(L);  
     taucs_supernodal_factor_free(L);  
     nnz = tm->colptr[tm->n];  
     SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, tm->n + 1));  
     Memcpy(INTEGER(GET_SLOT(val, Matrix_pSym)), tm->colptr, tm->n + 1);  
73      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));      SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
     Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), tm->rowind, nnz);  
74      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));      SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
75      Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), tm->values.d, nnz);      SET_SLOT(val, Matrix_DSym, allocVector(REALSXP, n));
76      dims = INTEGER(GET_SLOT(val, Matrix_DimSym));      info = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent,
77      dims[0] = dims[1] = tm->n;                           INTEGER(GET_SLOT(val, Matrix_iSym)),
78      taucs_dccs_free(tm);                           REAL(GET_SLOT(val, Matrix_xSym)),
79      UNPROTECT(1);                           REAL(GET_SLOT(val, Matrix_DSym)),
80      return set_factorization(x, val, "Cholesky");                           (piv) ? P : (int *)NULL,
81                             (piv) ? Pinv : (int *)NULL);
82        if (info != n)
83            error(_("Leading minor of size %d (possibly after permutation) is indefinite"),
84                  info + 1);
85        if (piv) {
86            Free(Pinv); Free(Ax); Free(Ai); Free(Ap);
87        }
88        UNPROTECT(lo ? 2 : 1);
89        return set_factors(xorig, val, "Cholesky");
90  }  }
91    
92  SEXP sscMatrix_matrix_solve(SEXP a, SEXP b)  SEXP dsCMatrix_matrix_solve(SEXP a, SEXP b, SEXP classed)
93  {  {
94      SEXP Chol = get_factorization(a, "Cholesky"),      int cl = asLogical(classed);
95          val = PROTECT(duplicate(b));      SEXP Chol = get_factors(a, "Cholesky"), perm,
96            bdP = cl ? GET_SLOT(b, Matrix_DimSym) : getAttrib(b, R_DimSymbol),
97            val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix")));
98      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),      int *adims = INTEGER(GET_SLOT(a, Matrix_DimSym)),
99          *bdims = INTEGER(getAttrib(b, R_DimSymbol)),          *bdims = INTEGER(bdP), *Li, *Lp, j, piv;
100          j, n = adims[1];      int n = adims[1], ncol = bdims[1];
101      taucs_ccs_matrix* tm;      double *Lx, *D, *in = REAL(cl ? GET_SLOT(b, Matrix_xSym) : b),
102      double *in = REAL(b), *out = REAL(val) ;          *out = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, n * ncol)),
103            *tmp = (double *) NULL;
104      if (!(isReal(b) && isMatrix(b)))  
105          error("Argument b must be a numeric matrix");      if (!cl && !(isReal(b) && isMatrix(b)))
106      if (*adims != *bdims || bdims[1] < 1 || *adims < 1)          error(_("Argument b must be a numeric matrix"));
107          error("Dimensions of system to be solved are inconsistent");      if (*adims != *bdims || ncol < 1 || *adims < 1)
108      if (Chol == R_NilValue) Chol = sscMatrix_chol(a, ScalarLogical(1.));          error(_("Dimensions of system to be solved are inconsistent"));
109      tm = csc_taucs_ptr(Chol, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_TRIANGULAR);      if (Chol == R_NilValue) Chol = dsCMatrix_chol(a, ScalarLogical(1));
110      if (!length(GET_SLOT(Chol, Matrix_permSym))) {      perm = GET_SLOT(Chol, Matrix_permSym);
111          for (j = 0; j < bdims[1]; j++, in += n, out += n) {      piv = length(perm);
112              int errcode = taucs_dccs_solve_llt(tm, out, in);      if (piv) tmp = Calloc(n, double);
113              if (errcode)      Li = INTEGER(GET_SLOT(Chol, Matrix_iSym));
114                  error("taucs_solve returned error code %d for column %d",      Lp = INTEGER(GET_SLOT(Chol, Matrix_pSym));
115                        errcode, j + 1);      Lx = REAL(GET_SLOT(Chol, Matrix_xSym));
116          }      D = REAL(GET_SLOT(Chol, Matrix_DSym));
117      } else {      for (j = 0; j < ncol; j++, in += n, out += n) {
118          int *iperm = INTEGER(GET_SLOT(Chol, Matrix_ipermSym));          if (piv) R_ldl_perm(n, out, in, INTEGER(perm));
119          double *tmpIn = (double *) R_alloc(n, sizeof(double)),          else Memcpy(out, in, n);
120              *tmpOut = (double *) R_alloc(n, sizeof(double));          R_ldl_lsolve(n, out, Lp, Li, Lx);
121            R_ldl_dsolve(n, out, D);
122          for (j = 0; j < bdims[1]; j++, in += n, out += n) {          R_ldl_ltsolve(n, out, Lp, Li, Lx);
123              int errcode, i;          if (piv) R_ldl_permt(n, out, Memcpy(tmp, out, n), INTEGER(perm));
                                 /* permute y */  
             for (i = 0; i < n; i++) tmpIn[iperm[i]] = in[i];  
                                 /* solve */  
             errcode = taucs_dccs_solve_llt(tm, tmpOut, tmpIn);  
             if (errcode)  
                 error("taucs_solve returned error code %d for column %d",  
                       errcode, j + 1);  
                                 /* inverse permute b */  
             for (i = 0; i < n; i++) out[i] = tmpOut[iperm[i]];  
         }  
124      }      }
125        if (piv) Free(tmp);
126      UNPROTECT(1);      UNPROTECT(1);
127      return val;      return val;
128  }  }
129    
130  SEXP sscMatrix_inverse_factor(SEXP A)  SEXP dsCMatrix_inverse_factor(SEXP A)
131  {  {
132      return mat_from_taucs(taucs_ccs_factor_xxt(      return R_NilValue;          /* FIXME: Write this function. */
         csc_taucs_ptr(A, TAUCS_DOUBLE | TAUCS_LOWER | TAUCS_SYMMETRIC)));  
133  }  }
134    
135  SEXP ssc_transpose(SEXP x)  SEXP ssc_transpose(SEXP x)
136  {  {
137      SEXP      SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))),
         ans = PROTECT(NEW_OBJECT(MAKE_CLASS("sscMatrix"))),  
138          islot = GET_SLOT(x, Matrix_iSym);          islot = GET_SLOT(x, Matrix_iSym);
139      int nnz = length(islot),      int nnz = length(islot), *adims,
         *adims = INTEGER(GET_SLOT(ans, Matrix_DimSym)),  
140          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));          *xdims = INTEGER(GET_SLOT(x, Matrix_DimSym));
141    
142        adims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
143      adims[0] = xdims[1]; adims[1] = xdims[0];      adims[0] = xdims[1]; adims[1] = xdims[0];
144      if (toupper(CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0]) == 'U')      if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'U')
145          SET_SLOT(ans, Matrix_uploSym, ScalarString(mkChar("L")));          SET_SLOT(ans, Matrix_uploSym, mkString("L"));
146      SET_SLOT(ans, Matrix_pSym, allocVector(INTSXP, xdims[0] + 1));      else
147      SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nnz));          SET_SLOT(ans, Matrix_uploSym, mkString("U"));
148      SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nnz));      csc_compTr(xdims[0], xdims[1], nnz,
149      csc_components_transpose(xdims[0], xdims[1], nnz,                 INTEGER(GET_SLOT(x, Matrix_pSym)), INTEGER(islot),
                              INTEGER(GET_SLOT(x, Matrix_pSym)),  
                              INTEGER(islot),  
150                               REAL(GET_SLOT(x, Matrix_xSym)),                               REAL(GET_SLOT(x, Matrix_xSym)),
151                               INTEGER(GET_SLOT(ans, Matrix_pSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, xdims[0] + 1)),
152                               INTEGER(GET_SLOT(ans, Matrix_iSym)),                 INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nnz)),
153                               REAL(GET_SLOT(ans, Matrix_xSym)));                 REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nnz)));
154      UNPROTECT(1);      UNPROTECT(1);
155      return ans;      return ans;
156  }  }
157    
158    SEXP dsCMatrix_to_dgTMatrix(SEXP x)
159    {
160        SEXP
161            ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
162            islot = GET_SLOT(x, Matrix_iSym),
163            pslot = GET_SLOT(x, Matrix_pSym);
164        int *ai, *aj, *iv = INTEGER(islot),
165            j, jj, nnz = length(islot), nout,
166            n = length(pslot) - 1,
167            *p = INTEGER(pslot), pos;
168        double *ax, *xv = REAL(GET_SLOT(x, Matrix_xSym));
169    
170        /* increment output count by number of off-diagonals */
171        nout = nnz;
172        for (j = 0; j < n; j++) {
173            int p2 = p[j+1];
174            for (jj = p[j]; jj < p2; jj++) {
175                if (iv[jj] != j) nout++;
176            }
177        }
178        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(x, Matrix_DimSym)));
179        SET_SLOT(ans, Matrix_iSym, allocVector(INTSXP, nout));
180        ai = INTEGER(GET_SLOT(ans, Matrix_iSym));
181        SET_SLOT(ans, Matrix_jSym, allocVector(INTSXP, nout));
182        aj = INTEGER(GET_SLOT(ans, Matrix_jSym));
183        SET_SLOT(ans, Matrix_xSym, allocVector(REALSXP, nout));
184        ax = REAL(GET_SLOT(ans, Matrix_xSym));
185        pos = 0;
186        for (j = 0; j < n; j++) {
187            int p2 = p[j+1];
188            for (jj = p[j]; jj < p2; jj++) {
189                int ii = iv[jj];
190                double xx = xv[jj];
191    
192                ai[pos] = ii; aj[pos] = j; ax[pos] = xx; pos++;
193                if (ii != j) {
194                    aj[pos] = ii; ai[pos] = j; ax[pos] = xx; pos++;
195                }
196            }
197        }
198        UNPROTECT(1);
199        return ans;
200    }
201    
202    SEXP dsCMatrix_ldl_symbolic(SEXP x, SEXP doPerm)
203    {
204        SEXP Ax, Dims = GET_SLOT(x, Matrix_DimSym),
205            ans = PROTECT(allocVector(VECSXP, 3)), tsc;
206        int i, n = INTEGER(Dims)[0], nz, nza,
207            *Ap, *Ai, *Lp, *Li, *Parent,
208            doperm = asLogical(doPerm),
209            *P = (int *) NULL, *Pinv = (int *) NULL;
210    
211    
212        if (CHAR(asChar(GET_SLOT(x, Matrix_uploSym)))[0] == 'L') {
213            x = PROTECT(ssc_transpose(x));
214        } else {
215            x = PROTECT(duplicate(x));
216        }
217        Ax = GET_SLOT(x, Matrix_xSym);
218        nza = length(Ax);
219        Ap = INTEGER(GET_SLOT(x, Matrix_pSym));
220        Ai = INTEGER(GET_SLOT(x, Matrix_iSym));
221        if (doperm) {
222            int *perm, *iperm = Calloc(n, int);
223    
224            SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, n));
225            perm = INTEGER(VECTOR_ELT(ans, 2));
226            ssc_metis_order(n, Ap, Ai, perm, iperm);
227            ssc_symbolic_permute(n, 1, iperm, Ap, Ai);
228            Free(iperm);
229        }
230        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
231        Parent = INTEGER(VECTOR_ELT(ans, 0));
232        SET_VECTOR_ELT(ans, 1, NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
233        tsc = VECTOR_ELT(ans, 1);
234        SET_SLOT(tsc, Matrix_uploSym, mkString("L"));
235        SET_SLOT(tsc, Matrix_diagSym, mkString("U"));
236        SET_SLOT(tsc, Matrix_DimSym, Dims);
237        SET_SLOT(tsc, Matrix_pSym, allocVector(INTSXP, n + 1));
238        Lp = INTEGER(GET_SLOT(tsc, Matrix_pSym));
239        R_ldl_symbolic(n, Ap, Ai, Lp, Parent, P, Pinv);
240        nz = Lp[n];
241        SET_SLOT(tsc, Matrix_iSym, allocVector(INTSXP, nz));
242        Li = INTEGER(GET_SLOT(tsc, Matrix_iSym));
243        SET_SLOT(tsc, Matrix_xSym, allocVector(REALSXP, nz));
244        for (i = 0; i < nza; i++) REAL(Ax)[i] = 0.00001;
245        for (i = 0; i < n; i++) REAL(Ax)[Ap[i+1]-1] = 10000.;
246        i = R_ldl_numeric(n, Ap, Ai, REAL(Ax), Lp, Parent, Li,
247                        REAL(GET_SLOT(tsc, Matrix_xSym)),
248                        (double *) R_alloc(n, sizeof(double)), /* D */
249                        P, Pinv);
250        UNPROTECT(2);
251        return ans;
252    }
253    
254    SEXP dsCMatrix_metis_perm(SEXP x)
255    {
256        SEXP pSlot = GET_SLOT(x, Matrix_pSym),
257            ans = PROTECT(allocVector(VECSXP, 2));
258        int n = length(pSlot) - 1;
259    
260        SET_VECTOR_ELT(ans, 0, allocVector(INTSXP, n));
261        SET_VECTOR_ELT(ans, 1, allocVector(INTSXP, n));
262        ssc_metis_order(n,
263                        INTEGER(pSlot),
264                        INTEGER(GET_SLOT(x, Matrix_iSym)),
265                        INTEGER(VECTOR_ELT(ans, 0)),
266                        INTEGER(VECTOR_ELT(ans, 1)));
267        UNPROTECT(1);
268        return ans;
269    }
270    

Legend:
Removed from v.10  
changed lines
  Added in v.683

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