SCM

SCM Repository

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

Diff of /pkg/src/dgCMatrix.c

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

revision 2298, Fri Oct 17 13:45:12 2008 UTC revision 2299, Fri Oct 17 16:07:38 2008 UTC
# Line 272  Line 272 
272  SEXP dgCMatrix_SPQR(SEXP Ap, SEXP order)  SEXP dgCMatrix_SPQR(SEXP Ap, SEXP order)
273  {  {
274      SEXP ans = PROTECT(allocVector(VECSXP, 4));      SEXP ans = PROTECT(allocVector(VECSXP, 4));
275      CHM_SP A = AS_CHM_SP(Ap), Al, Q, R;      CHM_SP A = AS_CHM_SP(Ap), Q, R;
276      UF_long *Ali, *Alp, *E, rank;      int *E, rank;
     int annz, *Ai = (int*)(A->i), *App = (int*)(A->p);  
     double *Ax = (double*)(A->x), *Alx;  
   
     annz = cholmod_nnz(A, &c);  
   
     /* incorporate this into a utility */  
     if (!(Al = cholmod_l_allocate_sparse(A->nrow, A->ncol,  
                                          (size_t) annz, A->sorted,  
                                          A->packed, A->stype,  
                                          A->xtype, &cl)))  
         error(_("cholmod_l_allocate_sparse error, status = "));  
     Ali = (UF_long*)(Al->i); Alp = (UF_long*)(Al->p); Alx = (double*)(Al->x);  
     for (int j = 0; j <= A->ncol; j++) Alp[j] = (UF_long)(App[j]);  
     for (int p = 0; p < annz; p++) {  
         Ali[p] = (UF_long)(Ai[p]);  
         Alx[p] = Ax[p];  
     }  
     /* end of utility */  
277    
278      if ((rank = SuiteSparseQR_C_QR(asInteger(order),      if ((rank = SuiteSparseQR_C_QR(asInteger(order),
279                                     SPQR_DEFAULT_TOL, 0,                                     SPQR_DEFAULT_TOL, 0,
280                                     A, &Q, &R, &E, &cl)) == -1)                                     A, &Q, &R, &E, &c)) == -1)
281          error(_("SuiteSparseQR_C_QR returned an error code"));          error(_("SuiteSparseQR_C_QR returned an error code"));
282      SET_VECTOR_ELT(ans, 0,      SET_VECTOR_ELT(ans, 0,
283                     chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue));                     chm_sparse_to_SEXP(Q, 0, 0, 0, "", R_NilValue));
284      SET_VECTOR_ELT(ans, 1,      SET_VECTOR_ELT(ans, 1,
285                     chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue));                     chm_sparse_to_SEXP(R, 0, 0, 0, "", R_NilValue));
286      cholmod_l_free_sparse(&Al, &cl);      cholmod_l_free_sparse(&R, &c);
287      cholmod_l_free_sparse(&R, &cl);      cholmod_l_free_sparse(&Q, &c);
     cholmod_l_free_sparse(&Q, &cl);  
288      if (E) {      if (E) {
         int *Er;  
289          SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, A->ncol));          SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, A->ncol));
290          Er = INTEGER(VECTOR_ELT(ans, 2));          Memcpy(INTEGER(VECTOR_ELT(ans, 2)), E, A->ncol);
         for (int i = 0; i < A->ncol; i++) Er[i] = (int) E[i];  
291          Free(E);          Free(E);
292      } else SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 0));      } else SET_VECTOR_ELT(ans, 2, allocVector(INTSXP, 0));
293      SET_VECTOR_ELT(ans, 3, ScalarInteger((int)rank));      SET_VECTOR_ELT(ans, 3, ScalarInteger(rank));
294      UNPROTECT(1);      UNPROTECT(1);
295      return ans;      return ans;
296  }  }
# Line 433  Line 412 
412          error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix"));          error(_("dgCMatrix_cholsol requires a 'short, wide' rectangular matrix"));
413      if (cy->nrow != cx->ncol)      if (cy->nrow != cx->ncol)
414          error(_("Dimensions of system to be solved are inconsistent"));          error(_("Dimensions of system to be solved are inconsistent"));
415      rhs = cholmod_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c);      rhs = cholmod_l_allocate_dense(cx->nrow, 1, cx->nrow, CHOLMOD_REAL, &c);
416      if (!(cholmod_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c)))      if (!(cholmod_l_sdmult(cx, 0 /* trans */, one, zero, cy, rhs, &c)))
417          error(_("cholmod_sdmult error"));          error(_("cholmod_l_sdmult error"));
418      L = cholmod_analyze(cx, &c);      L = cholmod_l_analyze(cx, &c);
419      if (!cholmod_factorize(cx, L, &c))      if (!cholmod_l_factorize(cx, L, &c))
420          error(_("cholmod_factorize failed: status %d, minor %d from ncol %d"),          error(_("cholmod_l_factorize failed: status %d, minor %d from ncol %d"),
421                c.status, L->minor, L->n);                c.status, L->minor, L->n);
422  /* FIXME: Do this in stages so an "effects" vector can be calculated */  /* FIXME: Do this in stages so an "effects" vector can be calculated */
423      if (!(cAns = cholmod_solve(CHOLMOD_A, L, rhs, &c)))      if (!(cAns = cholmod_l_solve(CHOLMOD_A, L, rhs, &c)))
424          error(_("cholmod_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"),          error(_("cholmod_l_solve (CHOLMOD_A) failed: status %d, minor %d from ncol %d"),
425                c.status, L->minor, L->n);                c.status, L->minor, L->n);
426      SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0));      SET_VECTOR_ELT(ans, 0, chm_factor_to_SEXP(L, 0));
427      SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow));      SET_VECTOR_ELT(ans, 1, allocVector(REALSXP, cx->nrow));
# Line 451  Line 430 
430      SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow));      SET_VECTOR_ELT(ans, 2, allocVector(REALSXP, cx->nrow));
431      Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(rhs->x), cx->nrow);      Memcpy(REAL(VECTOR_ELT(ans, 1)), (double*)(rhs->x), cx->nrow);
432    
433      cholmod_free_factor(&L, &c);      cholmod_l_free_factor(&L, &c);
434      cholmod_free_dense(&rhs, &c);      cholmod_l_free_dense(&rhs, &c);
435      cholmod_free_dense(&cAns, &c);      cholmod_l_free_dense(&cAns, &c);
436      UNPROTECT(1);      UNPROTECT(1);
437      return ans;      return ans;
438  }  }

Legend:
Removed from v.2298  
changed lines
  Added in v.2299

R-Forge@R-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business University of Wisconsin - Madison Powered By FusionForge