Handling Sparse Matrices

最新推荐文章于 2025-12-02 19:45:00 发布

转载最新推荐文章于 2025-12-02 19:45:00 发布 · 984 阅读

文章标签：

#matlab #arrays #input #equals #pointers #output

本文介绍了一个MATLAB MEX文件示例，该示例演示了如何将完整的二维双精度矩阵转换为稀疏矩阵。通过调整稀疏度百分比参数，可以确保所有非零元素都被正确地复制到新的稀疏矩阵中。

Handling Sparse Matrices

The MATLAB API provides a set of functions that allow you to create and manipulate sparse arrays from within your MEX-files. These API routines access and manipulate ir and jc, two of the parameters associated with sparse arrays. For more information on how MATLAB stores sparse arrays, refer to the section, The MATLAB Array.

This example illustrates how to populate a sparse matrix.

/*
 * =============================================================
 * fulltosparse.c
 * This example demonstrates how to populate a sparse
 * matrix.  For the purpose of this example, you must pass in a
 * non-sparse 2-dimensional argument of type double.
 *
 * Comment: You might want to modify this MEX-file so that you can 
 * use it to read large sparse data sets into MATLAB.
 *
 * This is a MEX-file for MATLAB.  
 * Copyright (c) 1984-2000 The MathWorks, Inc.
 * =============================================================
 */

/* $ Revision: 1.5 $ */

#include <math.h> /* Needed for the ceil() prototype. */
#include "mex.h"

/* If you are using a compiler that equates NaN to be zero, you 
 * must compile this example using the flag  -DNAN_EQUALS_ZERO.
 * For example:
 *
 *     mex -DNAN_EQUALS_ZERO fulltosparse.c
 *
 * This will correctly define the IsNonZero macro for your C
 * compiler.
 */

#if defined(NAN_EQUALS_ZERO)
#define IsNonZero(d) ((d) != 0.0 || mxIsNaN(d))
#else
#define IsNonZero(d) ((d) != 0.0)
#endif

void mexFunction(
        int nlhs,       mxArray *plhs[],
        int nrhs, const mxArray *prhs[]
        )
{
  /* Declare variables. */
  int j,k,m,n,nzmax,*irs,*jcs,cmplx,isfull;
  double *pr,*pi,*si,*sr;
  double percent_sparse;

  /* Check for proper number of input and output arguments. */    
  if (nrhs != 1) {
    mexErrMsgTxt("One input argument required.");
  } 
  if (nlhs > 1) {
    mexErrMsgTxt("Too many output arguments.");
  }

  /* Check data type of input argument. */
  if (!(mxIsDouble(prhs[0]))) {
    mexErrMsgTxt("Input argument must be of type double.");
  } 

  if (mxGetNumberOfDimensions(prhs[0]) != 2) {
    mexErrMsgTxt("Input argument must be two dimensional/n");
  }

  /* Get the size and pointers to input data. */
  m  = mxGetM(prhs[0]);
  n  = mxGetN(prhs[0]);
  pr = mxGetPr(prhs[0]);
  pi = mxGetPi(prhs[0]);
  cmplx = (pi == NULL ? 0 : 1);

  /* Allocate space for sparse matrix. 
   * NOTE:  Assume at most 20% of the data is sparse.  Use ceil
   * to cause it to round up. 
   */

  percent_sparse = 0.2;
  nzmax = (int)ceil((double)m*(double)n*percent_sparse);

  plhs[0] = mxCreateSparse(m,n,nzmax,cmplx);
  sr  = mxGetPr(plhs[0]);
  si  = mxGetPi(plhs[0]);
  irs = mxGetIr(plhs[0]);
  jcs = mxGetJc(plhs[0]);
    
  /* Copy nonzeros. */
  k = 0; 
  isfull = 0;
  for (j = 0; (j < n); j++) {
    int i;
    jcs[j] = k;
    for (i = 0; (i < m); i++) {
      if (IsNonZero(pr[i]) || (cmplx && IsNonZero(pi[i]))) {

        /* Check to see if non-zero element will fit in 
         * allocated output array.  If not, increase
         * percent_sparse by 10%, recalculate nzmax, and augment
         * the sparse array.
         */
        if (k >= nzmax) {
          int oldnzmax = nzmax;
          percent_sparse += 0.1;
          nzmax = (int)ceil((double)m*(double)n*percent_sparse);

          /* Make sure nzmax increases atleast by 1. */
          if (oldnzmax == nzmax) 
            nzmax++;

          mxSetNzmax(plhs[0], nzmax); 
          mxSetPr(plhs[0], mxRealloc(sr, nzmax*sizeof(double)));
          if (si != NULL)
          mxSetPi(plhs[0], mxRealloc(si, nzmax*sizeof(double)));
          mxSetIr(plhs[0], mxRealloc(irs, nzmax*sizeof(int)));

          sr  = mxGetPr(plhs[0]);
          si  = mxGetPi(plhs[0]);
          irs = mxGetIr(plhs[0]);
        }
        sr[k] = pr[i];
        if (cmplx) {
          si[k] = pi[i];
        }
        irs[k] = i;
        k++;
      }
    }
    pr += m;
    pi += m;
  }
  jcs[n] = k;
}

At the MATLAB prompt, entering

full = eye(5)
full =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     0     1

creates a full, 5-by-5 identity matrix. Using fulltosparse on the full matrix produces the corresponding sparse matrix.

spar = fulltosparse(full)
spar =
   (1,1)        1
   (2,2)        1
   (3,3)        1
   (4,4)        1
   (5,5)        1