SCM

SCM Repository

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

Diff of /pkg/Matrix/src/dense.c

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

revision 3019, Sat Oct 11 20:51:53 2014 UTC revision 3020, Tue Oct 14 16:14:02 2014 UTC
# Line 340  Line 340 
340  SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)  SEXP dense_to_symmetric(SEXP x, SEXP uplo, SEXP symm_test)
341  /* Class of result will be [dln]syMatrix */  /* Class of result will be [dln]syMatrix */
342  {  {
343    /*== FIXME: allow  uplo = NA   and then behave a bit like symmpart():
344     *== -----  would use the *dimnames* to determine U or L   (??)
345     */
346    
347      int symm_tst = asLogical(symm_test);      int symm_tst = asLogical(symm_test);
348      SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));      SEXP dx = PROTECT(dup_mMatrix_as_geMatrix(x));
349      SEXP ans, dns;      SEXP ans, dns, nms_dns;
350      const char *cl = class_P(dx);      const char *cl = class_P(dx);
351      /* same as in ..._geMatrix() above:*/      /* same as in ..._geMatrix() above:*/
352      enum dense_enum { ddense, ldense, ndense      enum dense_enum { ddense, ldense, ndense
# Line 376  Line 380 
380              CHECK_SYMMETRIC              CHECK_SYMMETRIC
381          }          }
382      }      }
383    #   undef CHECK_SYMMETRIC
384    
385        ans = PROTECT(NEW_OBJECT(MAKE_CLASS( M_type == ddense ? "dsyMatrix" :
386                                            (M_type == ldense ? "lsyMatrix" :
387                                             "nsyMatrix"))));
388    
389    
390    // --- FIXME: Use MK_SYMMETRIC_DIMNAMES_AND_RETURN  from below -- with "uplo" argument
391    
392        /* need _symmetric_ dimnames */
393      dns = GET_SLOT(dx, Matrix_DimNamesSym);      dns = GET_SLOT(dx, Matrix_DimNamesSym);
394      if(!equal_string_vectors(VECTOR_ELT(dns,0),      if(!equal_string_vectors(VECTOR_ELT(dns,0),
395                               VECTOR_ELT(dns,1))) {                               VECTOR_ELT(dns,1))) {
         /* need _symmetric_ dimnames */  
396          if(*CHAR(asChar(uplo)) == 'U')          if(*CHAR(asChar(uplo)) == 'U')
397              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));
398          else          else
399              SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));              SET_VECTOR_ELT(dns,1, VECTOR_ELT(dns,0));
400      }      }
401        if(!isNull(nms_dns = getAttrib(dns, R_NamesSymbol)) &&
402           !R_compute_identical(STRING_ELT(nms_dns, 0),
403                                STRING_ELT(nms_dns, 1), 15)) { // names(dimnames(.)) :
404            if(*CHAR(asChar(uplo)) == 'U')
405                SET_STRING_ELT(nms_dns, 0, STRING_ELT(nms_dns,1));
406            else
407                SET_STRING_ELT(nms_dns, 1, STRING_ELT(nms_dns,0));
408            setAttrib(dns, R_NamesSymbol, nms_dns);
409        }
410    
411      /* Copy dx to ans;      /* Copy dx to ans;
412       * Because slots of dx are freshly allocated and dx will not be       * Because slots of dx are freshly allocated and dx will not be
413       * used, we use the slots themselves and don't duplicate */       * used, we use the slots themselves and don't duplicate */
     ans = PROTECT(NEW_OBJECT(MAKE_CLASS( M_type == ddense ? "dsyMatrix" :  
                                         (M_type == ldense ? "lsyMatrix" :  
                                          "nsyMatrix"))));  
