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 99, Sat Apr 17 16:58:38 2004 UTC revision 100, Sat Apr 17 17:03:59 2004 UTC
# Line 1201  Line 1201 
1201              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);              F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1202              if (info) error("DPOTRF returned error code %d", info);              if (info) error("DPOTRF returned error code %d", info);
1203              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);              F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1204              if (info) error("DPOTRF returned error code %d", info);              if (info) error("DPOTRI returned error code %d", info);
1205          }          }
1206          status[0] = status[1] = 0;          status[0] = status[1] = 0;
1207      }      }
# Line 1209  Line 1209 
1209      return R_NilValue;      return R_NilValue;
1210  }  }
1211    
1212    SEXP ssclme_gradient(SEXP x, SEXP REMLp)
1213    {
1214        SEXP
1215            RZXsl = GET_SLOT(x, Matrix_RZXSym),
1216            ans = PROTECT(duplicate(GET_SLOT(x, Matrix_OmegaSym))),
1217            ncsl = GET_SLOT(x, Matrix_ncSym),
1218            bVar = GET_SLOT(x, Matrix_bVarSym);
1219        int
1220            *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1221            *dims = INTEGER(getAttrib(RZXsl, R_DimSymbol)),
1222            *nc = INTEGER(ncsl),
1223            REML = asLogical(REMLp),
1224            i, info,
1225            n = dims[0],
1226            nf = length(ncsl) - 2,
1227            nobs = nc[nf + 1],
1228            p,
1229            pp1 = dims[1];
1230        double
1231            *RZX = REAL(RZXsl),
1232            *b,
1233            alpha,
1234            one = 1.;
1235    
1236        p = pp1 - 1;
1237        b = RZX + p * n;
1238        ssclme_invert(x);
1239        for (i = 0; i < nf; i++) {
1240            int ki = Gp[i+1] - Gp[i],
1241                nci = nc[i],
1242                mi = ki/nci;
1243            double
1244                *vali = REAL(VECTOR_ELT(ans, i));
1245    
1246            F77_CALL(dpotrf)("U", &nci, vali, &nci, &info);
1247            if (info)
1248                error("DPOTRF returned error code %d for component %d of Omega",
1249                      info, i + 1);
1250            F77_CALL(dpotri)("U", &nci, vali, &nci, &info);
1251            if (info)
1252                error("DPOTRI returned error code %d for component %d of Omega",
1253                      info, i + 1);
1254            alpha = (double) -mi;
1255            F77_CALL(dsyrk)("U", "N", &nci, &ki,
1256                            &one, REAL(VECTOR_ELT(bVar, i)), &nci,
1257                            &alpha, vali, &nci);
1258            alpha = ((double)(REML?(nobs-p):nobs));
1259            F77_CALL(dsyrk)("U", "N", &nci, &mi,
1260                            &alpha, b + Gp[i], &nci,
1261                            &one, vali, &nci);
1262            if (REML) {
1263                int j;
1264                for (j = 0; j < p; j++) {
1265                    F77_CALL(dsyrk)("U", "N", &nci, &mi,
1266                                    &one, RZX + Gp[i] + j*n, &nci,
1267                                    &one, vali, &nci);
1268                }
1269            }
1270        }
1271        UNPROTECT(1);
1272        return ans;
1273    }
1274    
1275  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)  SEXP ssclme_fitted(SEXP x, SEXP facs, SEXP mmats)
1276  {  {
1277      SEXP val, b;      SEXP val, b;
# Line 1286  Line 1349 
1349      UNPROTECT(2);      UNPROTECT(2);
1350      return val;      return val;
1351  }  }
1352    

Legend:
Removed from v.99  
changed lines
  Added in v.100

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