SCM

SCM Repository

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

Diff of /pkg/src/ssclme.c

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

revision 125, Sun Apr 25 11:59:37 2004 UTC revision 157, Sun May 9 14:26:36 2004 UTC
# Line 1  Line 1 
1  #include "ssclme.h"  #include "ssclme.h"
2    
 /**  
  * Check for a nested series of grouping factors in the sparse,  
  *  symmetric representation of the pairwise cross-tabulations.  
  *  
  * @param n size of pairwise cross-tabulation matrix  
  * @param nf number of groups of columns in pairwise cross-tabulation  
  * @param upper non-zero if the upper triangle is stored  
  * @param Ap array of pointers to columns  
  * @param Ai row indices  
  * @param Gp array of pointers to groups  
  *  
  * @return 0 for non-nested groups, 1 for nested groups  
  */  
 static  
 int ctab_isNested(int n, int nf, int upper,  
                   const int Ap[], const int Ai[], const int Gp[])  
 {  
     if (nf > 1) {  /* single factor always nested */  
         int  i;  
         if (upper) {  
             int *nnz = (int *) R_alloc(n, sizeof(int)), nz = Ap[n];  
                                 /* count number of nonzeros in each row */  
             for (i = 0; i < n; i++) nnz[i] = 0;  
             for (i = 0; i < nz; i++) nnz[Ai[i]]++;  
             for (i = 0; i < nf; i++) {  
                 int j, p2 = Gp[i+1], target = nf - i;  
                 for (j = Gp[i]; j < p2; j++) {  
                     if (nnz[j] != target) return 0;  
                 }  
             }  
         } else {                /* lower triangle - the easy case */  
             for (i = 0; i < nf; i++) {  
                 int j, p2 = Gp[i+1], target = nf - i;  
                 for (j = Gp[i]; j < p2; j++) {  
                     if ((Ap[j+1] - Ap[j]) != target)  
                         return 0;  
                 }  
             }  
         }  
     }  
     return 1;  
 }  
   
 /**  
  * Determine if a fill-reducing permutation is needed for the pairwise  
  * cross-tabulation matrix.  If so, determine such a permutation  
  * (using Metis) then separate the groups.  
  *  
  * @param ctab pointer to a pairwise cross-tabulation object  
  *  
  * @return pointer to an integer R vector.  
  */  
   
 SEXP ctab_permute(SEXP ctab)  
 {  
     SEXP val, GpSl = GET_SLOT(ctab, Matrix_GpSym);  
     int *Ai = INTEGER(GET_SLOT(ctab, Matrix_iSym)),  
         *Ap = INTEGER(GET_SLOT(ctab, Matrix_pSym)),  
         *Gp = INTEGER(GpSl),  
         *perm,  
         *work,  
         i,  
         j,  
         n = INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],  
         nf = length(GpSl) - 1,  
         pos;  
   
     if (ctab_isNested(n, nf, 1, Ap, Ai, Gp))  
         return allocVector(INTSXP, 0);  
     val =  allocVector(INTSXP, n);  
     perm = INTEGER(val);  
     work = (int *) R_alloc(n, sizeof(int));  
     ssc_metis_order(n, Ap, Ai, work, perm);     /* perm gets inverse perm */  
     /* work now contains desired permutation but with groups scrambled */  
   
     /* copy work into perm preserving the order of the groups */  
     pos = 0;            /* position in new permutation */  
     for (i = 0; i < nf; i++) {  
         for (j = 0; j < n; j++) {  
             int jj = work[j];  
             if (Gp[i] <= jj && jj < Gp[i+1]) {  
                 perm[pos] = jj;  
                 pos++;  
             }  
         }  
     }  
     return val;  
 }  
   