414      SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));      SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));
415      SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));      SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));
416      SET_SLOT(ans, Matrix_DimNamesSym, dns);      SET_SLOT(ans, Matrix_DimNamesSym, dns);
# Line 405  Line 423 
423  SEXP ddense_symmpart(SEXP x)  SEXP ddense_symmpart(SEXP x)
424  /* Class of the value will be dsyMatrix */  /* Class of the value will be dsyMatrix */
425  {  {
426      SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x));      SEXP dx = dup_mMatrix_as_dgeMatrix(x);
427      int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];      int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
428    
429      if(n != adims[1]) {      if(n != adims[1]) {
         UNPROTECT(1);  
430          error(_("matrix is not square! (symmetric part)"));          error(_("matrix is not square! (symmetric part)"));
431          return R_NilValue; /* -Wall */          return R_NilValue; /* -Wall */
432      } else {      } else {
433          SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns;          PROTECT(dx);
434            SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsyMatrix"))), dns, nms_dns;
435          double *xx = REAL(GET_SLOT(dx, Matrix_xSym));          double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
         int i,j;  
436    
437          /* only need to assign the *upper* triangle (uplo = "U");          /* only need to assign the *upper* triangle (uplo = "U");
438           * noting that diagonal remains unchanged */           * noting that diagonal remains unchanged */
439          for (j = 0; j < n; j++) {          for (int j = 0; j < n; j++) {
440              for (i = 0; i < j; i++) {              for (int i = 0; i < j; i++) {
441                  xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;                  xx[j * n + i] = (xx[j * n + i] + xx[i * n + j]) / 2.;
442              }              }
443          }          }
444    
445          dns = GET_SLOT(dx, Matrix_DimNamesSym);  #       define MK_SYMMETRIC_DIMNAMES_AND_RETURN                         \
446          if(!equal_string_vectors(VECTOR_ELT(dns,0),                                                                          \
447                                   VECTOR_ELT(dns,1))) {          dns = GET_SLOT(dx, Matrix_DimNamesSym);                         \
448              /* need _symmetric_ dimnames */          int J = 1;                                                      \
449              SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* ==> uplo = "U" */          if(!equal_string_vectors(VECTOR_ELT(dns,0),                     \
450          }                                   VECTOR_ELT(dns,1))) {                  \
451                /* _symmetric_ dimnames: behave as symmDN(*, col=TRUE) */   \
452          /* Copy dx to ans;              if(isNull(VECTOR_ELT(dns, J)))                              \
453           * Because slots of dx are freshly allocated and dx will not be                  J = !J;                                                 \
454           * used, we use the slots themselves and don't duplicate */              SET_VECTOR_ELT(dns, !J, VECTOR_ELT(dns, J));                \
455            } /* names(dimnames(.)):*/                                      \
456          SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));          if(!isNull(nms_dns = getAttrib(dns, R_NamesSymbol)) &&          \
457          SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));             !R_compute_identical(STRING_ELT(nms_dns, 0),                 \
458          SET_SLOT(ans, Matrix_DimNamesSym, dns);                                  STRING_ELT(nms_dns, 1), 15)) {          \
459          SET_SLOT(ans, Matrix_uploSym,     mkString("U"));              SET_STRING_ELT(nms_dns, !J, STRING_ELT(nms_dns, J));        \
460                setAttrib(dns, R_NamesSymbol, nms_dns);                     \
461            }                                                               \
462                                                                            \
463            /* Copy dx to ans;                                              \
464             * Because slots of dx are freshly allocated and dx will not be \
465             * used, we use the slots themselves and don't duplicate */     \
466                                                                            \
467            SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));   \
468            SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym)); \
469            SET_SLOT(ans, Matrix_DimNamesSym, dns);                         \
470            SET_SLOT(ans, Matrix_uploSym,     mkString("U"));               \
471                                                                            \
472            UNPROTECT(2);                                                   \
473            return ans
474    
475          UNPROTECT(2);          MK_SYMMETRIC_DIMNAMES_AND_RETURN;
         return ans;  
476      }      }
477  }  }
478    
479  SEXP ddense_skewpart(SEXP x)  SEXP ddense_skewpart(SEXP x)
480  /* Class of the value will be dgeMatrix */  /* Class of the value will be dgeMatrix */
481  {  {
482      SEXP dx = PROTECT(dup_mMatrix_as_dgeMatrix(x));      SEXP dx = dup_mMatrix_as_dgeMatrix(x);
483      int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];      int *adims = INTEGER(GET_SLOT(dx, Matrix_DimSym)), n = adims[0];
484    
485      if(n != adims[1]) {      if(n != adims[1]) {
         UNPROTECT(1);  
486          error(_("matrix is not square! (skew-symmetric part)"));          error(_("matrix is not square! (skew-symmetric part)"));
487          return R_NilValue; /* -Wall */          return R_NilValue; /* -Wall */
488      } else {      } else {
489          SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns;          PROTECT(dx);
490            SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dns, nms_dns;
491          double *xx = REAL(GET_SLOT(dx, Matrix_xSym));          double *xx = REAL(GET_SLOT(dx, Matrix_xSym));
         int i,j;  
492    
493          for (j = 0; j < n; j++) {          for (int j = 0; j < n; j++) {
494              xx[j * n + j] = 0.;              xx[j * n + j] = 0.;
495              for (i = 0; i < j; i++) {              for (int i = 0; i < j; i++) {
496                  double s = (xx[j * n + i] - xx[i * n + j]) / 2.;                  double s = (xx[j * n + i] - xx[i * n + j]) / 2.;
497                  xx[j * n + i] =  s;                  xx[j * n + i] =  s;
498                  xx[i * n + j] = -s;                  xx[i * n + j] = -s;
499              }              }
500          }          }
501    
502          dns = GET_SLOT(dx, Matrix_DimNamesSym);          MK_SYMMETRIC_DIMNAMES_AND_RETURN;
         if(!equal_string_vectors(VECTOR_ELT(dns,0),  
                                  VECTOR_ELT(dns,1))) {  
             /* need _symmetric_ dimnames */  
             SET_VECTOR_ELT(dns,0, VECTOR_ELT(dns,1));/* uplo = "U" */  
         }  
   
         /* Copy dx to ans;  
          * Because slots of dx are freshly allocated and dx will not be  
          * used, we use the slots themselves and don't duplicate */  
   
         SET_SLOT(ans, Matrix_xSym,        GET_SLOT(dx, Matrix_xSym));  
         SET_SLOT(ans, Matrix_DimSym,      GET_SLOT(dx, Matrix_DimSym));  
         SET_SLOT(ans, Matrix_DimNamesSym, dns);  
         SET_SLOT(ans, Matrix_uploSym,     mkString("U"));  
   
         UNPROTECT(2);  
         return ans;  
503      }      }
504  }  }

Legend:
Removed from v.3019  
changed lines
  Added in v.3020

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