text: /* ========================================================================== */
  /* === Core/cholmod_sparse ================================================== */
  /* ========================================================================== */
  
  /* -----------------------------------------------------------------------------
   * CHOLMOD/Core Module.  Version 1.0.  Copyright (C) 2005, Univ. of Florida.
   * Author: Timothy A. Davis
   * The CHOLMOD/Core Module is licensed under Version 2.1 of the GNU
   * Lesser General Public License.  See lesser.txt for a text of the license.
   * CHOLMOD is also available under other licenses; contact authors for details.
   * http://www.cise.ufl.edu/research/sparse
   * -------------------------------------------------------------------------- */
  
  /* Core utility routines for the cholmod_sparse object:
   *
   * A sparse matrix is held in compressed column form.  In the basic type
   * ("packed", which corresponds to a MATLAB sparse matrix), an n-by-n matrix
   * with nz entries is held in three arrays: p of size n+1, i of size nz, and x
   * of size nz.  Row indices of column j are held in i [p [j] ... p [j+1]-1] and
   * in the same locations in x.  There may be no duplicate entries in a column.
   * Row indices in each column may be sorted or unsorted (CHOLMOD keeps track).
   *
   * Primary routines:
   * -----------------
   * cholmod_allocate_sparse	allocate a sparse matrix
   * cholmod_free_sparse		free a sparse matrix
   *
   * Secondary routines:
   * -------------------
   * cholmod_reallocate_sparse	change the size (# entries) of sparse matrix
   * cholmod_nnz			number of nonzeros in a sparse matrix
   * cholmod_speye		sparse identity matrix
   * cholmod_spzeros		sparse zero matrix
   * cholmod_copy_sparse		create a copy of a sparse matrix
   *
   * All xtypes are supported (pattern, real, complex, and zomplex)
   */
  
  #include "cholmod_core.h"
  #include "cholmod_internal.h"
  
  
  /* ========================================================================== */
  /* === cholmod_allocate_sparse ============================================== */
  /* ========================================================================== */
  
  /* Allocate space for a matrix.  A->i and A->x are not initialized.  A->p
   * (and A->nz if A is not packed) are set to zero, so a matrix containing no
   * entries (all zero) is returned.  See also cholmod_spzeros.
   *
   * workspace: none
   */
  
  cholmod_sparse *CHOLMOD(allocate_sparse)
  (
      /* ---- input ---- */
      size_t nrow,	/* # of rows of A */
      size_t ncol,	/* # of columns of A */
      size_t nzmax,	/* max # of nonzeros of A */
      int sorted,		/* TRUE if columns of A sorted, FALSE otherwise */
      int packed,		/* TRUE if A will be packed, FALSE otherwise */
      int stype,		/* stype of A */
      int xtype,		/* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
      /* --------------- */
      cholmod_common *Common
  )
  {
      cholmod_sparse *A ;
      Int *Ap, *Anz ;
      size_t nzmax0 ;
      Int j ;
  
      /* ---------------------------------------------------------------------- */
      /* get inputs */
      /* ---------------------------------------------------------------------- */
  
      RETURN_IF_NULL_COMMON (NULL) ;
      if (stype != 0 && nrow != ncol)
      {
  	ERROR (CHOLMOD_INVALID, "rectangular matrix with stype != 0 invalid") ;
  	return (NULL) ;
      }
      if (xtype < CHOLMOD_PATTERN || xtype > CHOLMOD_ZOMPLEX)
      {
  	ERROR (CHOLMOD_INVALID, "xtype invalid") ;
  	return (NULL) ;
      }
      Common->status = CHOLMOD_OK ;
  
      /* ---------------------------------------------------------------------- */
      /* allocate header */
      /* ---------------------------------------------------------------------- */
  
      A = CHOLMOD(malloc) (1, sizeof (cholmod_sparse), Common) ;
      if (Common->status < CHOLMOD_OK)
      {
  	return (NULL) ;	    /* out of memory */
      }
      PRINT1 (("cholmod_allocate_sparse %d-by-%d nzmax %d sorted %d packed %d"
  		" xtype %d\n", nrow, ncol, nzmax, sorted, packed, xtype)) ;
  
      nzmax = MAX (1, nzmax) ;
  
      A->nrow = nrow ;
      A->ncol = ncol ;
      A->nzmax = nzmax ;
      A->packed = packed ;    /* default is packed (A->nz not present) */
      A->stype = stype ;
      A->itype = ITYPE ;
      A->xtype = xtype ;
      A->dtype = DTYPE ;
  
      A->nz = NULL ;
      A->p = NULL ;
      A->i = NULL ;
      A->x = NULL ;
      A->z = NULL ;
  
      /* A 1-by-m matrix always has sorted columns */
      A->sorted = (nrow <= 1) ? TRUE : sorted ;
  
      /* ---------------------------------------------------------------------- */
      /* allocate the matrix itself */
      /* ---------------------------------------------------------------------- */
  
      /* allocate O(ncol) space */
      A->p = CHOLMOD(malloc) (ncol+1, sizeof (Int), Common) ;
      if (!packed)
      {
  	A->nz = CHOLMOD(malloc) (ncol, sizeof (Int), Common) ;
      }
  
      /* allocate O(nz) space */
      nzmax0 = 0 ;
      CHOLMOD(realloc_multiple) (nzmax, 1, xtype, &(A->i), NULL, &(A->x), &(A->z),
  	    &nzmax0, Common) ;
  
      if (Common->status < CHOLMOD_OK)
      {
  	CHOLMOD(free_sparse) (&A, Common) ;
  	return (NULL) ;	    /* out of memory */
      }
  
      /* ---------------------------------------------------------------------- */
      /* initialize A->p and A->nz so that A is an empty matrix */
      /* ---------------------------------------------------------------------- */
  
      Ap = A->p ;
      for (j = 0 ; j <= (Int) ncol ; j++)
      {
  	Ap [j] = 0 ;
      }
      if (!packed)
      {
  	Anz = A->nz ;
  	for (j = 0 ; j < (Int) ncol ; j++)
  	{
  	    Anz [j] = 0 ;
  	}
      }
      return (A) ;
  }
  
  
  /* ========================================================================== */
  /* === cholmod_free_sparse ================================================== */
  /* ========================================================================== */
  
  /* free a sparse matrix
   *
   * workspace: none
   */
  
  int CHOLMOD(free_sparse)
  (
      /* ---- in/out --- */
      cholmod_sparse **AHandle,	/* matrix to deallocate, NULL on output */
      /* --------------- */
      cholmod_common *Common
  )
  {
      Int n, nz ;
      cholmod_sparse *A ;
  
      RETURN_IF_NULL_COMMON (FALSE) ;
  
      if (AHandle == NULL)
      {
  	/* nothing to do */
  	return (TRUE) ;
      }
      A = *AHandle ;
      if (A == NULL)
      {
  	/* nothing to do */
  	return (TRUE) ;
      }
      n = A->ncol ;
      nz = A->nzmax ;
      A->p  = CHOLMOD(free) (n+1, sizeof (Int), A->p,  Common) ;
      A->i  = CHOLMOD(free) (nz,  sizeof (Int), A->i,  Common) ;
      A->nz = CHOLMOD(free) (n,   sizeof (Int), A->nz, Common) ;
  
      switch (A->xtype)
      {
  	case CHOLMOD_REAL:
  	    A->x = CHOLMOD(free) (nz, sizeof (double), A->x,  Common) ;
  	    break ;
  
  	case CHOLMOD_COMPLEX:
  	    A->x = CHOLMOD(free) (nz, 2*sizeof (double), A->x,  Common) ;
  	    break ;
  
  	case CHOLMOD_ZOMPLEX:
  	    A->x = CHOLMOD(free) (nz, sizeof (double), A->x,  Common) ;
  	    A->z = CHOLMOD(free) (nz, sizeof (double), A->z,  Common) ;
  	    break ;
      }
  
      *AHandle = CHOLMOD(free) (1, sizeof (cholmod_sparse), (*AHandle), Common) ;
      return (TRUE) ;
  }
  
  
  /* ========================================================================== */
  /* === cholmod_reallocate_sparse ============================================ */
  /* ========================================================================== */
  
  /* Change the size of A->i and A->x, or allocate them if their current size
   * is zero.  A->x is not modified if A->xtype is CHOLMOD_PATTERN.
   * 
   * workspace: none
   */
  
  int CHOLMOD(reallocate_sparse)
  (
      /* ---- input ---- */
      size_t nznew,	/* new # of entries in A */
      /* ---- in/out --- */
      cholmod_sparse *A,	/* matrix to reallocate */
      /* --------------- */
      cholmod_common *Common
  )
  {
  
      /* ---------------------------------------------------------------------- */
      /* get inputs */
      /* ---------------------------------------------------------------------- */
  
      RETURN_IF_NULL_COMMON (FALSE) ;
      RETURN_IF_NULL (A, FALSE) ;
      RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
      Common->status = CHOLMOD_OK ;
      PRINT1 (("realloc matrix %d to %d, xtype: %d\n",
  		A->nzmax, nznew, A->xtype)) ;
  
      /* ---------------------------------------------------------------------- */
      /* resize the matrix */
      /* ---------------------------------------------------------------------- */
  
      CHOLMOD(realloc_multiple) (MAX (1,nznew), 1, A->xtype, &(A->i), NULL,
  	    &(A->x), &(A->z), &(A->nzmax), Common) ;
  
      return (Common->status == CHOLMOD_OK) ;
  }
  
  
  /* ========================================================================== */
  /* === cholmod_speye ======================================================== */
  /* ========================================================================== */
  
  /* Return a sparse identity matrix. */
  
  cholmod_sparse *CHOLMOD(speye)
  (
      /* ---- input ---- */
      size_t nrow,	/* # of rows of A */
      size_t ncol,	/* # of columns of A */
      int xtype,		/* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
      /* --------------- */
      cholmod_common *Common
  )
  {
      double *Ax, *Az ;
      cholmod_sparse *A ;
      Int *Ap, *Ai ;
      Int j, n ;
  
      /* ---------------------------------------------------------------------- */
      /* get inputs */
      /* ---------------------------------------------------------------------- */
  
      RETURN_IF_NULL_COMMON (NULL) ;
      Common->status = CHOLMOD_OK ;
  
      /* ---------------------------------------------------------------------- */
      /* allocate the matrix */
      /* ---------------------------------------------------------------------- */
  
      n = MIN (nrow, ncol) ;
      A = CHOLMOD(allocate_sparse) (nrow, ncol, n, TRUE, TRUE, 0, xtype,
  	    Common) ;
  
      if (Common->status < CHOLMOD_OK)
      {
  	return (NULL) ;	    /* out of memory or inputs invalid */
      }
  
      /* ---------------------------------------------------------------------- */
      /* create the identity matrix */
      /* ---------------------------------------------------------------------- */
  
      Ap = A->p ;
      Ai = A->i ;
      Ax = A->x ;
      Az = A->z ;
  
      for (j = 0 ; j < n ; j++)
      {
  	Ap [j] = j ;
      }
      for (j = n ; j <= ((Int) ncol) ; j++)
      {
  	Ap [j] = n ;
      }
      for (j = 0 ; j < n ; j++)
      {
  	Ai [j] = j ;
      }
  
      switch (xtype)
      {
  	case CHOLMOD_REAL:
  	    for (j = 0 ; j < n ; j++)
  	    {
  		Ax [j] = 1 ;
  	    }
  	    break ;
  
  	case CHOLMOD_COMPLEX:
  	    for (j = 0 ; j < n ; j++)
  	    {
  		Ax [2*j  ] = 1 ;
  		Ax [2*j+1] = 0 ;
  	    }
  	    break ;
  
  	case CHOLMOD_ZOMPLEX:
  	    for (j = 0 ; j < n ; j++)
  	    {
  		Ax [j] = 1 ;
  	    }
  	    for (j = 0 ; j < n ; j++)
  	    {
  		Az [j] = 0 ;
  	    }
  	    break ;
      }
  
      return (A) ;
  }
  
  
  /* ========================================================================== */
  /* === cholmod_spzeros ====================================================== */
  /* ========================================================================== */
  
  /* Return a sparse zero matrix. */
  
  cholmod_sparse *CHOLMOD(spzeros)
  (
      /* ---- input ---- */
      size_t nrow,	/* # of rows of A */
      size_t ncol,	/* # of columns of A */
      size_t nzmax,	/* max # of nonzeros of A */
      int xtype,		/* CHOLMOD_PATTERN, _REAL, _COMPLEX, or _ZOMPLEX */
      /* --------------- */
      cholmod_common *Common
  )
  {
  
      /* ---------------------------------------------------------------------- */
      /* get inputs */
      /* ---------------------------------------------------------------------- */
  
      RETURN_IF_NULL_COMMON (NULL) ;
      Common->status = CHOLMOD_OK ;
  
      /* ---------------------------------------------------------------------- */
      /* allocate the matrix */
      /* ---------------------------------------------------------------------- */
  
      return (CHOLMOD(allocate_sparse) (nrow, ncol, nzmax, TRUE, TRUE, 0, xtype,
  	    Common)) ;
  }
  
  
  /* ========================================================================== */
  /* === cholmod_nnz ========================================================== */
  /* ========================================================================== */
  
  /* Return the number of entries in a sparse matrix.
   *
   * workspace: none
   * integer overflow cannot occur, since the matrix is already allocated.
   */
  
  long CHOLMOD(nnz)
  (
      /* ---- input ---- */
      cholmod_sparse *A,
      /* --------------- */
      cholmod_common *Common
  )
  {
      Int *Ap, *Anz ;
      size_t nz ;
      Int j, ncol ;
  
      /* ---------------------------------------------------------------------- */
      /* get inputs */
      /* ---------------------------------------------------------------------- */
  
      RETURN_IF_NULL_COMMON (EMPTY) ;
      RETURN_IF_NULL (A, EMPTY) ;
      RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, EMPTY) ;
      Common->status = CHOLMOD_OK ;
  
      /* ---------------------------------------------------------------------- */
      /* return nnz (A) */
      /* ---------------------------------------------------------------------- */
  
      ncol = A->ncol ;
      if (A->packed)
      {
  	Ap = A->p ;
  	RETURN_IF_NULL (Ap, EMPTY) ;
  	nz = Ap [ncol] ;
      }
      else
      {
  	Anz = A->nz ;
  	RETURN_IF_NULL (Anz, EMPTY) ;
  	nz = 0 ;
  	for (j = 0 ; j < ncol ; j++)
  	{
  	    nz += MAX (0, Anz [j]) ;
  	}
      }
      return (nz) ;
  }
  
  
  /* ========================================================================== */
  /* === cholmod_copy_sparse ================================================== */
  /* ========================================================================== */
  
  /* C = A.  Create an exact copy of a sparse matrix, with one exception.
   * Entries in unused space are not copied (they might not be initialized,
   * and copying them would cause program checkers such as purify and
   * valgrind to complain).  The xtype of the resulting matrix C is the same as
   * the xtype of the input matrix A.
   *
   * See also Core/cholmod_copy, which copies a matrix with possible changes
   * in stype, presence of diagonal entries, pattern vs. numerical values,
   * real and/or imaginary parts, and so on.
   */
  
  cholmod_sparse *CHOLMOD(copy_sparse)
  (
      /* ---- input ---- */
      cholmod_sparse *A,	/* matrix to copy */
      /* --------------- */
      cholmod_common *Common
  )
  {
      double *Ax, *Cx, *Az, *Cz ;
      Int *Ap, *Ai, *Anz, *Cp, *Ci, *Cnz ;
      cholmod_sparse *C ;
      Int p, pend, j, ncol, packed, nzmax, nz, xtype ;
  
      /* ---------------------------------------------------------------------- */
      /* check inputs */
      /* ---------------------------------------------------------------------- */
  
      RETURN_IF_NULL_COMMON (NULL) ;
      RETURN_IF_NULL (A, NULL) ;
      RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
      if (A->stype != 0 && A->nrow != A->ncol)
      {
  	ERROR (CHOLMOD_INVALID, "rectangular matrix with stype != 0 invalid") ;
  	return (NULL) ;
      }
      Common->status = CHOLMOD_OK ;
      ASSERT (CHOLMOD(dump_sparse) (A, "A original", Common) >= 0) ;
  
      /* ---------------------------------------------------------------------- */
      /* get inputs */
      /* ---------------------------------------------------------------------- */
  
      ncol = A->ncol ;
      nzmax = A->nzmax ;
      packed = A->packed ;
      Ap = A->p ;
      Ai = A->i ;
      Ax = A->x ;
      Az = A->z ;
      Anz = A->nz ;
      xtype = A->xtype ;
  
      /* ---------------------------------------------------------------------- */
      /* allocate the copy */
      /* ---------------------------------------------------------------------- */
  
      C = CHOLMOD(allocate_sparse) (A->nrow, A->ncol, A->nzmax, A->sorted,
  	    A->packed, A->stype, A->xtype, Common) ;
  
      if (Common->status < CHOLMOD_OK)
      {
  	return (NULL) ;	    /* out of memory */
      }
  
      Cp = C->p ;
      Ci = C->i ;
      Cx = C->x ;
      Cz = C->z ;
      Cnz = C->nz ;
  
      /* ---------------------------------------------------------------------- */
      /* copy the matrix */
      /* ---------------------------------------------------------------------- */
  
      for (j = 0 ; j <= ncol ; j++)
      {
  	Cp [j] = Ap [j] ;
      }
  
      if (packed)
      {
  	nz = Ap [ncol] ;
  	for (p = 0 ; p < nz ; p++)
  	{
  	    Ci [p] = Ai [p] ;
  	}
  
  	switch (xtype)
  	{
  	    case CHOLMOD_REAL:
  		for (p = 0 ; p < nz ; p++)
  		{
  		    Cx [p] = Ax [p] ;
  		}
  		break ;
  
  	    case CHOLMOD_COMPLEX:
  		for (p = 0 ; p < 2*nz ; p++)
  		{
  		    Cx [p] = Ax [p] ;
  		}
  		break ;
  
  	    case CHOLMOD_ZOMPLEX:
  		for (p = 0 ; p < nz ; p++)
  		{
  		    Cx [p] = Ax [p] ;
  		    Cz [p] = Az [p] ;
  		}
  		break ;
  	}
  
      }
      else
      {
  
  	for (j = 0 ; j < ncol ; j++)
  	{
  	    Cnz [j] = Anz [j] ;
  	}
  
  	switch (xtype)
  	{
  	    case CHOLMOD_PATTERN:
  		for (j = 0 ; j < ncol ; j++)
  		{
  		    p = Ap [j] ;
  		    pend = p + Anz [j] ;
  		    for ( ; p < pend ; p++)
  		    {
  			Ci [p] = Ai [p] ;
  		    }
  		}
  		break ;
  
  	    case CHOLMOD_REAL:
  		for (j = 0 ; j < ncol ; j++)
  		{
  		    p = Ap [j] ;
  		    pend = p + Anz [j] ;
  		    for ( ; p < pend ; p++)
  		    {
  			Ci [p] = Ai [p] ;
  			Cx [p] = Ax [p] ;
  		    }
  		}
  		break ;
  
  	    case CHOLMOD_COMPLEX:
  		for (j = 0 ; j < ncol ; j++)
  		{
  		    p = Ap [j] ;
  		    pend = p + Anz [j] ;
  		    for ( ; p < pend ; p++)
  		    {
  			Ci [p] = Ai [p] ;
  			Cx [2*p  ] = Ax [2*p  ] ;
  			Cx [2*p+1] = Ax [2*p+1] ;
  		    }
  		}
  		break ;
  
  	    case CHOLMOD_ZOMPLEX:
  		for (j = 0 ; j < ncol ; j++)
  		{
  		    p = Ap [j] ;
  		    pend = p + Anz [j] ;
  		    for ( ; p < pend ; p++)
  		    {
  			Ci [p] = Ai [p] ;
  			Cx [p] = Ax [p] ;
  			Cz [p] = Az [p] ;
  		    }
  		}
  		break ;
  	}
      }
  
      /* ---------------------------------------------------------------------- */
      /* return the result */
      /* ---------------------------------------------------------------------- */
  
      ASSERT (CHOLMOD(dump_sparse) (C, "C copy", Common) >= 0) ;
      return (C) ;
  }
  