matthiasm@8: 
matthiasm@8:   m  = mxGetM(prhs[0]);
matthiasm@8:   n  = mxGetN(prhs[0]);
matthiasm@8:   pr = mxGetPr(prhs[0]);
matthiasm@8:   pi = mxGetPi(prhs[0]);
matthiasm@8:   cmplx = (pi == NULL ? 0 : 1);
matthiasm@8: 
matthiasm@8:   /* Allocate space for sparse matrix. 
matthiasm@8:    * NOTE:  Assume at most 20% of the data is sparse.  Use ceil
matthiasm@8:    * to cause it to round up. 
matthiasm@8:    */
matthiasm@8: 
matthiasm@8:   percent_sparse = 0.2;
matthiasm@8:   nzmax = (int)ceil((double)m*(double)n*percent_sparse);
matthiasm@8: 
matthiasm@8:   plhs[0] = mxCreateSparse(m,n,nzmax,cmplx);
matthiasm@8:   sr  = mxGetPr(plhs[0]);
matthiasm@8:   si  = mxGetPi(plhs[0]);
matthiasm@8:   irs = mxGetIr(plhs[0]);
matthiasm@8:   jcs = mxGetJc(plhs[0]);
matthiasm@8:     
matthiasm@8:   /* Copy nonzeros. */
matthiasm@8:   k = 0; 
matthiasm@8:   isfull = 0;
matthiasm@8:   for (j = 0; (j < n); j++) {
matthiasm@8:     int i;
matthiasm@8:     jcs[j] = k;
matthiasm@8:     for (i = 0; (i < m); i++) {
matthiasm@8:       if (IsNonZero(pr[i]) || (cmplx && IsNonZero(pi[i]))) {
matthiasm@8: 
matthiasm@8:         /* Check to see if non-zero element will fit in 
matthiasm@8:          * allocated output array.  If not, increase
matthiasm@8:          * percent_sparse by 10%, recalculate nzmax, and augment
matthiasm@8:          * the sparse array.
matthiasm@8:          */
matthiasm@8:         if (k >= nzmax) {
matthiasm@8:           int oldnzmax = nzmax;
matthiasm@8:           percent_sparse += 0.1;
matthiasm@8:           nzmax = (int)ceil((double)m*(double)n*percent_sparse);
matthiasm@8: 
matthiasm@8:           /* Make sure nzmax increases atleast by 1. */
matthiasm@8:           if (oldnzmax == nzmax) 
matthiasm@8:             nzmax++;
matthiasm@8: 
matthiasm@8:           mxSetNzmax(plhs[0], nzmax); 
matthiasm@8:           mxSetPr(plhs[0], mxRealloc(sr, nzmax*sizeof(double)));
matthiasm@8:           if (si != NULL)
matthiasm@8:           mxSetPi(plhs[0], mxRealloc(si, nzmax*sizeof(double)));
matthiasm@8:           mxSetIr(plhs[0], mxRealloc(irs, nzmax*sizeof(int)));
matthiasm@8: 
matthiasm@8:           sr  = mxGetPr(plhs[0]);
matthiasm@8:           si  = mxGetPi(plhs[0]);
matthiasm@8:           irs = mxGetIr(plhs[0]);
matthiasm@8:         }
matthiasm@8:         sr[k] = pr[i];
matthiasm@8:         if (cmplx) {
matthiasm@8:           si[k] = pi[i];
matthiasm@8:         }
matthiasm@8:         irs[k] = i;
matthiasm@8:         k++;
matthiasm@8:       }
matthiasm@8:     }
matthiasm@8:     pr += m;
matthiasm@8:     pi += m;
matthiasm@8:   }
matthiasm@8:   jcs[n] = k;
matthiasm@8: }