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 108, Mon Apr 19 17:44:22 2004 UTC revision 111, Tue Apr 20 15:46:04 2004 UTC
# Line 432  Line 432 
432          for (k = j+1; k < nf; k++) { /* off-diagonals */          for (k = j+1; k < nf; k++) { /* off-diagonals */
433              int *fpk = INTEGER(VECTOR_ELT(facs, k)),              int *fpk = INTEGER(VECTOR_ELT(facs, k)),
434                  *Apk = Ap + Gp[k],                  *Apk = Ap + Gp[k],
435                  nck = nc[k];                  nck = nc[k],
436                    scalar = ncj == 1 && nck == 1;
437              double              double
438                  *Zk = Z[k];                  *Zk = Z[k], *work;
439                if (!scalar) work = Calloc(ncj * nck, double);
440              for (i = 0; i < nobs; i++) {              for (i = 0; i < nobs; i++) {
441                  int ii, ind = -1, fpji = fpj[i] - 1,                  int ii, ind = -1, fpji = fpj[i] - 1,
442                      row = Gp[j] + fpji * ncj,                      row = Gp[j] + fpji * ncj,
443                      fpki = fpk[i] - 1,                      fpki = fpk[i] - 1,
444                      lastind = Apk[fpki + 1];                      lastind = Apk[fpki*nck + 1];
445                  for (ii = Apk[fpki]; ii < lastind; ii++) {                  for (ii = Apk[fpki*nck]; ii < lastind; ii++) {
446                      if (Ai[ii] == row) {                      if (Ai[ii] == row) {
447                          ind = ii;                          ind = ii;
448                          break;                          break;
449                      }                      }
450                  }                  }
451                  if (ind < 0) error("logic error in ssclme_update_mm");                  if (ind < 0) error("logic error in ssclme_update_mm");
452                  if (Ncj || nck > 1) {                  if (scalar) {   /* update scalars directly */
                                 /* FIXME: run a loop to update */  
                     error("code not yet written");  
                 } else {        /* update scalars directly */  
453                      Ax[ind] += Zj[i] * Zk[i];                      Ax[ind] += Zj[i] * Zk[i];
454                    } else {
455                        int jj, offset = ind - Apk[fpki * nck];
456                        F77_CALL(dgemm)("T", "N", &ncj, &nck, &ione, &one,
457                                        Zj + i, &nobs, Zk + i, &nobs,
458                                        &zero, work, &ncj);
459                        for (jj = 0; jj < nck; jj++) {
460                            ind = Apk[fpki * nck + jj] + offset;
461                            if (Ai[ind] != row)
462                                error("logic error in ssclme_update_mm");
463                            for (ii = 0; ii < ncj; ii++) {
464                                Ax[ind++] += work[jj * ncj + ii];
465                            }
466                  }                  }
467              }              }
468          }          }
469                if (!scalar) Free(work);
470            }
471      }      }
472      Free(Z);      Free(Z);
473      ssclme_transfer_dimnames(x, facs, mmats);      ssclme_transfer_dimnames(x, facs, mmats);
# Line 721  Line 733 
733                          "Rank deficient variance matrix at group %d, level %d",                          "Rank deficient variance matrix at group %d, level %d",
734                          i + 1, j + 1);                          i + 1, j + 1);
735              }              }
736                Free(tmp);
737          }          }
738          return R_NilValue;          return R_NilValue;
739      }      }

Legend:
Removed from v.108  
changed lines
  Added in v.111

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