SCM

SCM Repository

[matrix] Diff of /branches/Matrix-mer2/src/lCholCMatrix.c
ViewVC logotype

Diff of /branches/Matrix-mer2/src/lCholCMatrix.c

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

revision 733, Tue May 17 17:18:30 2005 UTC revision 746, Wed May 18 17:03:17 2005 UTC
# Line 66  Line 66 
66   * where A is an lCholCMatrix object and B and C are lgCMatrix   * where A is an lCholCMatrix object and B and C are lgCMatrix
67   * objects.   * objects.
68   *   *
69     * An early check is done to see if A is the identity, in which case
70     * BIP is returned and bp is unmodified.
71     *
72   * @param side LFT or RGT   * @param side LFT or RGT
73   * @param transa TRN or NTR   * @param transa TRN or NTR
74   * @param m number of rows in B and C   * @param m number of rows in B and C
75   * @param n number of columns in B and C   * @param n number of columns in B and C
76   * @param Parent Parent array of A   * @param Parent Parent array of A
77   * @param bi array of row indices of B   * @param BIP pointer to the row indices for B
78   * @param bp array of column pointers of B   * @param bp array of column pointers for B (may be overwritten)
  * @param CIP pointer to the row indices for C  
  * @param cp array of column pointers for C  
79   *   *
80   * @return Pointer to the updated row indices for C.  The column   * @return Pointer to the updated row indices for C.  The column
81   * pointers may also be overwritten.   * pointers bp are overwritten with cp.
82   */   */
83  SEXP  SEXP
84  lCholClgCsm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m,  lCholClgCsm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m,
85              int n, const int Parent[], const int bi[], const int bp[],              int n, const int Parent[], SEXP BIP, int bp[])
             SEXP CIP, int cp[])  
