#### Methods for the virtual class 'CsparseMatrix' of sparse matrices stored in
####  "column compressed" format.
#### -- many more specific things are e.g. in ./dgCMatrix.R

setAs("CsparseMatrix", "TsparseMatrix",
      function(from)
          ## |-> cholmod_C -> cholmod_T -> chm_triplet_to_SEXP
          ## modified to support triangular (../src/Csparse.c)
          .Call(Csparse_to_Tsparse, from, is(from, "triangularMatrix")))


## special cases (when a specific "to" class is specified)
setAs("dgCMatrix", "dgTMatrix",
      function(from) .Call(Csparse_to_Tsparse, from, FALSE))

setAs("dsCMatrix", "dsTMatrix",
      function(from) .Call(Csparse_to_Tsparse, from, FALSE))

setAs("dsCMatrix", "dgCMatrix",
      function(from) .Call(Csparse_symmetric_to_general, from))

for(prefix in c("d", "l", "n"))
    setAs(paste(prefix,"sCMatrix",sep=''), "generalMatrix",
	  function(from) .Call(Csparse_symmetric_to_general, from))

setAs("dtCMatrix", "dtTMatrix",
      function(from) .Call(Csparse_to_Tsparse, from, TRUE))

setAs("CsparseMatrix", "denseMatrix",
      function(from) {
	  ## |-> cholmod_C -> cholmod_dense -> chm_dense_to_dense
	  cld <- getClassDef(class(from))
	  if (extends(cld, "generalMatrix"))
	      .Call(Csparse_to_dense, from)
	  else {
	      ## Csparse_to_dense  loses symmetry and triangularity properties.
	      ## With suitable changes to chm_dense_to_SEXP (../src/chm_common.c)
	      ## we could do this in C code -- or do differently in C {FIXME!}
	      if (extends(cld, "triangularMatrix") && from@diag == "U")
		  from <- .Call(Csparse_diagU2N, from)
	      as(.Call(Csparse_to_dense, from), # -> "[dln]geMatrix"
		 paste(.M.kind(from, cld),
		       .dense.prefixes[.M.shape(from, cld)], "Matrix", sep=''))
	  }
      })

## special cases (when a specific "to" class is specified)
setAs("dgCMatrix", "dgeMatrix",
      function(from) .Call(Csparse_to_dense, from))

## cholmod_sparse_to_dense converts symmetric storage to general
## storage so symmetric classes are ok for conversion to matrix.
## unit triangular needs special handling
setAs("CsparseMatrix", "matrix",
      function(from) {
          ## |-> cholmod_C -> cholmod_dense -> chm_dense_to_matrix
          if (is(from, "triangularMatrix") && from@diag == "U")
              from <- .Call(Csparse_diagU2N, from)
          .Call(Csparse_to_matrix, from)
      })

setAs("CsparseMatrix", "symmetricMatrix",
      function(from) {
	  if(isSymmetric(from)) {
	      isTri <- is(from, "triangularMatrix")# i.e. effectively *diagonal*
	      if (isTri && from@diag == "U")
		  from <- .Call(Csparse_diagU2N, from)
	      .Call(Csparse_general_to_symmetric, from,
		    uplo = if(isTri) from@uplo else "U")
	  } else
	  stop("not a symmetric matrix; consider forceSymmetric() or symmpart()")
      })


.validateCsparse <- function(x, sort.if.needed = FALSE)
    .Call(Csparse_validate2, x, sort.if.needed)
##-> to be used in sparseMatrix(.), e.g. --- but is unused currently
## NB: 'sort.if.needed' is called 'maybe_modify' in C -- so be careful

### Some group methods:

setMethod("Math",
	  signature(x = "CsparseMatrix"),
	  function(x) {
	      f0 <- callGeneric(0.)
	      if(is0(f0)) {
		  ## sparseness, symm., triang.,... preserved
                  cl <- class(x)
                  has.x <- !extends(cl, "nsparseMatrix")
                  ## has.x  <==> *not* nonzero-pattern == "nMatrix"
                  if(has.x) {
                      type <- storage.mode(x@x)
                      r <- callGeneric(x@x)
                  } else { ## nsparseMatrix
                      type <- ""
		      r <- rep.int(as.double(callGeneric(TRUE)),
				   switch(.sp.class(cl),
					  CsparseMatrix = length(x@i),
					  TsparseMatrix = length(x@i),
					  RsparseMatrix = length(x@j)))
                  }
		  if(type == storage.mode(r)) {
		      x@x <- r
		      x
		  } else { ## e.g. abs( <lgC> ) --> integer Csparse
		      ## FIXME: when we have 'i*' classes, use them here:
		      rx <- new(sub("^.", "d", cl))
		      rx@x <- as.double(r)
		      ## result is "same"
		      sNams <- slotNames(cl)
		      for(nm in sNams[sNams != "x"])
			  slot(rx, nm) <- slot(x, nm)
		      rx
		  }
	      } else { ## no sparseness
		  callGeneric(as_dense(x))
	      }
	  }) ## {Math}


### Subsetting -- basic things (drop = "missing") are done in ./Matrix.R
### ---------- "["  and (currently) also ./sparseMatrix.R

subCsp_cols <- function(x, j, drop)
{
    ## x[ , j, drop=drop]   where we know that	x is Csparse*
    dn <- x@Dimnames
    jj <- intI(j, n = x@Dim[2], dn[[2]], give.dn = FALSE)
    r <- .Call(Csparse_submatrix, x, NULL, jj)
    if(!is.null(n <- dn[[1]])) r@Dimnames[[1]] <- n
    if(!is.null(n <- dn[[2]])) r@Dimnames[[2]] <- n[jj+1L]
    if(drop && any(r@Dim == 1L)) drop(as(r, "matrix")) else r
}

subCsp_rows <- function(x, i, drop)# , cl = getClassDef(class(x))
{
    ## x[ i,  drop=drop]   where we know that  x is Csparse*
    dn <- x@Dimnames
    ii <- intI(i, n = x@Dim[1], dn[[1]], give.dn = FALSE)
    r <- .Call(Csparse_submatrix, x, ii, NULL)
    if(!is.null(n <- dn[[1]])) r@Dimnames[[1]] <- n[ii+1L]
    if(!is.null(n <- dn[[2]])) r@Dimnames[[2]] <- n
    if(drop && any(r@Dim == 1L)) drop(as(r, "matrix")) else r
}

subCsp_ij <- function(x, i, j, drop)
{
    ## x[i, j, drop=drop]   where we know that	x is Csparse*
    d <- x@Dim
    dn <- x@Dimnames
    ## Take care that	x[i,i]	for symmetricM* stays symmetric
    i.eq.j <- identical(i,j) # < want fast check
    ii <- intI(i, n = d[1], dn[[1]], give.dn = FALSE)
    jj <- if(i.eq.j && d[1] == d[2]) ii
    else intI(j, n = d[2], dn[[2]], give.dn = FALSE)
    r <- .Call(Csparse_submatrix, x, ii, jj)
    if(!is.null(n <- dn[[1]])) r@Dimnames[[1]] <- n[ii + 1L]
    if(!is.null(n <- dn[[2]])) r@Dimnames[[2]] <- n[jj + 1L]
    if(!i.eq.j) {
	if(drop && any(r@Dim == 1L)) drop(as(r, "matrix")) else r
    } else { ## i == j
	if(drop) drop <- any(r@Dim == 1L)
	if(drop)
	    drop(as(r, "matrix"))
	else if(extends((cx <- getClassDef(class(x))),
                        "symmetricMatrix")) ## TODO? make more efficient:
	    .gC2sym(r, uplo = x@uplo)## preserve uplo !
	else if(extends(cx, "triangularMatrix") && !is.unsorted(ii))
	    as(r, "triangularMatrix")
	else r
    }
}

setMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
			 drop = "logical"),
	  function (x, i,j, ..., drop) {
	      na <- nargs()
	      Matrix.msg("Csp[i,m,l] : nargs()=",na, .M.level = 2)
	      if(na == 4)
		  subCsp_rows(x, i, drop=drop)
	      else if(na == 3)
		  .M.vectorSub(x, i) # as(x, "TsparseMatrix")[i, drop=drop]
	      else ## should not happen
		  stop("Matrix-internal error in <CsparseM>[i,,d]; please report")
	  })

setMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
			 drop = "logical"),
	  function (x,i,j, ..., drop) {
	      Matrix.msg("Csp[m,i,l] : nargs()=",nargs(), .M.level = 2)
	      subCsp_cols(x, j, drop=drop)
	  })

setMethod("[", signature(x = "CsparseMatrix",
			 i = "index", j = "index", drop = "logical"),
	  function (x, i, j, ..., drop) {
	      Matrix.msg("Csp[i,i,l] : nargs()=",nargs(), .M.level = 2)
	      subCsp_ij(x, i, j, drop=drop)
	  })




## workhorse for "[<-" -- both for d* and l*  C-sparse matrices :
## ---------     -----  FIXME(2): keep in sync with replTmat() in ./Tsparse.R
replCmat <- function (x, i, j, ..., value)
{
    di <- dim(x)
    dn <- dimnames(x)
    iMi <- missing(i)
    jMi <- missing(j)
    spV <- is(value, "sparseVector")
    na <- nargs()
    Matrix.msg("replCmat[x,i,j,.., val] : nargs()=", na,"; ",
	       if(iMi | jMi) sprintf("missing (i,j) = (%d,%d)", iMi,jMi),
	       .M.level = 2)
    if(na == 3) { ## "vector (or 2-col) indexing"  M[i] <- v
	x <- as(x, "TsparseMatrix")
	x[i] <- value # may change class e.g. from dtT* to dgT*
	clx <- sub(".Matrix$", "CMatrix", class(x))
	return(if(any(is0(x@x))) ## drop all values that "happen to be 0"
	       drop0(x, is.Csparse=FALSE) else as_CspClass(x, clx))
    }
    ## nargs() == 4 :

    lenV <- length(value)
    i1 <- if(iMi) 0:(di[1] - 1L) else .ind.prep2(i, 1, di, dn)
    i2 <- if(jMi) 0:(di[2] - 1L) else .ind.prep2(j, 2, di, dn)
    dind <- c(length(i1), length(i2)) # dimension of replacement region
    lenRepl <- prod(dind)
    if(lenV == 0) {
        if(lenRepl != 0)
            stop("nothing to replace with")
        else return(x)
    }
    ## else: lenV := length(value)	 is > 0
    if(lenRepl %% lenV != 0)
	stop("number of items to replace is not a multiple of replacement length")
    if(lenV > lenRepl)
	stop("too many replacement values")

    clx <- class(x)
    clDx <- getClassDef(clx) # extends() , is() etc all use the class definition

    ## keep "symmetry" if changed here:
    x.sym <- extends(clDx, "symmetricMatrix")
    if(x.sym) { ## only half the indices are there..
	## using array() for large dind is a disaster...
	mkArray <- if(spV) # TODO: room for improvement
	    function(v, dim) spV2M(v, dim[1],dim[2]) else array
	x.sym <-
	    (dind[1] == dind[2] && all(i1 == i2) &&
	     (lenRepl == 1 || lenV == 1 ||
	      isSymmetric(mkArray(value, dim=dind))))
	## x.sym : result is *still* symmetric
	x <- .Call(Csparse_symmetric_to_general, x) ## but do *not* redefine clx!
    }
    else if(extends(clDx, "triangularMatrix")) {
        xU <- x@uplo == "U"
	r.tri <- ((any(dind == 1) || dind[1] == dind[2]) &&
		  if(xU) max(i1) <= min(i2) else max(i2) <= min(i1))
	if(r.tri) { ## result is *still* triangular
            if(any(i1 == i2)) # diagonal will be changed
                x <- diagU2N(x) # keeps class (!)
	}
	else { # go to "generalMatrix" and continue
	    x <- as(x, paste(.M.kind(x), "gCMatrix", sep='')) ## & do not redefine clx!
	}
    }

    if(extends(clDx, "dMatrix")) {
	has.x <- TRUE
	x <- .Call(Csparse_subassign,
                   if(clx == "dgCMatrix")x else as(x, "dgCMatrix"),
                   i1, i2,
                   as(value, "dsparseVector"))
    }
    else {
    xj <- .Call(Matrix_expand_pointers, x@p)
    sel <- (!is.na(match(x@i, i1)) &
	    !is.na(match( xj, i2)))
    has.x <- "x" %in% slotNames(clDx)# === slotNames(x),
    ## has.x  <==> *not* nonzero-pattern == "nMatrix"

    if(has.x && sum(sel) == lenRepl) { ## all entries to be replaced are non-zero:
        ## need indices instead of just 'sel', for, e.g.,  A[2:1, 2:1] <- v
	non0 <- cbind(match(x@i[sel], i1),
		      match(xj [sel], i2)) - 1L
	iN0 <- 1L + .Call(m_encodeInd, non0, di = dind, checkBounds = FALSE)

        has0 <-
            if(spV) length(value@i) < lenV else any(value[!is.na(value)] == 0)
        if(lenV < lenRepl)
            value <- rep(value, length = lenRepl)
	## Ideally we only replace them where value != 0 and drop the value==0
	## ones; FIXME: see Davis(2006) "2.7 Removing entries", p.16, e.g. use cs_dropzeros()
        ##       but really could be faster and write something like cs_drop_k(A, k)
	## v0 <- 0 == value
	## if (lenRepl == 1) and v0 is TRUE, the following is not doing anything
	##-  --> ./dgTMatrix.R	and its	 replTmat()
	## x@x[sel[!v0]] <- value[!v0]
        x@x[sel] <- as.vector(value[iN0])
	if(extends(clDx, "compMatrix") && length(x@factors)) # drop cashed ones
	    x@factors <- list()
        if(has0) x <- .Call(Csparse_drop, x, 0)

	return(if(x.sym) as_CspClass(x, clx) else x)
    }
    ## else go via Tsparse.. {FIXME: a waste! - we already have 'xj' ..}
    ## and inside  Tsparse... the above i1, i2,..., sel  are *all* redone!
## Happens too often:
##     Matrix.msg("wasteful C -> T -> C in replCmat(x,i,j,v) for <sparse>[i,j] <- v")
    x <- as(x, "TsparseMatrix")
    if(missing(i))
	x[ ,j] <- value
    else if(missing(j))
	x[i, ] <- value
    else
	x[i,j] <- value
    if(extends(clDx, "compMatrix") && length(x@factors)) # drop cashed ones
	x@factors <- list()
}# else{ not using new memory-sparse  code
    if(has.x && any(is0(x@x))) ## drop all values that "happen to be 0"
	as_CspClass(drop0(x), clx)
    else as_CspClass(x, clx)
} ## replCmat

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
                                value = "replValue"),
                 replCmat)

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
                                value = "replValue"),
                 replCmat)

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "index",
				value = "replValue"),
                 replCmat)

### When the RHS 'value' is  a sparseVector, now can use  replCmat  as well
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "missing", j = "index",
				value = "sparseVector"),
		 replCmat)

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "missing",
				value = "sparseVector"),
		 replCmat)

setReplaceMethod("[", signature(x = "CsparseMatrix", i = "index", j = "index",
				value = "sparseVector"),
		 replCmat)

## A[ ij ] <- value,  where ij is (i,j) 2-column matrix
setReplaceMethod("[", signature(x = "CsparseMatrix", i = "matrix", j = "missing",
				value = "replValue"),
		 function(x, i, value)
		 ## goto Tsparse modify and convert back:
		 as(.TM.repl.i.mat(as(x, "TsparseMatrix"), i=i, value=value),
		    "CsparseMatrix"))
## more in ./Matrix.R


setMethod("t", signature(x = "CsparseMatrix"),
	  function(x) .Call(Csparse_transpose, x, is(x, "triangularMatrix")))


## NB: have extra tril(), triu() methods for symmetric ["dsC" and "lsC"] and
## NB: for all triangular ones, where the latter may 'callNextMethod()' these:
setMethod("tril", "CsparseMatrix",
	  function(x, k = 0, ...) {
	      k <- as.integer(k[1])
	      dd <- dim(x); sqr <- dd[1] == dd[2]
	      stopifnot(-dd[1] <= k, k <= dd[1]) # had k <= 0
	      r <- .Call(Csparse_band, x, -dd[1], k)
	      ## return "lower triangular" if k <= 0
	      if(sqr && k <= 0)
		  as(r, paste(.M.kind(x), "tCMatrix", sep='')) else r
	  })

setMethod("triu", "CsparseMatrix",
	  function(x, k = 0, ...) {
	      k <- as.integer(k[1])
	      dd <- dim(x); sqr <- dd[1] == dd[2]
	      stopifnot(-dd[1] <= k, k <= dd[1]) # had k >= 0
	      r <- .Call(Csparse_band, x, k, dd[2])
	      ## return "upper triangular" if k >= 0
	      if(sqr && k >= 0)
		  as(r, paste(.M.kind(x), "tCMatrix", sep='')) else r
	  })

setMethod("band", "CsparseMatrix",
	  function(x, k1, k2, ...) {
	      k1 <- as.integer(k1[1])
	      k2 <- as.integer(k2[1])
	      dd <- dim(x); sqr <- dd[1] == dd[2]
	      stopifnot(-dd[1] <= k1, k1 <= k2, k2 <= dd[2])
	      r <- .Call(Csparse_band, diagU2N(x), k1, k2)
	      if(sqr && k1 * k2 >= 0) ## triangular
		  as(r, paste(.M.kind(x), "tCMatrix", sep=''))
	      else if (k1 < 0  &&  k1 == -k2  && isSymmetric(x)) ## symmetric
		  as(r, paste(.M.kind(x), "sCMatrix", sep=''))
	      else
		  r
	  })

setMethod("diag", "CsparseMatrix",
	  function(x, nrow, ncol) {
              ## "FIXME": could be more efficient; creates new ..CMatrix:
	      dm <- .Call(Csparse_band, diagU2N(x), 0, 0)
	      dlen <- min(dm@Dim)
	      ind1 <- dm@i + 1L	# 1-based index vector
	      if (is(dm, "nMatrix")) {
		  val <- rep.int(FALSE, dlen)
		  val[ind1] <- TRUE
	      }
	      else if (is(dm, "lMatrix")) {
		  val <- rep.int(FALSE, dlen)
		  val[ind1] <- as.logical(dm@x)
	      }
	      else {
		  val <- rep.int(0, dlen)
		  ## cMatrix not yet active but for future expansion
		  if (is(dm, "cMatrix")) val <- as.complex(val)
		  val[ind1] <- dm@x
	      }
	      val
	  })

setMethod("writeMM", "CsparseMatrix",
	  function(obj, file, ...)
          .Call(Csparse_MatrixMarket, obj, as.character(file)))

setMethod("Cholesky", signature(A = "CsparseMatrix"),
	  function(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...)
	  Cholesky(as(A, "symmetricMatrix"),
		   perm=perm, LDL=LDL, super=super, Imult=Imult, ...))

## TODO (in ../TODO for quite a while .....):
setMethod("Cholesky", signature(A = "nsparseMatrix"),
	  function(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...)
	  stop("Cholesky(<nsparse...>) -> *symbolic* factorization -- not yet implemented"))