SCM

SCM Repository

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

Diff of /pkg/src/bCrosstab.c

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

revision 409, Sun Dec 19 02:19:58 2004 UTC revision 410, Tue Dec 28 14:51:42 2004 UTC
# Line 1  Line 1 
1  #include "bCrosstab.h"  #include "bCrosstab.h"
2    /* TODO:
3     * Watch distinction between ZtZ and ZZpO slots.
4     * Add code to check for a dense nnz before calculating a fill-reducing permutation.
5     */
6    
7    /**
8     * Allocate a 3-dimensional array
9     *
10     * @param TYP The R Type code (e.g. INTSXP)
11     * @param nr number of rows
12     * @param nc number of columns
13     * @param nf number of faces
14     *
15     * @return A 3-dimensional array of the indicated dimensions and type
16     */
17    static
18    SEXP alloc3Darray(int TYP, int nr, int nc, int nf)
19    {
20        SEXP val, dd = PROTECT(allocVector(INTSXP, 3));
21    
22        INTEGER(dd)[0] = nr; INTEGER(dd)[1] = nc; INTEGER(dd)[2] = nf;
23        val = allocArray(TYP, dd);
24        UNPROTECT(1);
25        return val;
26    }
27    
28  /**  /**
29   * Calculate the zero-based index in a row-wise packed lower   * Calculate the zero-based index in a row-wise packed lower
# Line 51  Line 76 
76  }  }
77    
78  /**  /**
79   * Create a named list of length n   * Create a named vector of type TYP
80   *   *
81   * @param n length of list to return   * @param TYP a vector SEXP type (e.g. REALSXP)
82   * @param names names of list elements   * @param names names of list elements with null string appended
83   *   *
84   * @return pointer to a named list of length n   * @return pointer to a named vector of type TYP
85   */   */
86  static SEXP  static SEXP
87  make_named_list(int n, char **names)  make_named(int TYP, char **names)
88  {  {
89      SEXP ans = PROTECT(allocVector(VECSXP, n)),      SEXP ans, nms;
90          nms = PROTECT(allocVector(STRSXP, n));      int i, n;
     int i;  
91    
92        for (n = 0; strlen(names[n]) > 0; n++) {}
93        ans = PROTECT(allocVector(TYP, n));
94        nms = PROTECT(allocVector(STRSXP, n));
95      for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i]));      for (i = 0; i < n; i++) SET_STRING_ELT(nms, i, mkChar(names[i]));
96      setAttrib(ans, R_NamesSymbol, nms);      setAttrib(ans, R_NamesSymbol, nms);
97      UNPROTECT(2);      UNPROTECT(2);
# Line 96  Line 123 
123  /**  /**
124   * Replace the structure of C by the structure of CA   * Replace the structure of C by the structure of CA
125   *   *
126   * @param A a unit lower tscMatrix   * @param A a unit lower triangular cscBlocked object
127   * @param C cscMatrix to be updated   * @param C a cscBlocked object to be updated
128   */   */
129  static void  static void
130  symbolic_right_unit_mm(SEXP A, SEXP C)  symbolic_right_unit_mm(SEXP A, SEXP C)
# Line 148  Line 175 
175          }          }
176      }      }
177      if (ncp[anc] > cp[anc]) {      if (ncp[anc] > cp[anc]) {
178          int *nci, nnz = ncp[anc], pos = 0;          int *dims, *nci, nnz = ncp[anc], pos = 0;
179          double *ncx;          double *ncx;
180    
181          SET_SLOT(C, Matrix_iSym, allocVector(INTSXP,  nnz));          SET_SLOT(C, Matrix_iSym, allocVector(INTSXP,  nnz));
182          nci = INTEGER(GET_SLOT(C, Matrix_iSym));          nci = INTEGER(GET_SLOT(C, Matrix_iSym));
183          SET_SLOT(C, Matrix_xSym, allocVector(REALSXP,  nnz));          dims = INTEGER(getAttrib(GET_SLOT(C, Matrix_xSym), R_DimSymbol));
184            SET_SLOT(C, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
185          ncx = REAL(GET_SLOT(C, Matrix_xSym));          ncx = REAL(GET_SLOT(C, Matrix_xSym));
186          for (i = 0; i < nnz; i++) ncx[i] = 1.;          for (i = 0; i < nnz; i++) ncx[i] = 1.;
187                                  /* As Diana Krall said, "Just do it again." */                                  /* As Diana Krall said, "Just do it again." */
# Line 217  Line 245 
245      if (!extra) return;      if (!extra) return;
246      {      {
247          int pos, nnz = tp[tnc];          int pos, nnz = tp[tnc];
248          int ntot = nnz + extra;          int ntot = nnz + extra, tnr;
249          int *Ai = Calloc(ntot, int),          int *Ai = Calloc(ntot, int),
250              *Ti = Calloc(ntot, int),              *Ti = Calloc(ntot, int),
251              *Tj = Calloc(ntot, int),              *Tj = Calloc(ntot, int),
252              *Dims = INTEGER(GET_SLOT(tb, Matrix_DimSym));              *dims;
253          double *Ax;          double *Ax;
254    
255          Memcpy(Ti, ti, nnz);    /* make a copy of the row indices */          Memcpy(Ti, ti, nnz);    /* make a copy of the row indices */
# Line 243  Line 271 
271                  }                  }
272              }              }
273          }          }
274          triplet_to_col(Dims[0], Dims[1], ntot, Ti, Tj, (double *) NULL,          /* Pass nlev instead of doing this.  The dimensions are nlev[i], nlev[k] */
275            /* Determine maximum row index in T */
276            tnr = -1; for (jj = 0; jj < ntot; jj++) if (Ti[jj] > tnr) tnr = Ti[jj];
277            tnr++;                  /* increment by 1 to get number of rows */
278            triplet_to_col(tnr, tnc, ntot, Ti, Tj, (double *) NULL,
279                         tp, Ai, (double *) NULL);                         tp, Ai, (double *) NULL);
280          nnz = tp[tnc];          nnz = tp[tnc];
281          SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));          SET_SLOT(tb, Matrix_iSym, allocVector(INTSXP, nnz));
282          Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);          Memcpy(INTEGER(GET_SLOT(tb, Matrix_iSym)), Ai, nnz);
283          SET_SLOT(tb, Matrix_xSym, allocVector(REALSXP, nnz));          dims = INTEGER(getAttrib(GET_SLOT(tb, Matrix_xSym), R_DimSymbol));
284            SET_SLOT(tb, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
285          Ax = REAL(GET_SLOT(tb, Matrix_xSym));          Ax = REAL(GET_SLOT(tb, Matrix_xSym));
286          for (j = 0; j < nnz; j++) Ax[j] = 1.;          for (j = 0; j < nnz; j++) Ax[j] = 1.;
287          Free(Ai); Free(Ti); Free(Tj);          Free(Ai); Free(Ti); Free(Tj);
# Line 256  Line 289 
289      }      }
290  }  }
291    
292    static int*
293    expand_column_pointers(int ncol, int nnz, const int mp[], int mj[])
294    {
295        int j;
296        for (j = 0; j < ncol; j++) {
297            int j2 = mp[j+1], jj;
298            for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;
299        }
300        return mj;
301    }
302    
303  /**  /**
304   * Permute the levels of one of the grouping factors in a bCrosstab object   * Permute the levels of one of the grouping factors in a bCrosstab object
305   *   *
# Line 267  Line 311 
311   * @param pperm inverse of the permutation   * @param pperm inverse of the permutation
312   */   */
313  static void  static void
314  bCrosstab_permute(SEXP ctab, int nf, int jj,  bCrosstab_permute(SEXP ctab, int nf, int jj, const int nlev[],
315                    const int perm[], const int iperm[])                    const int perm[], const int iperm[])
316  {  {
317      SEXP trip, ipt, jpt;      int j;
     int j, nnz;  
   
318      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
319          int ind = (j < jj ? Lind(jj, j) : Lind(j, jj));          int ind = (j < jj ? Lind(jj, j) : Lind(j, jj)),
320          SEXP csc = VECTOR_ELT(ctab, ind);              ncol = (j < jj ? nlev[j] : nlev[jj]),
321          int *Dims = INTEGER(GET_SLOT(csc, Matrix_DimSym));              nrow = (j < jj ? nlev[jj] : nlev[j]);
322            SEXP cscb = VECTOR_ELT(ctab, ind),
323          trip = csc_to_triplet(csc);              cscbi = GET_SLOT(cscb, Matrix_iSym);
324          ipt = GET_SLOT(trip, Matrix_iSym);          int *cp = INTEGER(GET_SLOT(cscb, Matrix_pSym)),
325          nnz = length(ipt);              nnz = length(cscbi);
326          jpt = GET_SLOT(trip, Matrix_jSym);          double *cx = REAL(GET_SLOT(cscb, Matrix_xSym));
327          if (j <= jj) ind_permute(INTEGER(ipt), nnz, iperm);          int *mj = expand_column_pointers(ncol, nnz, cp, Calloc(nnz, int));
328          if (j >= jj) ind_permute(INTEGER(jpt), nnz, iperm);          int *mi = Memcpy(Calloc(nnz, int), INTEGER(cscbi), nnz);
329          if (j == jj)          double *mx = Memcpy(Calloc(nnz, double), cx, nnz);
330              make_upper_triangular(INTEGER(ipt), INTEGER(jpt), nnz);  
331          triplet_to_col(Dims[0], Dims[1], nnz, INTEGER(ipt),          if (j <= jj) ind_permute(mi, nnz, iperm);
332                         INTEGER(jpt),          if (j >= jj) ind_permute(mj, nnz, iperm);
333                         REAL(GET_SLOT(trip, Matrix_xSym)),          if (j == jj) make_upper_triangular(mi, mj, nnz);
334                         INTEGER(GET_SLOT(csc, Matrix_pSym)),          triplet_to_col(nrow, ncol, nnz, mi, mj, mx, cp, INTEGER(cscbi), cx);
335                         INTEGER(GET_SLOT(csc, Matrix_iSym)),          Free(mi); Free(mj); Free(mx);
                        REAL(GET_SLOT(csc, Matrix_xSym)));  
336      }      }
337  }  }
338    
# Line 326  Line 367 
367  }  }
368    
369  /**  /**
370   * Derive the structure of elements of an lmeRep object from a blocked   * Create and populate slots in an lmer object from the blocked crosstabulation.
  * crosstabulation.  
371   *   *
372   * @param bCtab Pointer to a list containing a blocked crosstabulation   * @param val Pointer to an lmer object
373   * object and the factors used to generate it.   */
374   *  void
375   * @return A list of the factor list, the blocked/blocked sparse unit  lmer_populate(SEXP val)
376   * lower triangular matrix L, the list of diagonal elements of  {
377   * L-inverse, the permutations and the Parent arrays.      SEXP D, L, Linv, Parent, cscB = MAKE_CLASS("cscBlocked"), ZZpO,
378   */          flist = GET_SLOT(val, Matrix_flistSym), perm, Omega,
379  SEXP          ZtZ = GET_SLOT(val, Matrix_ZtZSym);
380  bCrosstab_convert(SEXP bCtab)      int j, nf = length(flist);
381  {      int *nc = INTEGER(GET_SLOT(val, Matrix_ncSym)), *Gp,
382      char *anms[] = {"flist", "L", "Linv", "perm", "Parent"};          *nlev = Calloc(nf, int), npairs = (nf * (nf + 1))/2;
383      SEXP flist = VECTOR_ELT(bCtab, 0),      char *statnms[] = {"factored", "inverted", ""},
384          ctab = PROTECT(duplicate(VECTOR_ELT(bCtab, 1))),          *devnms[] = {"ML", "REML", ""};
385          ans = PROTECT(make_named_list(sizeof(anms)/sizeof(char *), anms));  
386      SEXP fcp, perm, L, Linv, Parent;      /* Allocate fixed-sized slots */
387      int ctbl = length(ctab), j, nf = length(flist);      SET_SLOT(val, Matrix_statusSym, make_named(LGLSXP, statnms ));
388        SET_SLOT(val, Matrix_devianceSym, make_named(REALSXP, devnms));
389      if (ctbl != (nf*(nf + 1))/2)      SET_SLOT(val, Matrix_devCompSym, allocVector(REALSXP, 4));
390          error("length of ctab = %d is not permisable", ctbl);      /* Copy ZtZ to ZZpO */
391      SET_VECTOR_ELT(ans, 0, duplicate(flist));      SET_SLOT(val, Matrix_ZZpOSym, duplicate(ZtZ));
392      SET_VECTOR_ELT(ans, 1, allocVector(VECSXP, ctbl));      ZZpO = GET_SLOT(val, Matrix_ZZpOSym);
393      SET_VECTOR_ELT(ans, 2, allocVector(VECSXP, nf));      /* Allocate slots that are lists of length nf */
394      SET_VECTOR_ELT(ans, 3, allocVector(VECSXP, nf));      SET_SLOT(val, Matrix_DSym, allocVector(VECSXP, nf));
395      SET_VECTOR_ELT(ans, 4, allocVector(VECSXP, nf));      D = GET_SLOT(val, Matrix_DSym);
396      fcp = VECTOR_ELT(ans, 0);      setAttrib(D, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
397      L = VECTOR_ELT(ans, 1);      SET_SLOT(val, Matrix_LinvSym, allocVector(VECSXP, nf));
398      setAttrib(L, R_NamesSymbol, duplicate(getAttrib(ctab, R_NamesSymbol)));      Linv = GET_SLOT(val, Matrix_LinvSym);
     Linv = VECTOR_ELT(ans, 2);  
399      setAttrib(Linv, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(Linv, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
400      perm = VECTOR_ELT(ans, 3);      SET_SLOT(val, Matrix_permSym, allocVector(VECSXP, nf));
401        perm = GET_SLOT(val, Matrix_permSym);
402      setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(perm, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
403      Parent = VECTOR_ELT(ans, 4);      SET_SLOT(val, Matrix_ParentSym, allocVector(VECSXP, nf));
404        Parent = GET_SLOT(val, Matrix_ParentSym);
405      setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));      setAttrib(Parent, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
406        SET_SLOT(val, Matrix_OmegaSym, allocVector(VECSXP, nf));
407        Omega = GET_SLOT(val, Matrix_OmegaSym);
408        setAttrib(Omega, R_NamesSymbol, duplicate(getAttrib(flist, R_NamesSymbol)));
409    
410        /* Allocate peculiar length slots */
411        SET_SLOT(val, Matrix_LSym, allocVector(VECSXP, npairs));
412        L = GET_SLOT(val, Matrix_LSym);
413        SET_SLOT(val, Matrix_GpSym, allocVector(INTSXP, nf + 1));
414        Gp = INTEGER(GET_SLOT(val, Matrix_GpSym));
415        Gp[0] = 0;
416        for (j = 0; j < nf; j++) {
417            nlev[j] = length(getAttrib(VECTOR_ELT(flist, j), R_LevelsSymbol));
418            Gp[j + 1] = Gp[j] + nc[j] * nlev[j];
419            SET_VECTOR_ELT(D, j, alloc3Darray(REALSXP, nc[j], nc[j], nlev[j]));
420            SET_VECTOR_ELT(Omega, j, allocMatrix(REALSXP, nc[j], nc[j]));
421        }
422        SET_SLOT(val, Matrix_XtXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));
423        SET_SLOT(val, Matrix_RXXSym, allocMatrix(REALSXP, nc[nf], nc[nf]));
424        Memcpy(REAL(GET_SLOT(val, Matrix_RXXSym)),
425               memset(REAL(GET_SLOT(val, Matrix_XtXSym)), 0,
426                      sizeof(double) * nc[nf] * nc[nf]),
427               nc[nf] * nc[nf]);
428        SET_SLOT(val, Matrix_ZtXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));
429        SET_SLOT(val, Matrix_RZXSym, allocMatrix(REALSXP, Gp[nf], nc[nf]));
430      for (j = 0; j < nf; j++) {      for (j = 0; j < nf; j++) {
431          int dind = Lind(j, j), i, k;          int dind = Lind(j, j), i, k;
432          SEXP ctd = VECTOR_ELT(ctab, dind); /* diagonal in crosstab */          SEXP ctd = VECTOR_ELT(ZZpO, dind); /* diagonal in crosstab */
433          SEXP Dimslot = GET_SLOT(ctd, Matrix_DimSym),          SEXP Linvj, Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),
             Linvj, Ljj, cpp = GET_SLOT(ctd, Matrix_pSym),  
434              cip = GET_SLOT(ctd, Matrix_iSym);              cip = GET_SLOT(ctd, Matrix_iSym);
435          int *Lp, *Linvp, *Perm, *cp = INTEGER(cpp),          int *Lp, *Linvp, *Perm, *cp = INTEGER(cpp),
436              *ci = INTEGER(cip), *parent,              *ci = INTEGER(cip), *dims, *parent,
437              ncj = length(cpp) - 1,              ncj = length(cpp) - 1,
438              nnz = length(cip);              nnz = length(cip);
439          double *dtmp;          double *dtmp;
# Line 378  Line 442 
442          parent = INTEGER(VECTOR_ELT(Parent, j));          parent = INTEGER(VECTOR_ELT(Parent, j));
443          SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));          SET_VECTOR_ELT(perm, j, allocVector(INTSXP, ncj));
444          Perm = INTEGER(VECTOR_ELT(perm, j));          Perm = INTEGER(VECTOR_ELT(perm, j));
445          SET_VECTOR_ELT(L, dind, NEW_OBJECT(MAKE_CLASS("cscMatrix")));          SET_VECTOR_ELT(L, dind, NEW_OBJECT(cscB));
446          Ljj = VECTOR_ELT(L, dind);          Ljj = VECTOR_ELT(L, dind);
447          SET_VECTOR_ELT(Linv, j, NEW_OBJECT(MAKE_CLASS("cscMatrix")));          SET_VECTOR_ELT(Linv, j, NEW_OBJECT(cscB));
448          Linvj = VECTOR_ELT(Linv, j);          Linvj = VECTOR_ELT(Linv, j);
         SET_SLOT(Ljj, Matrix_DimSym, duplicate(Dimslot));  
         SET_SLOT(Linvj, Matrix_DimSym, duplicate(Dimslot));  
         SET_SLOT(Ljj, Matrix_factorization, allocVector(VECSXP, 0));  
         SET_SLOT(Linvj, Matrix_factorization, allocVector(VECSXP, 0));  
449          SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));          SET_SLOT(Ljj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
450          SET_SLOT(Linvj, Matrix_pSym, allocVector(INTSXP, ncj + 1));          SET_SLOT(Linvj, Matrix_pSym, allocVector(INTSXP, ncj + 1));
451          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));          Lp = INTEGER(GET_SLOT(Ljj, Matrix_pSym));
452          Linvp = INTEGER(GET_SLOT(Linvj, Matrix_pSym));          Linvp = INTEGER(GET_SLOT(Linvj, Matrix_pSym));
453            dims = INTEGER(getAttrib(GET_SLOT(ctd, Matrix_xSym), R_DimSymbol));
454          if (nnz > ncj) {        /* calculate fill-reducing permutation */          if (nnz > ncj) {        /* calculate fill-reducing permutation */
455                SEXP fac = VECTOR_ELT(flist, j);
456                SEXP fcp = PROTECT(duplicate(fac));
457              int *iPerm = Calloc(ncj, int);              int *iPerm = Calloc(ncj, int);
458    
459              ssc_metis_order(ncj, cp, ci, Perm, iPerm);              ssc_metis_order(ncj, cp, ci, Perm, iPerm);
460                                  /* apply to the crosstabulation */                                  /* apply to the crosstabulation and ZZpO */
461              bCrosstab_permute(ctab, nf, j, Perm, iPerm);              bCrosstab_permute(ZtZ, nf, j, nlev, Perm, iPerm);
462                bCrosstab_permute(ZZpO, nf, j, nlev, Perm, iPerm);
463                                  /* apply to the factor */                                  /* apply to the factor */
464              factor_levels_permute(VECTOR_ELT(fcp, j),              factor_levels_permute(fac, fcp, Perm, iPerm);
                                   VECTOR_ELT(flist, j), Perm, iPerm);  
465                                  /* symbolic analysis to get Parent */                                  /* symbolic analysis to get Parent */
466              R_ldl_symbolic(ncj, cp, ci, Lp, parent,              R_ldl_symbolic(ncj, cp, ci, Lp, parent,
467                           (int *) NULL, (int *) NULL);                           (int *) NULL, (int *) NULL);
# Line 409  Line 472 
472              for (i = 0; i < ncj; i++) dtmp[cp[i+1] - 1] = 1.;              for (i = 0; i < ncj; i++) dtmp[cp[i+1] - 1] = 1.;
473              nnz = Lp[ncj];              nnz = Lp[ncj];
474              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, nnz));
475              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, nnz));              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
476              i = R_ldl_numeric(ncj, cp, ci, dtmp, Lp, parent,              i = R_ldl_numeric(ncj, cp, ci, dtmp, Lp, parent,
477                                 INTEGER(GET_SLOT(Ljj, Matrix_iSym)),                                 INTEGER(GET_SLOT(Ljj, Matrix_iSym)),
478                                 REAL(GET_SLOT(Ljj, Matrix_xSym)),                                 REAL(GET_SLOT(Ljj, Matrix_xSym)),
479                                 (double *) R_alloc(ncj, sizeof(double)), /* D */                                 (double *) R_alloc(ncj, sizeof(double)), /* D */
480                                 (int *) NULL, (int *) NULL);                                 (int *) NULL, (int *) NULL);
481              if (i < ncj) error("identity matrix is not positive definite?");              if (i < ncj) error("identity matrix is not positive definite?");
482              Free(iPerm);              Free(iPerm); UNPROTECT(1);
483          } else {          } else {
484              for (i = 0; i <= ncj; i++) {              for (i = 0; i <= ncj; i++) {
485                  Lp[i] = 0;                  Lp[i] = 0;
# Line 424  Line 487 
487                  Perm[i] = i;                  Perm[i] = i;
488              }              }
489              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));              SET_SLOT(Ljj, Matrix_iSym, allocVector(INTSXP, 0));
490              SET_SLOT(Ljj, Matrix_xSym, allocVector(REALSXP, 0));              SET_SLOT(Ljj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], 0));
491          }          }
492                                  /* Derive the diagonal block of Linv */                                  /* Derive the diagonal block of Linv */
493          nnz = parent_inv_ap(ncj, 0, parent, Linvp);          nnz = parent_inv_ap(ncj, 0, parent, Linvp);
494          SET_SLOT(Linvj, Matrix_iSym, allocVector(INTSXP, nnz));          SET_SLOT(Linvj, Matrix_iSym, allocVector(INTSXP, nnz));
495          parent_inv_ai(ncj, 0, parent, INTEGER(GET_SLOT(Linvj, Matrix_iSym)));          parent_inv_ai(ncj, 0, parent, INTEGER(GET_SLOT(Linvj, Matrix_iSym)));
496          SET_SLOT(Linvj, Matrix_xSym, allocVector(REALSXP, nnz));          SET_SLOT(Linvj, Matrix_xSym, alloc3Darray(REALSXP, dims[0], dims[1], nnz));
497          dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));          dtmp = REAL(GET_SLOT(Linvj, Matrix_xSym));
498          for (i = 0; i < nnz; i++) dtmp[i] = 1.;          for (i = 0; i < nnz; i++) dtmp[i] = 1.;
499          for (k = j+1; k < nf; k++) { /* Update other blocks in this column */          for (k = j+1; k < nf; k++) { /* Update other blocks in this column */
500              symbolic_right_unit_mm(Linvj, VECTOR_ELT(ctab, Lind(k,j)));              symbolic_right_unit_mm(Linvj, VECTOR_ELT(ZZpO, Lind(k,j)));
501          }          }
502          for (k = j+1; k < nf; k++) { /* Update remaining columns */          for (k = j+1; k < nf; k++) { /* Update remaining columns */
503              for (i = k; i < nf; i++) block_update(ctab, j, k, i);              for (i = k; i < nf; i++) block_update(ZZpO, j, k, i);
504          }          }
505          for (i= 0; i < j; i++) { /* copy blocks to the left */          for (i= 0; i < j; i++) { /* copy blocks to the left */
506              SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ctab, Lind(j,i))));              SET_VECTOR_ELT(L, Lind(j,i), duplicate(VECTOR_ELT(ZZpO, Lind(j,i))));
507          }          }
508      }      }
509      UNPROTECT(2);      Free(nlev);
     return ans;  
510  }  }
   
   
   
   

Legend:
Removed from v.409  
changed lines
  Added in v.410

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