86  {  {
87      int *ci = INTEGER(CIP), cnz, extra, i, j, pari, pos;      int *bi = INTEGER(BIP), bnz, extra, ident, j, pari, pos;
88        int nca = (transa == TRN) ? n : m;
89    
90        ident = 1;
91        for (j = 0; j < nca; j++)
92            if (Parent[j] >= 0) {
93                ident = 0;
94                break;
95            }
96        if (ident) return BIP;
97    
98      extra = 0;      extra = 0;
99      if (side == LFT) {      if (side == LFT) {
# Line 96  Line 105 
105              for (j = 0; j < n; j++) {              for (j = 0; j < n; j++) {
106                  int ii, ii2 = bp[j + 1];                  int ii, ii2 = bp[j + 1];
107                  for (ii = bp[j]; ii < ii2; ii++)                  for (ii = bp[j]; ii < ii2; ii++)
108                      for (pari = bi[ii]; pari >= 0; pari = Parent[pari])                      for (pari = Parent[bi[ii]]; pari >= 0;
109                          if (check_csc_index(cp, ci, pari, j, -1) < 0) extra++;                           pari = Parent[pari]) extra++;
110              }              }
 /* This is not quite right.  C may contain more elements than in the solution. */  
             if (!extra) return CIP;  
111    
112              cnz = cp[n];              bnz = bp[n];
113              ntot = cnz + extra;              ntot = bnz + extra;
114              Ti = Memcpy(Calloc(ntot, int), ci, cnz);              Ti = Memcpy(Calloc(ntot, int), bi, bnz);
115              Tj = expand_cmprPt(n, cp, Calloc(ntot, int));              Tj = expand_cmprPt(n, bp, Calloc(ntot, int));
116              Tci = Calloc(ntot, int);              Tci = Calloc(ntot, int);
117    
118              pos = cnz;              pos = bnz;
119              for (j = 0; j < n; j++) {              for (j = 0; j < n; j++) {
120                  int ii, ii2 = bp[j + 1];                  int ii, ii2 = bp[j + 1];
121                  for (ii = bp[j]; ii < ii2; ii++)                  for (ii = bp[j]; ii < ii2; ii++)
122                      for (pari = bi[ii]; pari >= 0; pari = Parent[pari])                      for (pari = Parent[bi[ii]]; pari >= 0;
123                          if (check_csc_index(cp, ci, pari, j, -1) < 0) {                           pari = Parent[pari]) {
124                              Ti[pos] = pari;                              Ti[pos] = pari;
125                              Tj[pos] = j;                              Tj[pos] = j;
126                              pos++;                              pos++;
127                          }                          }
128              }              }
129              triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,              triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
130                             cp, Tci, (double *) NULL);                             bp, Tci, (double *) NULL);
131    
132              cnz = cp[n];              bnz = bp[n];
133              CIP = PROTECT(allocVector(INTSXP, cnz));              BIP = PROTECT(allocVector(INTSXP, bnz));
134              Memcpy(INTEGER(CIP), Tci, cnz);              Memcpy(INTEGER(BIP), Tci, bnz);
135    
136              Free(Tci); Free(Ti); Free(Tj);              Free(Tci); Free(Ti); Free(Tj);
137              UNPROTECT(1);              UNPROTECT(1);
138              return CIP;              return BIP;
139          }          }
140      } else {      } else {
141          if (transa == TRN) {          if (transa == TRN) {
142              int *Tci, *Tj, *Ti, ntot;              int *Tci, *Tj, *Ti, ntot;
143              for (j = 0; j < n; j++) {              for (j = 0; j < n; j++) {
144                  int ii, ii2 = bp[j + 1];                  int ii, ii2 = bp[j + 1];
145                  for (pari = j; pari >= 0; pari = Parent[pari]) {                  for (pari = Parent[j]; pari >= 0; pari = Parent[pari])
146                      for (ii = bp[j]; ii < ii2; ii++)                      for (ii = bp[j]; ii < ii2; ii++) extra++;
                         if (check_csc_index(cp, ci, ii, pari, -1) < 0) extra++;  
147              }              }
 /* This is not quite right.  C may contain more elements than in the solution. */  
             if (!extra) return CIP;  
148    
149              cnz = cp[n];              bnz = bp[n];
150              ntot = cnz + extra;              ntot = bnz + extra;
151              Ti = Memcpy(Calloc(ntot, int), ci, cnz);              Ti = Memcpy(Calloc(ntot, int), bi, bnz);
152              Tj = expand_cmprPt(n, cp, Calloc(ntot, int));              Tj = expand_cmprPt(n, bp, Calloc(ntot, int));
153              Tci = Calloc(ntot, int);              Tci = Calloc(ntot, int);
154    
155              pos = cnz;              pos = bnz;
156              for (j = 0; j < n; j++) {              for (j = 0; j < n; j++) {
157                  int ii, ii2 = bp[j + 1];                  int ii, ii2 = bp[j + 1];
158                  for (pari = j; pari >= 0; pari = Parent[pari]) {                  for (pari = Parent[j]; pari >= 0; pari = Parent[pari]) {
159                      for (ii = bp[j]; ii < ii2; ii++)                      for (ii = bp[j]; ii < ii2; ii++) {
160                          if (check_csc_index(cp, ci, ii, pari, -1) < 0) {                          Ti[pos] = bi[ii];
                             Ti[pos] = ii;  
161                              Tj[pos] = pari;                              Tj[pos] = pari;
162                              pos++;                              pos++;
163                          }                          }
164              }              }
165                }
166              triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,              triplet_to_col(m, n, ntot, Ti, Tj, (double *) NULL,
167                             cp, Tci, (double *) NULL);                             bp, Tci, (double *) NULL);
168                bnz = bp[n];
169              cnz = cp[n];              BIP = PROTECT(allocVector(INTSXP, bnz));
170              CIP = PROTECT(allocVector(INTSXP, cnz));              Memcpy(INTEGER(BIP), Tci, bnz);
             Memcpy(INTEGER(CIP), Tci, cnz);  
171    
172              Free(Tci); Free(Ti); Free(Tj);              Free(Tci); Free(Ti); Free(Tj);
173              UNPROTECT(1);              UNPROTECT(1);
174              return CIP;              return BIP;
175          } else {          } else {
176              error(_("code not yet written"));              error(_("code not yet written"));
177              return R_NilValue;              return R_NilValue;
# Line 188  Line 191 
191      SET_SLOT(ans, Matrix_iSym,      SET_SLOT(ans, Matrix_iSym,
192               lCholClgCsm(LFT, NTR, bdims[0], bdims[1],               lCholClgCsm(LFT, NTR, bdims[0], bdims[1],
193                           INTEGER(GET_SLOT(a, Matrix_ParentSym)),                           INTEGER(GET_SLOT(a, Matrix_ParentSym)),
                          INTEGER(GET_SLOT(b, Matrix_iSym)),  
                          INTEGER(GET_SLOT(b, Matrix_pSym)),  
194                           GET_SLOT(ans, Matrix_iSym),                           GET_SLOT(ans, Matrix_iSym),
195                           INTEGER(GET_SLOT(ans, Matrix_pSym))));                           INTEGER(GET_SLOT(ans, Matrix_pSym))));
196      UNPROTECT(1);      UNPROTECT(1);

Legend:
Removed from v.733  
changed lines
  Added in v.746

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