SCM

SCM Repository

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

Diff of /pkg/src/ssclme.c

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

revision 18, Tue Mar 23 23:29:32 2004 UTC revision 19, Tue Mar 23 23:29:55 2004 UTC
# Line 957  Line 957 
957  }  }
958    
959  /**  /**
960     * Extract the upper triangles of the Omega matrices in the unconstrained
961     * parameterization.
962     * (These are not in any sense "coefficients" but the extractor is
963     * called coef for historical reasons.)
964     *
965     * @param x pointer to an ssclme object
966     *
967     * @return numeric vector of the values in the upper triangles of the
968     * Omega matrices
969     */
970    SEXP ssclme_coefUnc(SEXP x)
971    {
972        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
973        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
974            i, nf = length(Omega), vind;
975        SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
976        double *vv = REAL(val);
977    
978        vind = 0;
979        for (i = 0; i < nf; i++) {
980            int  nci = nc[i];
981            if (nci == 1) {
982                vv[vind++] = log(REAL(VECTOR_ELT(Omega, i))[0]);
983            } else {
984                int j, k, ncip1 = nci + 1, ncisq = nci * nci;
985                double *tmp = Memcpy(Calloc(ncisq, double),
986                                     REAL(VECTOR_ELT(Omega, i)), ncisq);
987                F77_CALL(dpotrf)("U", &nci, tmp, &nci, &j);
988                if (j)              /* should never happen */
989                    error("DPOTRF returned error code %d", j);
990                for (j = 0; j < nci; j++) {
991                    double diagj = tmp[j * ncip1];
992                    vv[vind++] = 2. * log(diagj);
993                    for (k = j + 1; k < nci; k++) {
994                        tmp[j + k * nci] /= diagj;
995                    }
996                }
997                for (j = 0; j < nci; j++) {
998                    for (k = j + 1; k < nci; k++) {
999                        vv[vind++] = tmp[j + k * nci];
1000                    }
1001                }
1002                Free(tmp);
1003            }
1004        }
1005        UNPROTECT(1);
1006        return val;
1007    }
1008    
1009    /**
1010     * Assign the upper triangles of the Omega matrices in the unconstrained
1011     * parameterization.
1012     *
1013     * @param x pointer to an ssclme object
1014     * @param coef pointer to an numeric vector of appropriate length
1015     *
1016     * @return R_NilValue
1017     */
1018    SEXP ssclme_coefGetsUnc(SEXP x, SEXP coef)
1019    {
1020        SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1021        int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1022            cind, i, nf = length(Omega),
1023            *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1024        double *cc = REAL(coef);
1025    
1026        if (length(coef) != coef_length(nf, nc) || !isReal(coef))
1027            error("coef must be a numeric vector of length %d",
1028                  coef_length(nf, nc));
1029        cind = 0;
1030        for (i = 0; i < nf; i++) {
1031            int nci = nc[i];
1032            if (nci == 1) {
1033                REAL(VECTOR_ELT(Omega, i))[0] = exp(cc[cind++]);
1034            } else {
1035                int odind = cind + nci, /* off-diagonal index */
1036                    j, k,
1037                    ncip1 = nci + 1,
1038                    ncisq = nci * nci;
1039                double
1040                    *omgi = REAL(VECTOR_ELT(Omega, i)),
1041                    *tmp = Calloc(ncisq, double),
1042                    diagj, one = 1.;
1043                                    /* LD in omgi and L' in tmp */
1044                memset(omgi, 0, sizeof(double) * ncisq);
1045                for (j = 0; j < nci; j++) {
1046                    omgi[j * ncip1] = diagj = exp(cc[cind++]);
1047                    for (k = j + 1; k < nci; k++) {
1048                        omgi[j*nci + k] = diagj * (tmp[k*nci + j] = cc[odind++]);
1049                    }
1050                }
1051                F77_CALL(dtrmm)("R", "U", "N", "U", &nci, &nci, &one,
1052                                tmp, &nci, omgi, &nci);
1053                Free(tmp);
1054                cind = odind;
1055            }
1056        }
1057        status[0] = status[1] = 0;
1058        return x;
1059    }
1060    
1061    /**
1062   * Assign the upper triangles of the Omega matrices.   * Assign the upper triangles of the Omega matrices.
1063   * (These are not in any sense "coefficients" but are   * (These are not in any sense "coefficients" but are
1064   * called coef for historical reasons.)   * called coef for historical reasons.)

Legend:
Removed from v.18  
changed lines
  Added in v.19

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