/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*//*             ********   ***                                 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 voidCompColMatVec_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 voidCompRowMatVec_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 voidCompRowMatVec_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 voidScaleRectangularArray_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 voidScaleRectangularArray_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