External Interfaces/API

Handling Sparse Matrices

The MATLAB API provides a set of functions that allow you to create and manipulate sparse matrices from within your MEX-files. There are special parameters associated with sparse matrices, namely ir, jc, and nzmax. For information on how to use these parameters and how MATLAB stores sparse matrices in general, refer to the section on The MATLAB Array.

Note    Sparse array indexing is zero-based, not one-based.

This example illustrates how to populate a sparse matrix.

• C==============================================================
  =
  C     fulltosparse.f
  C     Example for illustrating how to populate a sparse matrix.
  C     For the purpose of this example, you must pass in a
  C     non-sparse 2-dimensional argument of type real double.
  C
  C     NOTE: The subroutine loadsparse() is in the file called 
  C     loadsparse.f. 
  C
  C     This is a MEX-file for MATLAB.
  C     Copyright (c) 1984-2000 The MathWorks, Inc. 
  C     $Revision: 1.9 $
  C==============================================================
  =
  
  C     The gateway routine.
        subroutine mexFunction(nlhs, plhs, nrhs, prhs)
  C--------------------------------------------------------------
  -
  C     (pointer) Replace integer by integer*8 on the DEC Alpha
  C     64-bit platform
  
        integer plhs(*), prhs(*)
        integer mxGetPr, mxCreateSparse, mxGetIr, mxGetJc
        integer pr, sr, irs, jcs
  C--------------------------------------------------------------
  -
  
        integer nlhs, nrhs
        integer m, n, nzmax
        integer mxGetM, mxGetN, mxIsComplex, mxIsDouble
        integer loadsparse
  
  C     Check for proper number of arguments. 
        if (nrhs .ne. 1) then
           call mexErrMsgTxt('One input argument required.')
        endif
        if (nlhs .gt. 1) then
           call mexErrMsgTxt('Too many output arguments.')
        endif
  
  C     Check data type of input argument.
        if (mxIsDouble(prhs(1)) .eq. 0) then
           call mexErrMsgTxt('Input argument must be of type 
  double.')
        endif
        if (mxIsComplex(prhs(1)) .eq. 1) then
           call mexErrMsgTxt('Input argument must be real only')
        endif
  
  C     Get the size and pointers to input data.
        m = mxGetM(prhs(1))
        n = mxGetN(prhs(1))
        pr = mxGetPr(prhs(1))
  
  C     Allocate space.
  C     NOTE: Assume at most 20% of the data is sparse.
        nzmax = dble(m*n) *.20 + .5
  
  C     NOTE: The maximum number of non-zero elements cannot be less
  C     than the number of columns in the matrix.
        if (n .gt. nzmax) then
           nzmax = n
        endif
  
        plhs(1) = mxCreateSparse(m,n,nzmax,0)
        sr = mxGetPr(plhs(1))
        irs = mxGetIr(plhs(1))
        jcs = mxGetJc(plhs(1))
  
  C     Load the sparse data.
        if (loadsparse(%val(pr), %val(sr), %val(irs), %val(jcs),
       +m,n,nzmax) .eq. 1) then
           call mexPrintf('Truncating output, input is > 20%% 
  sparse')
        endif
        return
        end
  

This is the subroutine that fulltosparse calls to fill the mxArray with the sparse data.

• C==============================================================
  =
  C     loadsparse.f
  C     This is the subfunction called by fulltosparse that fills the 
  C     mxArray with the sparse data.  Your version of
  C     loadsparse can operate however you would like it to on the
  C     data.
  C
  C     This is a MEX-file for MATLAB.
  C     Copyright (c) 1984-2000 The MathWorks, Inc. 
  C     $Revision: 1.9 $
  C==============================================================
  =
  
  C     Load sparse data subroutine.
        function loadsparse(a, b, ir, jc, m, n, nzmax)
        integer nzmax, m, n
        integer ir(*), jc(*)
        real*8 a(*), b(*)
        integer i, j, k
  
  C     Copy nonzeros.
        k = 1
        do 100 j=1,n
  C        NOTE: Sparse indexing is zero based.
           jc(j) = k-1
           do 200 i=1,m
              if (a((j-1)*m+i).ne. 0.0) then
                 if (k .gt. nzmax) then
                    jc(n+1) = nzmax
                    loadsparse = 1
                    goto 300
                 endif
                 b(k) = a((j-1)*m+i)
  C              NOTE: Sparse indexing is zero based.
                 ir(k) = i-1
                 k = k+1
              endif
   200     continue
   100  continue
  C     NOTE: Sparse indexing is zero based.
        jc(n+1) = k-1
        loadsparse = 0
   300  return
        end
  

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
  

Dynamically Allocating Memory

Calling Functions from Fortran MEX-Files