3  static  static
4  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)  void ssclme_copy_ctab(int nf, const int nc[], SEXP ctab, SEXP ssc)
5  {  {
# Line 230  Line 141 
141      ssc = VECTOR_ELT(val, 0);      ssc = VECTOR_ELT(val, 0);
142                                  /* Pairwise cross-tabulation */                                  /* Pairwise cross-tabulation */
143      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));      ctab = PROTECT(sscCrosstab(facs, ScalarLogical(1)));
144      SET_VECTOR_ELT(val, 1, ctab_permute(ctab));      SET_VECTOR_ELT(val, 1, sscCrosstab_groupedPerm(ctab));
145      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */      if (length(VECTOR_ELT(val, 1)) > 0) {/* Fill-reducing permutation */
146          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],          ssc_symbolic_permute(INTEGER(GET_SLOT(ctab, Matrix_DimSym))[1],
147                               1, INTEGER(VECTOR_ELT(val, 1)),                               1, INTEGER(VECTOR_ELT(val, 1)),
# Line 369  Line 280 
280          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),          *Ap = INTEGER(GET_SLOT(x, Matrix_pSym)),
281          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),          *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
282          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),          *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
283            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym)),
284          i, j, k,          i, j, k,
285          ione = 1,          ione = 1,
286          nf = length(mmats) - 1,          nf = length(mmats) - 1,
# Line 471  Line 383 
383      }      }
384      Free(Z);      Free(Z);
385      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
386        status[0] = status[1] = 0;
387      return R_NilValue;      return R_NilValue;
388  }  }
389    
# Line 742  Line 655 
655      }      }
656      if (length(LIisl)) error("logic error in ssclme_ldl_inverse");      if (length(LIisl)) error("logic error in ssclme_ldl_inverse");
657      else {                  /* LIi and LIx are too big and not used */      else {                  /* LIi and LIx are too big and not used */
658          int *counts = Calloc(nzc, int), info, maxod = -1;          int info, maxod;
659          int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));          int *Parent = INTEGER(GET_SLOT(x, Matrix_ParentSym));
660          int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));          int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym));
661          double one = 1.0, zero = 0.;          double one = 1.0, zero = 0.;
662                                  /* determine maximum # of off-diagonals */  
663          for (i = nzc - 1; i >= 0; i--) { /* in a column of L^{-1} */          maxod = -1;         /* determine maximum # of off-diagonals */
664              counts[i] = (Parent[i] < 0) ? 0 : 1 + counts[Parent[i]];          for (i = 0; i < nzc; i++) { /* in a column of L^{-1} */
665              if (counts[i] > maxod) maxod = counts[i];              int ci = LIp[i+1] - LIp[i];
666                if (ci > maxod) maxod = ci;
667          }          }
         Free(counts);  
668    
669          for (i = 0; i < nf; i++) {          for (i = 0; i < nf; i++) {
670              int j, jj, k, kk, nci = nc[i], nr, p, p2, pp,              int j, jj, k, kk, nci = nc[i], nr, p, p2, pp,
671                  m = maxod + nci,                  m = maxod + 1,
672                  *ind = Calloc(m, int);                  *ind = Calloc(m, int);
673              double              double
674                  *tmp = Calloc(m * nci, double),                  *tmp = Calloc(m * nci, double),
675                  *mpt = REAL(VECTOR_ELT(bVar, i));                  *mpt = REAL(VECTOR_ELT(bVar, i));
676    
677              for (j = Gp[i]; j < Gp[i+1]; j += nci) {              for (j = Gp[i]; j < Gp[i+1]; j += nci) {
                 memset((void *) tmp, 0, sizeof(double) * m * nci);  
   
678                  kk = 0;         /* ind holds indices of non-zeros */                  kk = 0;         /* ind holds indices of non-zeros */
679                  jj = j;         /* in this block of columns */                  jj = j;         /* in this block of columns */
680                  while (jj >= 0) {                  while (jj >= 0) {
# Line 776  Line 687 
687                  for (k = 0; k < nci; k++) {                  for (k = 0; k < nci; k++) {
688                      double *ccol = tmp + k * nr;                      double *ccol = tmp + k * nr;
689    
690                      ccol[k] = 1.;                      for (kk = 0; kk < nr; kk++) ccol[kk] = 0.;
691                      kk = k;                      ccol[k] = 1.; /* initialize from unit diagonal */
692                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {                      for (jj = j + k; jj >= 0; jj = Parent[jj]) {
693                          p2 = Lp[jj+1];                          p2 = Lp[jj+1];
694                          pp = kk;                          pp = k;
695                          for (p = Lp[jj]; p < p2; p++) {                          for (p = Lp[jj]; p < p2; p++) {
696                              pp = ldl_update_ind(Li[p], pp, ind);                              pp = ldl_update_ind(Li[p], pp, ind);
697                              ccol[pp] -= Lx[p] * ccol[kk];                              ccol[pp] -= Lx[p] * ccol[jj];
698                          }                          }
699                      }                      }
700    
701    /*                  for (j = k + 1; j < nr; j++) {  */
702    /*                      jj = ind[j]; p2 = Lp[jj+1]; */
703    /*                      pp = k; */
704    /*                      for (p = Lp[jj]; p < p2; p++) { */
705    /*                          pp = ldl_update_ind(Li[p], pp, ind); */
706    /*                          ccol[pp] -= Lx[p]*ccol[j]; */
707    /*                      } */
708    /*                  } */
709                  }                  }
710                                  /* scale rows */  
711                  for (kk = 0; kk < nr; kk++) {                  for (kk = 0; kk < nr; kk++) { /* scale rows */
712                      for (k = 0; k < nci; k++) {                      for (k = 0; k < nci; k++) {
713                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];                          tmp[k * nr + kk] *= DIsqrt[ind[kk]];
714                      }                      }
# Line 1263  Line 1183 
1183          alpha,          alpha,
1184          one = 1.;          one = 1.;
1185      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));      SEXP ans = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));  
1186    
1187      ssclme_factor(x);      ssclme_factor(x);
1188      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {      if (!R_FINITE(REAL(GET_SLOT(x, Matrix_devianceSym))[0])) {

Legend:
Removed from v.125  
changed lines
  Added in v.157

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