spmm.cc

http://gams.cam.nist.gov/acmd/Staff/RPozo/sparselib++.html

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*             ********   ***                                 SparseLib++    */
/*          *******  **  ***       ***      ***                              */
/*           *****      ***     ******** ********                            */
/*            *****    ***     ******** ********              R. Pozo        */
/*       **  *******  ***   **   ***      ***                 K. Remington   */
/*        ********   ********                                 A. Lumsdaine   */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*                                                                           */
/*                                                                           */
/*                     SparseLib++ : Sparse Matrix Library                   */
/*                                                                           */
/*               National Institute of Standards and Technology              */
/*                        University of Notre Dame                           */
/*              Authors: R. Pozo, K. Remington, A. Lumsdaine                 */
/*                                                                           */
/*                                 NOTICE                                    */
/*                                                                           */
/* Permission to use, copy, modify, and distribute this software and         */
/* its documentation for any purpose and without fee is hereby granted       */
/* provided that the above notice appear in all copies and supporting        */
/* documentation.                                                            */
/*                                                                           */
/* Neither the Institutions (National Institute of Standards and Technology, */
/* University of Notre Dame) nor the Authors make any representations about  */
/* the suitability of this software for any purpose.  This software is       */
/* provided ``as is'' without expressed or implied warranty.                 */
/*                                                                           */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/*
 * Home Grown Sparse BLAS
 *
 * These are just a subset of the functions described in SPARKER
 * Working Note #3.
 *
 * Would be great if these could be templated some day
 *
 */

#include <stdlib.h>
#include <iostream>
#include "spblas.h"

#define _SpMatVal(_a,_lda,_row,_col) ((_a)[(_lda)*(_col)+(_row)])

static void CoordMatVec_float(int m, int n, int k, const float &alpha,
        const float *val, const int *indx, const int *jndx,
        const int &nnz,
        const float *b, int ldb, float *c, int ldc)
{
  int i, j;

  // To make the compiler happy
  if (k && m)
    ;

  // Frob these so we can use one-based indexing externally
  b -= 1;
  c -= 1;

  if (alpha == 1.0) {
    if (n == 1)
      for (j = 0; j < nnz; j++)
    c[indx[j]] += b[jndx[j]] * val[j];
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < nnz; j++)
      _SpMatVal(c, ldc, indx[j], i) += _SpMatVal(b, ldb, indx[j], i) * val[j];
  } else {
    if (n == 1)
      for (j = 0; j < nnz; j++)
    c[indx[j]] += alpha * b[jndx[j]] * val[j];
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < nnz; j++)
      _SpMatVal(c, ldc, indx[j], i) +=
        alpha * _SpMatVal(b, ldb, indx[j], i) * val[j];
  }
}

static void CoordMatVec_double(int m, int n, int k, const double &alpha,
        const double *val, const int *indx, const int *jndx,
        const int &nnz,
        const double *b, int ldb, double *c, int ldc)
{
  int i, j;

  // To make the compiler happy
  if (k && m)
    ;

  // Frob these so we can use one-based indexing externally
  b -= 1;
  c -= 1;

  if (alpha == 1.0) {
    if (n == 1)
      for (j = 0; j < nnz; j++)
    c[indx[j]] += b[jndx[j]] * val[j];
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < nnz; j++)
      _SpMatVal(c, ldc, indx[j], i) += _SpMatVal(b, ldb, indx[j], i) * val[j];
  } else {
    if (n == 1)
      for (j = 0; j < nnz; j++)
    c[indx[j]] += alpha * b[jndx[j]] * val[j];
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < nnz; j++)
      _SpMatVal(c, ldc, indx[j], i) +=
        alpha * _SpMatVal(b, ldb, indx[j], i) * val[j];
  }
}

static void
CompColMatVec_double(int m, int n, int k, const double &alpha,
            const double *val, const int *indx, const int *pntr,
            const double *b, int ldb, double *c, int ldc)
{
  int i, j, l;

  if (alpha == 0.0)
    return;

  // To make the compiler happy
  if (m)
    ;

  // Frob these so we can use one-based indexing externally
  c -= 1;
  val -= pntr[0];
  indx -= pntr[0];

  if (alpha == 1.0) {
    if (n == 1)
      for (i = 0; i < k; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[indx[j]] += b[i] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < k; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, indx[j], l) += _SpMatVal(b, ldb, i, l) * val[j];
  } else {
    if (n == 1)
      for (i = 0; i < k; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[indx[j]] += alpha * b[i] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < k; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, indx[j], l) +=
          alpha * _SpMatVal(b, ldb, i, l) * val[j];
  }
}

