calogrank <- function(ftime, fstatus, grp, cvt, strat=NULL) {
  call <- match.call()
  n0 <- length(ftime)
  cvt <- cbind(cvt)
  if (missing(strat)) strat <- rep(1,n0)
  if (length(fstatus) != n0 | length(grp) != n0 | nrow(cvt) != n0) stop("Data length mismatch")
  p <- ncol(cvt)
# if more than 3 covariates stop
  if (p > 3) stop("if more than 3 covariates use with principal components")

  ii <- !is.na(ftime) & !is.na(fstatus) & !is.na(grp) & !is.na(strat)
  for (i in 1:p) ii <- ii & !is.na(cvt[,i])
  ftime	<- ftime[ii]
  fstatus <- fstatus[ii]
  grp <- grp[ii]
  cvt <- cbind(cvt[ii,])
  strat <- strat[ii]
 
  n0 <- length(ftime)
  ng <- length(unique(grp))
  nstrat <- length(unique(strat))

  bb <- apply(cvt, 2, sd)/n0^0.26  

  rval <- survdiff(Surv(ftime, fstatus) ~ grp + strata(strat))
  rval$call <- call

  ord <- order(strat, ftime, -fstatus)
  osts <- fstatus[ord]
  ogrp <- as.numeric(as.factor(grp))[ord]
  ocov <- as.matrix(cvt[ord,])
  strat <- strat[ord]
  ustat <- rep(0, ng)
  uvar <- matrix(0, ng, ng)
  uexpected <- rep(0, ng*nstrat)
  l <- 0
  for(i in unique(strat)) {
    n <- sum(strat==i)
    zz <- .Fortran("uclrst",
                   as.integer(n),
                   as.integer(ng),
                   as.integer(p),
                   as.double(osts[strat==i]),
                   as.integer(ogrp[strat==i]),
                   as.double(ocov[strat==i,]),
                   a0=double(n),
                   a1=double(n*ng),
                   xi=double(p),
                   xj=double(p),
                   Vii=double(ng),
                   Vij=double(ng),
                   Vji=double(ng),
                   Vjj=double(ng),
                   Vidot=double(ng*n),
                   Vdotj=double(ng*n),
                   Vijij=double(ng*ng),
                   igrp=double(ng),
                   jgrp=double(ng),
                   lrmn=double(ng),
                   lrvar=double(ng*ng),
                   as.double(bb))
    ustat <- ustat + zz$lrmn
    uvar <- uvar + matrix(zz$lrvar,ng,ng)
    uexpected[l*ng + (1:ng)] <- zz$jgrp
    l <- l+1
  }
  chi <- sum(solve(uvar[-1,-1], ustat[-1]) * ustat[-1])
  rval$uexp <- uexpected
  rval$uchisq <- chi
  rval$ustat <- ustat
  rval$uvar <- uvar
  class(rval) <- "calogrank"
  rval
}

print.calogrank <- function(x, digits = max(options()$digits - 4, 3), ...) {

  saveopt <-options(digits=digits)
  on.exit(options(saveopt))

  if (!inherits(x, 'calogrank')) stop("Object is not the result of calogrank")

  if (!is.null(cl<- x$call)) {
    cat("Call:\n")
    dput(cl)
    cat("\n")
  }

  omit <- x$na.action
  if (length(omit)) cat("n=", sum(x$n), ", ", naprint(omit),
					  ".\n\n", sep='')

  if (is.matrix(x$obs)){
    otmp <- apply(x$obs,1,sum)
    etmp <- apply(x$exp,1,sum)
    uetmp <- apply(x$uexp,1,sum)
  }
  else {
    otmp <- x$obs
    etmp <- x$exp
    uetmp <- x$uexp
  }
  df <- (sum(1*(etmp>0))) -1
  temp <- cbind(x$n, otmp, etmp, uetmp, ((otmp-etmp)^2)/ diag(x$var),
					  ((otmp-uetmp)^2)/ diag(x$uvar))
  dimnames(temp) <- list(names(x$n), c("N", "Observed", "Expected",
                                "Exptd-adj","(O-E)^2/V", "(O-aE)^2/aV"))
  print(temp)
  cat("\n Unconditional chisq=", format(round(x$chisq,2)),
      " on", df, "degrees of freedom, p=",
      format(signif(1-pchisq(x$chisq, df),digits)), "\n")
  
  cat("      Adjusted chisq=", format(round(x$uchisq,2)),
      " on", df, "degrees of freedom, p=",
      format(signif(1-pchisq(x$uchisq, df),digits)), "\n")
  
  invisible(x)
}