SCM

SCM Repository

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

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

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

revision 1251, Sat Apr 15 13:08:26 2006 UTC revision 1262, Sun Apr 30 15:40:39 2006 UTC
# Line 260  Line 260 
260      UNPROTECT(1);      UNPROTECT(1);
261      return ans;      return ans;
262  }  }
263    
264    SEXP dtCMatrix_upper_solve(SEXP a)
265    {
266        SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
267        int lo = uplo_P(a)[0] == 'L', unit = diag_P(a)[0] == 'U',
268            n = INTEGER(GET_SLOT(a, Matrix_DimSym))[0], nnz,
269            *ai = INTEGER(GET_SLOT(a,Matrix_iSym)),
270            *ap = INTEGER(GET_SLOT(a, Matrix_pSym)), *bi,
271            *bp = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1));
272        int bnz = 10 * ap[n];         /* initial estimate of nnz in b */
273        int *ti = Calloc(bnz, int), j, nz;
274        double *ax = REAL(GET_SLOT(a, Matrix_xSym)), *tx = Calloc(bnz, double),
275            *tmp = Calloc(n, double);
276    
277        if (lo || (!unit)) error(_("Code written for unit upper triangular unit matrices"));
278        bp[0] = 0;
279        for (j = 0; j < n; j++) {
280            int i, i1 = ap[j + 1];
281            AZERO(tmp, n);
282            for (i = ap[j]; i < i1; i++) {
283                int ii = ai[i], k;
284                if (unit) tmp[ii] -= ax[i];
285                for (k = bp[ii]; k < bp[ii + 1]; k++) tmp[ti[k]] -= ax[i] * tx[k];
286            }
287            for (i = 0, nz = 0; i < n; i++) if (tmp[i]) nz++;
288            bp[j + 1] = bp[j] + nz;
289            if (bp[j + 1] > bnz) {
290                while (bp[j + 1] > bnz) bnz *= 2;
291                ti = Realloc(ti, bnz, int);
292                tx = Realloc(tx, bnz, double);
293            }
294            i1 = bp[j];
295            for (i = 0; i < n; i++) if (tmp[i]) {ti[i1] = i; tx[i1] = tmp[i]; i1++;}
296        }
297        nz = bp[n];
298        Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)), ti, nz);
299        Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)), tx, nz);
300        Free(tmp); Free(tx); Free(ti);
301        SET_SLOT(ans, Matrix_DimSym, duplicate(GET_SLOT(a, Matrix_DimSym)));
302        SET_SLOT(ans, Matrix_DimNamesSym, duplicate(GET_SLOT(a, Matrix_DimNamesSym)));
303        SET_SLOT(ans, Matrix_uploSym, duplicate(GET_SLOT(a, Matrix_uploSym)));
304        SET_SLOT(ans, Matrix_diagSym, duplicate(GET_SLOT(a, Matrix_diagSym)));
305        UNPROTECT(1);
306        return ans;
307    }

Legend:
Removed from v.1251  
changed lines
  Added in v.1262

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