THESE/codes/graphe/Ncut_9/spmtimesd.cpp

   1 /*================================================================
   2 * spmtimesd.c
   3 * This routine computes a sparse matrix times a diagonal matrix
   4 * whose diagonal entries are stored in a full vector.
   5 *
   6 * Examples:
   7 *     spmtimesd(m,d,[]) = diag(d) * m,
   8 *     spmtimesd(m,[],d) = m * diag(d)
   9 *     spmtimesd(m,d1,d2) = diag(d1) * m * diag(d2)
  10 *     m could be complex, but d is assumed real
  11 *
  12 * Stella X. Yu's first MEX function, Nov 9, 2001.
  13
  14 % test sequence:
  15
  16 m = 1000;
  17 n = 2000;
  18 a=sparse(rand(m,n));
  19 d1 = rand(m,1);
  20 d2 = rand(n,1);
  21 tic; b=spmtimesd(a,d1,d2); toc
  22 tic; bb = spdiags(d1,0,m,m) * a * spdiags(d2,0,n,n); toc
  23 e = (bb-b);
  24 max(abs(e(:)))
  25
  26 *=================================================================*/
  27
  28 # include "mex.h"
  29
  30 void mexFunction(
  31     int nargout,
  32     mxArray *out[],
  33     int nargin,
  34     const mxArray *in[]
  35 )
  36 {
  37     /* declare variables */
  38     int i, j, k, m, n, nzmax, xm, yn;
  39     mwIndex *pir, *pjc, *qir, *qjc;
  40     double *x, *y, *pr, *pi, *qr, *qi;
  41
  42     /* check argument */
  43     if (nargin != 3) {
  44         mexErrMsgTxt("Three input arguments required");
  45     }
  46     if (nargout>1) {
  47         mexErrMsgTxt("Too many output arguments.");
  48     }
  49     if (!(mxIsSparse(in[0]))) {
  50         mexErrMsgTxt("Input argument #1 must be of type sparse");
  51     }
  52     if ( mxIsSparse(in[1]) || mxIsSparse(in[2]) ) {
  53         mexErrMsgTxt("Input argument #2 & #3 must be of type full");
  54     }
  55
  56     /* computation starts */
  57     m = mxGetM(in[0]);
  58     n = mxGetN(in[0]);
  59     pr = mxGetPr(in[0]);
  60     pi = mxGetPi(in[0]);
  61     pir = mxGetIr(in[0]);
  62     pjc = mxGetJc(in[0]);
  63
  64     i = mxGetM(in[1]);
  65     j = mxGetN(in[1]);
  66     xm = ((i>j)? i: j);
  67
  68     i = mxGetM(in[2]);
  69     j = mxGetN(in[2]);
  70     yn = ((i>j)? i: j);
  71
  72     if ( xm>0 && xm != m) {
  73         mexErrMsgTxt("Row multiplication dimension mismatch.");
  74     }
  75     if ( yn>0 && yn != n) {
  76         mexErrMsgTxt("Column multiplication dimension mismatch.");
  77     }
  78
  79
  80     nzmax = mxGetNzmax(in[0]);
  81     mxComplexity cmplx = (pi==NULL ? mxREAL : mxCOMPLEX);
  82     out[0] = mxCreateSparse(m,n,nzmax,cmplx);
  83     if (out[0]==NULL) {
  84         mexErrMsgTxt("Not enough space for the output matrix.");
  85     }
  86
  87     qr = mxGetPr(out[0]);
  88     qi = mxGetPi(out[0]);
  89     qir = mxGetIr(out[0]);
  90     qjc = mxGetJc(out[0]);
  91
  92     /* left multiplication */
  93     x = mxGetPr(in[1]);
  94     if (yn==0) {
  95         for (j=0; j<n; j++) {
  96             qjc[j] = pjc[j];
  97             for (k=pjc[j]; k<pjc[j+1]; k++) {
  98                 i = pir[k];
  99                 qir[k] = i;
 100                  qr[k] = x[i] * pr[k];
 101                  if (cmplx) {
 102                      qi[k] = x[i] * pi[k];
 103                  }
 104             }
 105         }
 106         qjc[n] = k;
 107         return;
 108     }
 109
 110     /* right multiplication */
 111     y = mxGetPr(in[2]);
 112     if (xm==0) {
 113         for (j=0; j<n; j++) {
 114             qjc[j] = pjc[j];
 115             for (k=pjc[j]; k<pjc[j+1]; k++) {
 116                 qir[k] = pir[k];
 117                 qr[k]  = pr[k] * y[j];
 118                 if (cmplx) {
 119                     qi[k] = qi[k] * y[j];
 120                 }
 121             }
 122        }
 123        qjc[n] = k;
 124        return;
 125    }
 126
 127    /* both sides */
 128    for (j=0; j<n; j++) {
 129        qjc[j] = pjc[j];
 130        for (k=pjc[j]; k<pjc[j+1]; k++) {
 131            i = pir[k];
 132            qir[k]= i;
 133            qr[k] = x[i] * pr[k] * y[j];
 134            if (cmplx) {
 135                qi[k] = x[i] * qi[k] * y[j];
 136            }
 137        }
 138        qjc[n] = k;
 139    }
 140
 141 }