# SCM Repository

[matrix] Diff of /pkg/src/dgBCMatrix.c
 [matrix] / pkg / src / dgBCMatrix.c

# Diff of /pkg/src/dgBCMatrix.c

revision 448, Fri Jan 21 21:59:45 2005 UTC revision 451, Sun Jan 23 10:05:53 2005 UTC
# Line 177  Line 177
177
178      if (cdims[0] != cdims[1]) error("blocks in C must be square");      if (cdims[0] != cdims[1]) error("blocks in C must be square");
179      if (trans == TRN) {      if (trans == TRN) {
error("Code for trans == 'T' not yet written");
180                                  /* FIXME: Write this part */                                  /* FIXME: Write this part */
181            error("Code for trans == 'T' not yet written");
182      } else {      } else {
183          if (adims[0] != cdims[0])          if (adims[0] != cdims[0])
184              error("Inconsistent dimensions: A[%d,%d,%d], C[%d,%d,%d]",              error("Inconsistent dimensions: A[%d,%d,%d], C[%d,%d,%d]",
# Line 221  Line 221
221  }  }
222
223  /**  /**
224   * Create the LDL' decomposition of the positive definite symmetric   * Create the LD^{T/2}D^{1/2}L' decomposition of the positive definite
225   * cscBlocked matrix A (upper triangle stored) in L and D.   * symmetric cscBlocked matrix A (upper triangle stored) in L and
226     * D^{1/2}.  The notation D^{1/2} denotes the upper Cholesky factor of
227     * the positive definite positive definite block diagonal matrix D.
228     * The diagonal blocks are of size nci.
229   *   *
230   * @param A pointer to a cscBlocked object containing the upper   * @param A pointer to a cscBlocked object containing the upper
231   * triangle of a positive definite symmetric matrix.   * triangle of a positive definite symmetric matrix.
# Line 314  Line 317
317                                      &one, Y + Li[p]*ncisqr, &nci);                                      &one, Y + Li[p]*ncisqr, &nci);
318                  }                  }
319                  /* FIXME: Check that this is the correct order and transposition */                  /* FIXME: Check that this is the correct order and transposition */
320                  F77_CALL(dtrsm)("L", "U", "N", "N", &nci, &nci,                  F77_CALL(dtrsm)("R", "U", "N", "N", &nci, &nci,
321                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
322                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,
323                                  &one, Dx + k*ncisqr, &nci);                                  &one, Dx + k*ncisqr, &nci);
324                  F77_CALL(dtrsm)("L", "U", "T", "N", &nci, &nci,                  F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,
325                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
326                  Li[p] = k;      /* store L[k,i] in column form of L */                  Li[p] = k;      /* store L[k,i] in column form of L */
327                  Memcpy(Lx + p * ncisqr, Yi, ncisqr);                  Memcpy(Lx + p * ncisqr, Yi, ncisqr);
# Line 355  Line 358
358            enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,            enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
359            double alpha, SEXP A, double B[], int m, int n, int ldb)            double alpha, SEXP A, double B[], int m, int n, int ldb)
360  {  {
361      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP /* ApP = GET_SLOT(A, Matrix_pSym), */
362          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
363      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int /* *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)), */
364          *Ap = INTEGER(ApP),  /*      *Ap = INTEGER(ApP), */
365          *xdims = INTEGER(getAttrib(AxP, R_DimSymbol)),          *xdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
366          i, j, nb = length(ApP) - 1;          i, j/* , nb = length(ApP) - 1 */;
367
368      if (xdims[0] != xdims[1])      if (xdims[0] != xdims[1])
369          error("Argument A to cscb_trmm is not triangular");          error("Argument A to cscb_trmm is not triangular");
# Line 479  Line 482
482             enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,             enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
483             double alpha, SEXP A, SEXP B)             double alpha, SEXP A, SEXP B)
484  {  {
485      int      SEXP
486          iunit = (diag == UNT);  /*      ApP = GET_SLOT(A, Matrix_pSym), */
SEXP ApP = GET_SLOT(A, Matrix_pSym),
487          AxP = GET_SLOT(A, Matrix_xSym),          AxP = GET_SLOT(A, Matrix_xSym),
488          BpP = GET_SLOT(B, Matrix_pSym),  /*      , BpP = GET_SLOT(B, Matrix_pSym) */
489          BxP = GET_SLOT(B, Matrix_xSym);          BxP = GET_SLOT(B, Matrix_xSym);
490      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int
491          *Ap = INTEGER(ApP),  /*      *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)), */
492          *Bi = INTEGER(GET_SLOT(B, Matrix_iSym)),  /*      *Ap = INTEGER(ApP), */
493          *Bp = INTEGER(BpP),  /*      *Bi = INTEGER(GET_SLOT(B, Matrix_iSym)), */
494    /*      *Bp = INTEGER(BpP), */
495          *axdims = INTEGER(getAttrib(AxP, R_DimSymbol)),          *axdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
496          *bxdims = INTEGER(getAttrib(BxP, R_DimSymbol)),          *bxdims = INTEGER(getAttrib(BxP, R_DimSymbol))
497          ncbA = length(ApP) - 1,  /*      , ncbA = length(ApP) - 1 */
498          ncbB = length(BpP) - 1;  /*      , ncbB = length(BpP) - 1 */
499            ;
500      int i, nbx = bxdims[0] * bxdims[1] * bxdims[2];      int i, nbx = bxdims[0] * bxdims[1] * bxdims[2];
501
502      if (axdims[0] != axdims[1])      if (axdims[0] != axdims[1])
# Line 556  Line 560
560          *Bp = INTEGER(BpP),          *Bp = INTEGER(BpP),
561          *axdims = INTEGER(getAttrib(AxP, R_DimSymbol)),          *axdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
562          *bxdims = INTEGER(getAttrib(BxP, R_DimSymbol)),          *bxdims = INTEGER(getAttrib(BxP, R_DimSymbol)),
563          ncbA = length(ApP) - 1,  /*      ncbA = length(ApP) - 1, */
564          ncbB = length(BpP) - 1;          ncbB = length(BpP) - 1;
565      int i, j, nbx = bxdims[0] * bxdims[1] * bxdims[2];      int i, j, nbx = bxdims[0] * bxdims[1] * bxdims[2];
566      double *Ax = REAL(AxP), *Bx = REAL(BxP);      double *Ax = REAL(AxP), *Bx = REAL(BxP);
# Line 672  Line 676
676      }      }
677      error("Code not yet written");      error("Code not yet written");
678  }  }
679
680    SEXP cscBlocked_2cscMatrix(SEXP A)
681    {
682        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("cscMatrix"))),
683            ApP = GET_SLOT(A, Matrix_pSym),
684            AiP = GET_SLOT(A, Matrix_iSym),
685            AxP = GET_SLOT(A, Matrix_xSym);
686        int *Ai = INTEGER(AiP), *Ap = INTEGER(ApP), *Bi, *Bp, *Dims,
687            *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
688            ii, j, ncb = length(ApP) - 1, nnz, nrb;
689        int nc = adims[1], nr = adims[0];
690        int sz = nc * nr;
691        double *Ax = REAL(AxP), *Bx;
692
693        SET_SLOT(val, Matrix_factorization, allocVector(VECSXP, 0));
694        SET_SLOT(val, Matrix_DimSym, allocVector(INTSXP, 2));
695        Dims = INTEGER(GET_SLOT(val, Matrix_DimSym));
696        Dims[1] = ncb * adims[1];
697                                    /* find number of row blocks */
698        for (j = 0, nrb = -1; j < adims[2]; j++) if (Ai[j] > nrb) nrb = Ai[j];
699        Dims[0] = nrb * adims[0];
700        nnz = length(AxP);
701
702        if (nc == 1) {              /* x slot is in the correct order */
703            SET_SLOT(val, Matrix_pSym, duplicate(ApP));
704            SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
705            SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
706            Memcpy(REAL(GET_SLOT(val, Matrix_xSym)), Ax, nnz);
707            if (nr == 1) {
708                Memcpy(INTEGER(GET_SLOT(val, Matrix_iSym)), Ai, nnz);
709            } else {
710                Bi = INTEGER(GET_SLOT(val, Matrix_iSym));
711                Bp = INTEGER(GET_SLOT(val, Matrix_pSym));
712                for (j = 0; j <= ncb; j++) Bp[j] *= nr;
713                for (j = 0; j < adims[2]; j++) {
714                    for (ii = 0; ii < nr; ii++) {
715                        Bi[j * nr + ii] = Ai[j] * nr + ii;
716                    }
717                }
718            }
719        } else {
720            SET_SLOT(val, Matrix_pSym, allocVector(INTSXP, Dims[1] + 1));
721            Bp = INTEGER(GET_SLOT(val, Matrix_pSym));
722            SET_SLOT(val, Matrix_iSym, allocVector(INTSXP, nnz));
723            Bi = INTEGER(GET_SLOT(val, Matrix_iSym));
724            SET_SLOT(val, Matrix_xSym, allocVector(REALSXP, nnz));
725            Bx = REAL(GET_SLOT(val, Matrix_xSym));
726
727            Bp[0] = 0;
728            for (j = 0; j < ncb; j++) { /* Column blocks of A */
729                int i, i1 = Ap[j], i2 = Ap[j + 1], jj;
730                int nzbc = (i2 - i1) * nr; /* No. of non-zeroes in B column */
731
732                for (jj = 0; jj < nc; jj++) { /* column within blocks */
733                    int jb = nc * j + jj; /* Column number in B */
734
735                    Bp[jb] = i1 * sz + jj * nzbc;
736                    for (i = i1; i < i2; i++) { /* index in Ai and Ax */
737                        for (ii = 0; ii < adims[0]; ii++) { /* row within blocks */
738                            int ind = ii + (i - i1) * nr + Bp[jb];
739
740                            Bi[ind] = Ai[i] * sz + jj * nzbc + ii;
741                            Bx[ind] = Ax[i * sz + jj * nc + ii];
742                        }
743                    }
744                }
745            }
746        }
747        UNPROTECT(1);
748        return val;
749    }

Legend:
 Removed from v.448 changed lines Added in v.451

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: