Check for DimNames propagation in coercion and other operations.
------
- Report the problem in the Linux ldexp manual page.  The second and
  third calls in the Synopsis should be to ldexpf and ldexpl.

- [,] indexing: for sparse "works", but not yet for negative indices!

- consider moving alloc3Darray from ./src/Mutils.c to
  $(RSRC)/src/base/array.c
------
- provide methods for "dspMatrix" and "dppMatrix"!

- implement (more) methods for supporting "packed" (symmetric / triangular)
  matrices; particularly something like pack() and unpack()  [to/from our
  classes from/to "numeric"] --- have already man/unpack.Rd but no method yet!

  (have some dtr* <-> dtp*)

-----

- combine the C functions for multiplication by special forms and
  solution wrt special forms by using a 'right' argument and a
  'classed' argument.
   [done with dgeMatrix_matrix_mm();  not yet for other classes;
    and for _crossprod()]

- add more comprehensive examples / tests for Schur decomposition

- arithmetic for sparse matrices:
	     <sparseMatrix>  o  <same-dim-sparseMatrix>
  should return a sparse matrix  for at least "+" and "*" , also %%,
  and "/" and "%/%" at least when the RHS is non-zero a scalar.
  Challenge: nice implementation (``common non-0''; but Tsparse* is not uniq).

-----

- "Math2" , "Math", "Arith":
   keep triangular and symmetric Matrices when appropriate:
   particularly desirable for  "Math2": round(), signif()

- "d" <-> "l" coercion for all "[TCR]" sparse matrices is really trivial:
  "d" -> "l" : drops the 'x' slot
  "l" -> "d" : construct an 'x' slot of all '1'
  We currently have many of these conversions explicitly, e.g.
   setAs("dsTMatrix", "lsTMatrix",
      function(from) new("lsTMatrix", i = from@i, j = from@j, uplo = from@uplo,
                         Dim = from@Dim, Dimnames = from@Dimnames))
  but I would rather want to automatically construct all these coercion
  methods at once by a ``method constructor'', i.e.,
  for all  "dsparse*" -> "lsparse*" and vice versa.
  How can one do this {in a documented way} ?

- Think of constructing  setAs(...) calls automatically in order to
  basically enable all ``sensible'' as(fromMatrix, toMatrix)  calls,
  possibly using canCoerce(.)

- setAs(<Mcl>,  "[dln]Matrix") for <Mcl> in {Matrix or denseMatrix + sparseMatrix}

- When we have a packed matrix, it's a waste to go through "full" to "sparse":
  ==> implement
	setAs("dspMatrix", "sparseMatrix")
	setAs("dppMatrix", "sparseMatrix")
	setAs("dtpMatrix", "sparseMatrix")
  and the same for "lsp" , "ltp"  and  "nsp" , "ntp" !

- tcrossprod(x, y) : do provide methods for y != NULL
  calling Lapack's DGEMM for "dense"
  [2005-12-xx: done for dgeMatrix at least]

- BUGlet:  Shouldn't lose factorization here:
  h6 <- Hilbert(6); chol(h6) ; str(h6) # has factor
  str(H6 <- as(h6, "dspMatrix"))       # has lost factor
  ## and the same in a similar situation involving  "dpo", "dpp"

- Things like  M[upper.tri(M)] are not really most useful for  sparse
  matrices.  --> provide generic functions
  upperTriMatrix(), lowerTriMatrix()  both with argument  'diag = TRUE'
  (which can be set to FALSE of course) which are used to extract a
  triangle from an arbitrary sparse matrix and  return a  "dtCMatrix".

- Factorizations: LU done; also Schur()  for  *sparse*  Matrices.

- band(), triu(), tril() for *all* including "matrix", not just sparse matrices

- is.na() method for all our matrices [ ==> which(*, arr.ind=TRUE) might work ]

- use  .Call(Csparse_drop, M, tol) in more places,
  both with 'tol = 0.' to drop "values that happen to be 0" and for
  zapsmall() methods for Csparse*

- implement .Call(Csparse_scale, ....) interfacing to cholmod_scale()
  in src/CHOLMOD/Include/cholmod_matrixops.h : for another function
  specifically for multiplying a cholmod_sparse object by a diagonal matrix.
  Use it in %*% and [t]crossprod methods.

- chol() and determinant() should ``work'': proper result or "good" error
  message.

- make sure *all* group methods have (maybe "bail-out") setMethod for "Matrix".
  e.g. zapsmall(<pMatrix>) fails "badly"

- speedup: pass class definition to non0ind() [check all calls ..]

- sum(): implement methods which work for *all* our matrices.

- Implement  expand(.) for the Cholesky() results
  "dCHMsimpl" and  "dCHMsuper"  -- currently have no *decent* way to get at
  the matrix factors of the corresponding matrix factorization !!

- rbind(<sparse>, <dense>) does not work  (e.g. <dgC>, <dge>)

- make sure  M[FALSE, FALSE]  works for all Matrices
  {e.g. fails for M <- Diagonal(4)}

- <sparse> %*% <dense>  {also in crossprod/tcrossprod}  currently always
  returns <dense>, since --> Csparse_dense_prod --> cholmod_sdmult
  and that does only return dense.
  When the sparse matrix is very sparse, i.e. has many rows with only zero
  entries, it would make much sense to return sparse.

- sparse-symmetric + diagonal should stay sparse-symmetric
  (only stays sparse): Matrix(0, 4, 4) + Diagonal(4, 1:4)
  --> R/diagMatrix.R ('FIXME')
  but also R/Ops.R  to ensure  sp-sym. + sp-sym. |-> sp-sym.  etc

- Diagonal(n) %*% A ---  too slow!! --> ~/R/MM/Pkg-ex/Matrix/diag-Tamas-ex.R

- For a square sparse matrix 'b' {typically dgCMatrix or dgTMatrix},
  we'd want a function  "Mat_plus_t_Mat" <- function(b) {....}
  which computes the symmetric sparse matrix   b + t(b)
  in way that never works with size-doubled vectors from  b@i etc..

- ! <symmetricMatrix>  loses symmetry, both for dense and sparse matrices.
  !M  where M is "sparseMatrix", currently always gives dense. This only
  makes sense when M is ``really sparse''.

- msy <- as(matrix(c(2:1,1:2),2), "dsyMatrix"); str(msy)

  shows that the Cholesky factorization is computed ``too quickly''.
  Can be a big pain for largish matrices, when it is unneeded.

- example(Cholesky, echo=FALSE) ; cm <- chol(mtm); str(cm); str(mtm)

  shows that chol() does not seems to make use of an already
  present factorization and rather uses one with more '0' in x slot.

- diag(m) <- val    currently automatically works via  m[cbind(i,i)] <- val
  However,
  we need methods for 'diag<-' at least for diagonalMatrix,
  triangularMatrix, and probably also "dense*general*Matrix" since the
  above currently goes via "matrix" and back instead of using the 'x' slot
  directly.

- image(M, ..): Think about an optional smart option which keeps
   "0 |-> transparent" and allows colors to differentiate negative and
   positive entries.

- examples for solve( Cholesky(.), b, system = c("A", "LDLt"....))
  probably rather in man/CHMfactor-class.Rd than man/Cholesky.Rd

- (A + tr(A))/2  := the symmetric part of A, is needed in several
  circumstances; unfortunately it's not "smart" (preserving symmetry, ...)
  --> define a generic and methods for it!
  Names:  symPart(A) or  symMat(A) or symmetrize(A) or ... ?
  Googling around I found that Nick Higham has a GPL contributed Matlab
  toolbox where he uses  symmpart(A) := (A + A') /. 2
  {and  skewpart(A) := (A - A') /. 2}

- tr(A %*% B) {and even  tr(A %*% B %*% C) ...} are also needed
  frequently in some computations {conditional normal distr. ...}.
  Since this can be done faster than by
    sum(diag(A %*% B))  even for traditional matrices, e.g.
    	       sum(A * t(B)) or {even faster for "full" mat}
	       crossprod(as.vector(A), as.vector(B))
  and even more so for, e.g.  <sparse> %*% <dense>
  {used in Soeren's 'gR' computations},
  we should also provide a generic and methods.

- qr.R(qr(x)) may differ for the "same" matrix, depending on it being
  sparse or dense:
    "qr.R(<sparse>) may differ from qr.R(<dense>) because of permutations"

  This is not really acceptable and currently influences  rcond() as well.