static void CompColMatVec_float(int m, int n, int k, const float &alpha,
            const float *val, const int *indx, const int *pntr,
            const float *b, int ldb, float *c, int ldc)
{
  int i, j, l;

  if (alpha == 0.0)
    return;

  // To make the compiler happy
  if (m)
    ;

  // Frob these so we can use one-based indexing externally
  c -= 1;
  val -= pntr[0];
  indx -= pntr[0];

  if (alpha == 1.0) {
    if (n == 1)
      for (i = 0; i < k; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[indx[j]] += b[i] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < k; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, indx[j], l) += _SpMatVal(b, ldb, i, l) * val[j];
  } else {
    if (n == 1)
      for (i = 0; i < k; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[indx[j]] += alpha * b[i] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < k; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, indx[j], l) +=
          alpha * _SpMatVal(b, ldb, i, l) * val[j];
  }
}


static void
CompRowMatVec_double(int m, int n, int k, const double &alpha,
            const double *val, const int *indx, const int *pntr,
            const double *b, int ldb, double *c, int ldc)
{
  int i, j, l;

  if (alpha == 0.0)
    return;

  // To make the compiler happy
  if (m || k)
    ;

  // Frob these so we can use one-based indexing externally
  b -= 1;
  val -= pntr[0];
  indx -= pntr[0];

  if (alpha == 1.0) {
    if (n == 1)
      for (i = 0; i < m; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[i] += b[indx[j]] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < m; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, i, l) += _SpMatVal(b, ldb, indx[j], l) * val[j];
  } else {
    if (n == 1)
      for (i = 0; i < m; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[i] += alpha * b[indx[j]] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < m; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, i, l) +=
          alpha * _SpMatVal(b, ldb, indx[j], l) * val[j];
  }
}


static void
CompRowMatVec_float(int m, int n, int k, const float &alpha,
            const float *val, const int *indx, const int *pntr,
            const float *b, int ldb, float *c, int ldc)
{
  int i, j, l;

  if (alpha == 0.0)
    return;

  // To make the compiler happy
  if (m || k)
    ;

  // Frob these so we can use one-based indexing externally
  b -= 1;
  val -= pntr[0];
  indx -= pntr[0];

  if (alpha == 1.0) {
    if (n == 1)
      for (i = 0; i < m; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[i] += b[indx[j]] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < m; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, i, l) += _SpMatVal(b, ldb, indx[j], l) * val[j];
  } else {
    if (n == 1)
      for (i = 0; i < m; i++)
    for (j = pntr[i]; j < pntr[i+1]; j++)
      c[i] += alpha * b[indx[j]] * val[j];
    else
      for (l = 0; l < n; l++)
    for (i = 0; i < m; i++)
      for (j = pntr[i]; j < pntr[i+1]; j++)
        _SpMatVal(c, ldc, i, l) +=
          alpha * _SpMatVal(b, ldb, indx[j], l) * val[j];
  }
}

static void
ScaleRectangularArray_double(int m, int n, double *c, int ldc, 
    const double &beta)
{
  int i, j;

  if (beta == 1.0)
    return;

  if (beta == 0.0) {
    if (n == 1)
      for (j = 0; j < m; j++)
    c[j] = 0.0;
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      _SpMatVal(c, ldc, j, i) = 0.0;
  } else {
    if (n == 1)
      for (j = 0; j < m; j++)
    c[j] *= beta;
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      _SpMatVal(c, ldc, j, i) *= beta;
  }
}

static void
ScaleRectangularArray_float(int m, int n, float *c, int ldc, 
    const double &beta)
{
  int i, j;

  if (beta == 1.0)
    return;

  if (beta == 0.0) {
    if (n == 1)
      for (j = 0; j < m; j++)
    c[j] = 0.0;
    else
      for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      _SpMatVal(c, ldc, j, i) = 0.0;
  } else {
    if (n == 1)
      for (j = 0; j < m; j++)
    c[j] *= beta;
    else