SCM

SCM Repository

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

Diff of /pkg/src/Tsparse.c

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

revision 2119, Tue Mar 4 21:44:04 2008 UTC revision 2120, Tue Mar 4 21:44:41 2008 UTC
# Line 59  Line 59 
59                                /* diag = */ CHAR(STRING_ELT(diag, 0)),                                /* diag = */ CHAR(STRING_ELT(diag, 0)),
60                                GET_SLOT(x, Matrix_DimNamesSym));                                GET_SLOT(x, Matrix_DimNamesSym));
61  }  }
62    
63    SEXP Tsparse_diagU2N(SEXP x)
64    {
65        char *valid[] = {"dtTMatrix", /* 0 */
66                         "ltTMatrix", /* 1 */
67                         "ntTMatrix", /* 2 : no x slot */
68                         "ztTMatrix", /* 3 */
69                         ""};
70    /* #define xSXP(iTyp) ((iTyp == 0) ? REALSXP : ((iTyp == 1) ? LGLSXP : /\* else *\/ CPLXSXP)); */
71    /* #define xTYPE(iTyp) ((iTyp == 0) ? double : ((iTyp == 1) ? int : /\* else *\/ Rcomplex)); */
72        int ctype = Matrix_check_class(class_P(x), valid);
73    
74        if (ctype < 0 || *diag_P(x) != 'U') {
75            /* "trivially fast" when not triangular (<==> no 'diag' slot),
76               or not *unit* triangular */
77            return (x);
78        }
79        else { /* instead of going to Csparse -> Cholmod -> Csparse -> Tsparse, work directly: */
80            int i, n = INTEGER(GET_SLOT(x, Matrix_DimSym))[0],
81                nnz = length(GET_SLOT(x, Matrix_iSym)),
82                new_n = nnz + n;
83            SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS(class_P(x))));
84            int *islot = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, new_n)),
85                *jslot = INTEGER(ALLOC_SLOT(ans, Matrix_jSym, INTSXP, new_n));
86    
87            slot_dup(ans, x, Matrix_DimSym);
88            SET_DimNames(ans, x);
89            slot_dup(ans, x, Matrix_uploSym);
90            SET_SLOT(ans, Matrix_diagSym, mkString("N"));
91    
92            /* Build the new i- and j- slots : first copy the current : */
93            Memcpy(islot, INTEGER(GET_SLOT(x, Matrix_iSym)), nnz);
94            Memcpy(jslot, INTEGER(GET_SLOT(x, Matrix_jSym)), nnz);
95            /* then, add the new (i,j) slot entries: */
96            for(i = 0; i < n; i++) {
97                islot[i + nnz] = i;
98                jslot[i + nnz] = i;
99            }
100    
101            /* build the new x-slot : */
102            switch(ctype) {
103            case 0: { /* "d" */
104                double *x_new = REAL(ALLOC_SLOT(ans, Matrix_xSym,
105                                                REALSXP, new_n));
106                Memcpy(x_new, REAL(GET_SLOT(x, Matrix_xSym)), nnz);
107                for(i = 0; i < n; i++) /* add  x[i,i] = 1. */
108                    x_new[i + nnz] = 1.;
109                break;
110            }
111            case 1: { /* "l" */
112                int *x_new = LOGICAL(ALLOC_SLOT(ans, Matrix_xSym,
113                                                LGLSXP, new_n));
114                Memcpy(x_new, LOGICAL(GET_SLOT(x, Matrix_xSym)), nnz);
115                for(i = 0; i < n; i++) /* add  x[i,i] = 1 (= TRUE) */
116                    x_new[i + nnz] = 1;
117                break;
118            }
119            case 2: /* "n" */
120                    /* nothing to do here */
121                break;
122    
123            case 3: { /* "z" */
124                Rcomplex *x_new = COMPLEX(ALLOC_SLOT(ans, Matrix_xSym,
125                                                     CPLXSXP, new_n));
126                Memcpy(x_new, COMPLEX(GET_SLOT(x, Matrix_xSym)), nnz);
127                for(i = 0; i < n; i++) /* add  x[i,i] = 1 (= TRUE) */
128                    x_new[i + nnz] = (Rcomplex) {1., 0.};
129                break;
130            }
131    
132            }/* switch() */
133    
134            UNPROTECT(1);
135            return ans;
136        }
137    }

Legend:
Removed from v.2119  
changed lines
  Added in v.2120

root@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