SCM

SCM Repository

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

Diff of /pkg/src/cscBlocked.c

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

revision 422, Tue Jan 4 23:45:15 2005 UTC revision 434, Thu Jan 13 19:36:12 2005 UTC
# Line 26  Line 26 
26      return ScalarLogical(1);      return ScalarLogical(1);
27  }  }
28    
29    static int*
30    expand_column_pointers(int ncol, const int mp[], int mj[])
31    {
32        int j;
33        for (j = 0; j < ncol; j++) {
34            int j2 = mp[j+1], jj;
35            for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;
36        }
37        return mj;
38    }
39    
40    static R_INLINE
41    int Tind(const int rowind[], const int colptr[], int i, int j)
42    {
43        int k, k2 = colptr[j + 1];
44        for (k = colptr[j]; k < k2; k++)
45            if (rowind[k] == i) return k;
46        error("row %d and column %d not defined in rowind and colptr",
47              i, j);
48        return -1;                  /* -Wall */
49    }
50    
51  /**  /**
52   * Perform one of the matrix operations   * Perform one of the matrix operations
53   *  C := alpha*op(A)*B + beta*C   *  C := alpha*op(A)*B + beta*C
# Line 289  Line 311 
311          j, k;          j, k;
312      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;
313      int scalar = (adims[0] == 1 && adims[1] == 1),      int scalar = (adims[0] == 1 && adims[1] == 1),
314            anc = length(ApP) - 1,
315          asz = adims[0] * adims[1],          asz = adims[0] * adims[1],
316          csz = cdims[0] * cdims[1];          csz = cdims[0] * cdims[1];
317    
# Line 313  Line 336 
336          if (beta != 1.)          if (beta != 1.)
337              for (j = 0; j < csz * cdims[2]; j++) Cx[j] *= beta;              for (j = 0; j < csz * cdims[2]; j++) Cx[j] *= beta;
338                                  /* individual products */                                  /* individual products */
339          for (j = 0; j < adims[2]; j++) {          for (j = 0; j < anc; j++) {
340              int k, kk, k2 = Ap[j+1];              int k, kk, k2 = Ap[j+1];
341              for (k = Ap[j]; k < k2; k++) {              for (k = Ap[j]; k < k2; k++) {
342                  int ii = Ai[k], K = Ind(ii, ii, Cp, Ci);                  int ii = Ai[k], K = Ind(ii, ii, Cp, Ci);
# Line 357  Line 380 
380      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
381          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
382      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
383          diag, j, n = length(ApP) - 1;          diag, info, j, k, n = length(ApP) - 1;
384      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
385          *Ap = INTEGER(ApP),          *Ap = INTEGER(ApP),
386          *Li = INTEGER(GET_SLOT(L, Matrix_iSym)),          *Li = INTEGER(GET_SLOT(L, Matrix_iSym)),
387          *Lp = INTEGER(GET_SLOT(L, Matrix_pSym));          *Lp = INTEGER(GET_SLOT(L, Matrix_pSym)), nci = adims[0];
388      double *Lx = REAL(GET_SLOT(L, Matrix_xSym)),      double *Lx = REAL(GET_SLOT(L, Matrix_xSym)),
389          *Ax = REAL(AxP), *Dx = REAL(D);          *Ax = REAL(AxP), *Dx = REAL(D), minus1 = -1., one = 1.;
390    
391        if (adims[1] != nci || nci < 1)
392            error("cscb_ldl: dim(A) is [%d, %d, %d]", adims[0], adims[1], adims[2]);
393      for (j = 0, diag = 1; j < n; j++) { /* check for trivial structure */      for (j = 0, diag = 1; j < n; j++) { /* check for trivial structure */
394          if (Parent[j] >= 0) {diag = 0; break;}          if (Parent[j] >= 0) {diag = 0; break;}
395      }      }
396      if (diag) {      if (diag) {
397          Memcpy(REAL(D), Ax, adims[0] * adims[1] * adims[2]);          int ncisqr = nci * nci;
398            Memcpy(Dx, Ax, ncisqr * n);
399            for (j = 0; j < n; j++) { /* form D_i^{1/2} */
400                F77_CALL(dpotrf)("U", &nci, Dx + j * ncisqr, &nci, &k);
401                if (k) error("D[ , , %d], level %d, is not positive definite",
402                             j + 1);
403            }
404          return n;          return n;
405      }      }
406      if (adims[0] == 1 && adims[1] == 1) {      if (nci == 1) {
407          return R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent, Li, Lx, Dx,          k = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent, Li, Lx, Dx,
408                               (int *) NULL, (int *) NULL);                               (int *) NULL, (int *) NULL);
409      } else {          if (k < n) error("cscb_ldl: error code %k from R_ldl_numeric", k);
410            for (j = 0; j < n; j++) Dx[j] = sqrt(Dx[j]);
411          error("code for nontrivial blocked L not yet written");          return n;
412        } else {               /* Copy of ldl_numeric from the LDL package
413                                * modified for blocked sparse matrices */
414            int i, k, p, p2, len, nci = adims[0], ncisqr = adims[0]*adims[0], top;
415            int *Lnz = Calloc(n, int),
416                *Pattern = Calloc(n, int),
417                *Flag = Calloc(n, int);
418            double *Y = Calloc(n * ncisqr, double), *Yi = Calloc(ncisqr, double);
419    
420            for (k = 0; k < n; k++) {
421                /* compute nonzero Pattern of kth row of L, in topological order */
422                AZERO(Y + k*ncisqr, ncisqr); /* Y[,,0:k] is now all zero */
423                top = n;            /* stack for pattern is empty */
424                Flag[k] = k;        /* mark node k as visited */
425                Lnz[k] = 0;         /* count of nonzeros in column k of L */
426                p2 = Ap[k+1];
427                for (p = Ap[k]; p < p2; p++) {
428                    i = Ai[p];      /* get A[i,k] */
429                    if (i <= k) {   /* [i,k] in upper triangle? Should always be true */
430                                    /* copy A(i,k) into Y */
431                        Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);
432                        /* follow path from i to root of etree, stop at flagged node */
433                        for (len = 0; Flag[i] != k; i = Parent[i]) {
434                            Pattern[len++] = i; /* L[k,i] is nonzero */
435                            Flag[i] = k; /* mark i as visited */
436                        }
437                        while (len > 0) { /* push path on top of stack */
438                            Pattern[--top] = Pattern[--len];
439                        }
440                    }
441                }
442                /* Pattern [top ... n-1] now contains nonzero pattern of L[,k] */
443                /* compute numerical values in kth row of L (a sparse triangular solve) */
444                Memcpy(Dx + k * ncisqr, Y + k * ncisqr, ncisqr); /* get D[,,k] */
445                AZERO(Y + k*ncisqr, ncisqr); /* clear Y[,,k] */
446                for (; top < n; top++) {
447                    i = Pattern[top];
448                    Memcpy(Yi, Y + i*ncisqr, ncisqr); /* copy Y[,,i] */
449                    AZERO(Y + i*ncisqr, ncisqr); /* clear Y[,,i] */
450                    p2 = Lp[i] + Lnz[i];
451                    for (p = Lp[i]; p < p2; p++) {
452                        F77_CALL(dgemm)("N", "N", &nci, &nci, &nci, &minus1,
453                                        Lx + p*ncisqr, &nci, Yi, &nci,
454                                        &one, Y + Li[p]*ncisqr, &nci);
455                    }
456                    /* FIXME: Check that this is the correct order and transposition */
457                    F77_CALL(dtrsm)("L", "U", "N", "N", &nci, &nci,
458                                    &one, Dx + i*ncisqr, &nci, Yi, &nci);
459                    F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,
460                                    &one, Dx + k*ncisqr, &nci);
461                    F77_CALL(dtrsm)("L", "U", "T", "N", &nci, &nci,
462                                    &one, Dx + i*ncisqr, &nci, Yi, &nci);
463                    Li[p] = k;      /* store L[k,i] in column form of L */
464                    Memcpy(Lx + p * ncisqr, Yi, ncisqr);
465                    Lnz[i]++;       /* increment count of nonzeros in col i */
466                }
467                F77_CALL(dpotrf)("U", &nci, Dx + k*ncisqr, &nci, &info);
468                if (info) {
469                    Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);
470                    return k;           /* failure, D[,,k] is not positive definite */
471                }
472            }
473            Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);
474            return n;       /* success, diagonal of D is all nonzero */
475      }      }
476      return -1;                  /* keep -Wall happy */      return -1;                  /* keep -Wall happy */
477  }  }
# Line 426  Line 520 
520  }  }
521    
522  /**  /**
523     * Solve a triangular system of the form op(A)*X = alpha*B where A
524     * is a cscBlocked triangular matrix and B is a dense matrix.
525     *
526     * @param uplo 'U' or 'L' for upper or lower
527     * @param trans 'T' or 'N' for transpose or no transpose
528     * @param diag 'U' or 'N' for unit diagonal or non-unit
529     * @param A pointer to a triangular cscBlocked object
530     * @param B pointer to the contents of the matrix B
531     * @param m number of rows in B
532     * @param n number of columns in B
533     * @param ldb leading dimension of B as declared in the calling function
534     */
535    void cscb_trsm(char uplo, char transa, char diag,
536                   double alpha, SEXP A, double B[], int m, int n, int ldb)
537    {
538        int iup = (uplo == 'U' || uplo == 'u'),
539            itr = (transa == 'T' || transa == 't'),
540            iunit = (diag == 'U' || diag == 'u');
541        SEXP ApP = GET_SLOT(A, Matrix_pSym),
542            AxP = GET_SLOT(A, Matrix_xSym);
543        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
544            *Ap = INTEGER(ApP),
545            *xdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
546            i, j, nb = length(ApP) - 1;
547        double *Ax = REAL(GET_SLOT(A, Matrix_xSym));
548    
549        if (xdims[0] != xdims[1])
550            error("Argument A to cscb_trsm is not triangular");
551        if (ldb < m || ldb <= 0 || n <= 0)
552            error("cscb_trsm: inconsistent dims m = %d, n = %d, ldb = %d",
553                  m, n, ldb);
554        if (m != (nb * xdims[0]))
555            error("cscb_trsm: inconsistent dims m = %d, A[%d,%d,]x%d",
556                  m, xdims[0], xdims[1], xdims[2]);
557        if (alpha != 1.0) {
558            for (j = 0; j < n; j++) { /* scale by alpha */
559                for (i = 0; i < m; i++)
560                    B[i + j * ldb] *= alpha;
561            }
562        }
563        if (iunit) {
564            if (xdims[2] < 1) return; /* A is the identity */
565            if (xdims[0] == 1) {
566                if (iup) error("Code for upper triangle not yet written");
567                if (itr) {
568                    for (j = 0; j < n; j++)
569                        R_ldl_ltsolve(m, B + j * ldb, Ap, Ai, Ax);
570                } else {
571                    for (j = 0; j < n; j++)
572                        R_ldl_lsolve(m, B + j * ldb, Ap, Ai, Ax);
573                }
574                return;
575            }
576            error("Code for non-scalar cscBlocked objects not yet written");
577        }
578        error("Code for non-unit cases of cscb_trsm not yet written");
579    }
580    
581    /**
582   * Perform one of the operations B := alpha*op(A)*B or   * Perform one of the operations B := alpha*op(A)*B or
583   * B := alpha*B*op(A) where A and B are both cscBlocked.   * B := alpha*B*op(A) where A and B are both cscBlocked.
584   *   *
# Line 470  Line 623 
623  }  }
624    
625  /**  /**
626     * Expand a column of a compressed, sparse, column-oriented matrix.
627     *
628     * @param dest array to hold the result
629     * @param m number of rows in the matrix
630     * @param j index (0-based) of column to expand
631     * @param Ap array of column pointers
632     * @param Ai array of row indices
633     * @param Ax array of non-zero values
634     *
635     * @return dest
636     */
637    static
638    double *expand_column(double *dest, int m, int j,
639                          const int Ap[], const int Ai[], const double Ax[])
640    {
641        int k, k2 = Ap[j + 1];
642    
643        for (k = 0; k < m; k++) dest[k] = 0.;
644        for (k = Ap[j]; k < k2; k++) dest[Ai[k]] = Ax[k];
645        return dest;
646    }
647    
648    /**
649     * Solve one of the systems op(A)*X = alpha*B or
650     * X*op(A) = alpha*B where A cscBlocked triangular and B is cscBlocked.
651     *
652     * @param side 'L' or 'R' for left or right
653     * @param uplo 'U' or 'L' for upper or lower
654     * @param transa 'T' or 'N' for transpose or no transpose
655     * @param diag 'U' or 'N' for unit diagonal or non-unit
656     * @param alpha scalar multiplier
657     * @param A pointer to a triangular cscBlocked object
658     * @param B pointer to a general cscBlocked matrix
659     */
660    void cscb_trcbsm(char side, char uplo, char transa, char diag,
661                     double alpha, SEXP A, const int Parent[], SEXP B)
662    {
663        int ileft = (side == 'L' || side == 'l'),
664            iup = (uplo == 'U' || uplo == 'u'),
665            itr = (transa == 'T' || transa == 't'),
666            iunit = (diag == 'U' || diag == 'u');
667        SEXP ApP = GET_SLOT(A, Matrix_pSym),
668            AxP = GET_SLOT(A, Matrix_xSym),
669            BpP = GET_SLOT(B, Matrix_pSym),
670            BxP = GET_SLOT(B, Matrix_xSym);
671        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
672            *Ap = INTEGER(ApP),
673            *Bi = INTEGER(GET_SLOT(B, Matrix_iSym)),
674            *Bp = INTEGER(BpP),
675            *axdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
676            *bxdims = INTEGER(getAttrib(BxP, R_DimSymbol)),
677            ncbA = length(ApP) - 1,
678            ncbB = length(BpP) - 1;
679        int i, j, nbx = bxdims[0] * bxdims[1] * bxdims[2];
680        double *Ax = REAL(AxP), *Bx = REAL(BxP);
681    
682        if (axdims[0] != axdims[1])
683            error("Argument A to cscb_trcbm is not triangular");
684        if (alpha != 1.0) {
685            for (i = 0; i < nbx; i++) { /* scale by alpha */
686                REAL(BxP)[i] *= alpha;
687            }
688        }
689        if (iunit && axdims[2] < 1) return; /* A is the identity */
690        if (iunit && axdims[0] == 1) { /* can use R_ldl code */
691            if (!ileft && itr) {    /* case required for lmer */
692                int *BTp, nnz = bxdims[2], nrbB;
693                int *tmp = expand_column_pointers(ncbB, Bp, Calloc(nnz, int));
694                int *BTi = Calloc(nnz, int);
695                double *BTx = Calloc(nnz, double), *rhs;
696    
697                                    /* transpose B */
698                for (i = 0, nrbB = -1; i < nnz; i++) if (Bi[i] > nrbB) nrbB = Bi[i];
699                BTp = Calloc(nrbB, int);
700                triplet_to_col(ncbB, nrbB, nnz, tmp, Bi, Bx, BTp, BTi, BTx);
701                                    /* sanity check */
702                if (BTp[nrbB] != nnz) error("cscb_trcbsm: transpose operation failed");
703                Free(tmp);
704                                    /* Solve one column at a time */
705                rhs = Calloc(ncbB, double);
706                AZERO(Bx, nnz);     /* zero the result */
707                for (i = 0; i < nrbB; i++) {
708                    R_ldl_lsolve(ncbB, expand_column(rhs, ncbB, i, BTp, BTi, BTx),
709                                 Ap, Ai, Ax);
710                    for (j = 0; j < ncbB; j++) { /* write non-zeros in sol'n into B */
711                        if (BTx[j]) Bx[Tind(Bi, Bp, j, i)] = BTx[j];
712                    }
713                    Free(rhs); Free(BTp); Free(BTx); Free(BTi);
714                }
715            }
716            error("cscb_trcbsm: method not yet written");
717        }
718        error("cscb_trcbsm: method not yet written");
719    }
720    
721    /**
722   * Perform one of the matrix-matrix operations   * Perform one of the matrix-matrix operations
723   *      C := alpha*op(A)*op(B) + beta*C   *      C := alpha*op(A)*op(B) + beta*C
724   * on compressed, sparse, blocked matrices.   * on compressed, sparse, blocked matrices.

Legend:
Removed from v.422  
changed lines
  Added in v.434

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