/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*//*             ********   ***                                 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)])/* * int &m   Number of rows in matrix c * * int &n   Number of columns in matrix c * * int &k   Number of rows in matrix b * * Assume diagonal elements are in proper place * * unitd = 1    D = I * unitd = 2    left (row scaling) * unitd = 3    right (column scaling) */static voidCopyRectangularArray_double(int m, int n,             const double *b, int ldb, double *c, int ldc){  int i, l;  if (b == c)    return;  if (n == 1)    for (i = 0; i < m; i++)      c[i] = b[i];  else    for (l = 0; l < n; l++)      for (i = 0; i < m; i++)    _SpMatVal(c, ldc, i, l) = _SpMatVal(b, ldb, i, l);}static voidCopyRectangularArray_float(int m, int n,             const float *b, int ldb, float *c, int ldc){  int i, l;  if (b == c)    return;  if (n == 1)    for (i = 0; i < m; i++)      c[i] = b[i];  else    for (l = 0; l < n; l++)      for (i = 0; i < m; i++)    _SpMatVal(c, ldc, i, l) = _SpMatVal(b, ldb, i, l);}static voidCompCol_LowerUnitSolve_double(int m, int n, int unitd, const double *dv,       double alpha, const double *val, const int *indx, const int *pntr,      const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = alpha * c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_LowerUnitSolve_float(int m, int n, int unitd, const float *dv,     float alpha, const float *val, const int *indx, const int *pntr,    const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = alpha * c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_LowerDiagSolve_double(int m, int n, int unitd, const double *dv, double alpha,               const double *val, const int *indx, const int *pntr,               const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = c[i+pntr[0]] / val[pntr[i]];    c[i+pntr[0]] = z;    for (j = pntr[i] + 1; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];      _SpMatVal(c, ldc, i+pntr[0], l) = z;      for (j = pntr[i] + 1; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = alpha * c[i+pntr[0]] / val[pntr[i]];    c[i+pntr[0]] = z;    for (j = pntr[i] + 1; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];      _SpMatVal(c, ldc, i, l) = z;      for (j = pntr[i] + 1; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_LowerDiagSolve_float(int m, int n, int unitd, const float *dv,                 float alpha,               const float *val, const int *indx, const int *pntr,               const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = c[i+pntr[0]] / val[pntr[i]];    c[i+pntr[0]] = z;    for (j = pntr[i] + 1; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];      _SpMatVal(c, ldc, i+pntr[0], l) = z;      for (j = pntr[i] + 1; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = alpha * c[i+pntr[0]] / val[pntr[i]];    c[i+pntr[0]] = z;    for (j = pntr[i] + 1; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];      _SpMatVal(c, ldc, i, l) = z;      for (j = pntr[i] + 1; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_UpperUnitSolve_double(int m, int n, int unitd, const double *dv,                 double alpha,               const double *val, const int *indx, const int *pntr,               const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];    indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = alpha * c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_UpperUnitSolve_float(int m, int n, int unitd, const float *dv,                 float alpha,               const float *val, const int *indx, const int *pntr,               const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];    indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = alpha * c[i+pntr[0]];    for (j = pntr[i]; j < pntr[i+1]; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);      for (j = pntr[i]; j < pntr[i+1]; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_UpperDiagSolve_double(int m, int n, int unitd, const double *dv,             double alpha,               const double *val, const int *indx, const int *pntr,               const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {