text: /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
  /*             ********   ***                                 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 void
  CopyRectangularArray_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 void
  CopyRectangularArray_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 void
  CompCol_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 void
  CompCol_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 void
  CompCol_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 void
  CompCol_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 void
  CompCol_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 void
  CompCol_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 void
  CompCol_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 )
      {
      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]] / val[pntr[i + 1] - 1];
          c[i + pntr[0]] = z;
          for( j = pntr[i]; j < pntr[i + 1] - 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) / val[pntr[i + 1] - 1];
            _SpMatVal(c, ldc, i + pntr[0], l) = z;
            for( j = pntr[i]; j < pntr[i + 1] - 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]] / val[pntr[i + 1] - 1];
          c[i + pntr[0]] = z;
          for( j = pntr[i]; j < pntr[i + 1] - 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) / val[pntr[i + 1] - 1];
            _SpMatVal(c, ldc, i + pntr[0], l) = z;
            for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
              {
              _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
              }
            }
          }
        }
      }
  }
  
  static void
  CompCol_UpperDiagSolve_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]] / val[pntr[i + 1] - 1];
          c[i + pntr[0]] = z;
          for( j = pntr[i]; j < pntr[i + 1] - 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) / val[pntr[i + 1] - 1];
            _SpMatVal(c, ldc, i + pntr[0], l) = z;
            for( j = pntr[i]; j < pntr[i + 1] - 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]] / val[pntr[i + 1] - 1];
          c[i + pntr[0]] = z;
          for( j = pntr[i]; j < pntr[i + 1] - 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) / val[pntr[i + 1] - 1];
            _SpMatVal(c, ldc, i + pntr[0], l) = z;
            for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
              {
              _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
              }
            }
          }
        }
      }
  }
  
  static void
  CompRow_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = b[i] - z;
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + 1] = alpha * (b[i] - z);
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
  }
  
  static void
  CompRow_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = b[i] - z;
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + 1] = alpha * (b[i] - z);
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
  }
  
  static void
  CompRow_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = (b[i] - z) / val[pntr[i + 1] - 1];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i + 1] - 1];
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = alpha * (b[i] - z) / val[pntr[i + 1] - 1];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i + 1] - 1];
            }
          }
        }
      }
  }
  
  static void
  CompRow_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = (b[i] - z) / val[pntr[i + 1] - 1];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i + 1] - 1];
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = 0; i < m; i++ )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = alpha * (b[i] - z) / val[pntr[i + 1] - 1];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = 0; i < m; i++ )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1] - 1; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i + 1] - 1];
            }
          }
        }
      }
  }
  
  static void
  CompRow_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = b[i] - z;
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = alpha * (b[i] - z);
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
  }
  
  static void
  CompRow_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = b[i] - z;
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i]; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = alpha * (b[i] - z);
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i]; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
            }
          }
        }
      }
  }
  
  static void
  CompRow_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 )
      {
      std::cerr << "unitd != 1 not implemented" << "\n";
      exit(1);
      }
  
    if( alpha == 0.0 )
      {
      return;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = (b[i] - z) / val[pntr[i]];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = alpha * (b[i] - z) / val[pntr[i]];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
            }
          }
        }
      }
  }
  
  static void
  CompRow_UpperDiagSolve_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;
      }
  
    c -= 1;
    val -= pntr[0];
    indx -= pntr[0];
  
    if( alpha == 1.0 )
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = (b[i] - z) / val[pntr[i]];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
            }
          }
        }
      }
    else
      {
      if( n == 1 )
        {
        for( i = m - 1; i >= 0; i-- )
          {
          z = 0;
          for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
            {
            z += c[indx[j]] * val[j];
            }
          c[i + pntr[0]] = alpha * (b[i] - z) / val[pntr[i]];
          }
        }
      else
        {
        for( l = 0; l < n; l++ )
          {
          for( i = m - 1; i >= 0; i-- )
            {
            z = 0;
            for( j = pntr[i] + 1; j < pntr[i + 1]; j++ )
              {
              z += _SpMatVal(c, ldc, indx[j], l) * val[j];
              }
            _SpMatVal(c, ldc, i + pntr[0], l) =
              alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
            }
          }
        }
      }
  }
  
  /*
   * C <- alpha D A^{-1} B + beta C
   * C <- alpha A^{-1} D B + beta C
   */
  void F77NAME(scscsm)
    (const int & transa, const int & m, const int & n,
    const int & unitd, const float * dv, const float & alpha,
    const int descra[], const float * val,
    const int * indx, const int * pntr, const float * b, int & ldb,
    const float & beta, float * c, const int & ldc,
    float * work, const int & lwork)
    {
    if( descra[0] != 0 )
      {
      std::cerr << "Must have general matrix" << "\n";
      exit(1);
      }
  
    if( beta != 0.0 )
      {
      std::cerr << "beta != 0 not implemented" << "\n";
      exit(1);
      }
  
    // To make the compiler happy
    if( work && lwork )
      {
      ;
      }
  
    if( transa == 0 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompCol_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompCol_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      }
    else if( transa == 1 || transa == 2 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompRow_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompRow_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else
        {
        std::cerr << "Bad descra[1]" << "\n";
        exit(1);
        }
      }
    else
      {
      std::cerr << "Bad transa" << "\n";
      exit(1);
      }
    }
  
  void F77NAME(scsrsm)
    (const int & transa, const int & m, const int & n,
    const int & unitd, const float * dv, const float & alpha,
    const int descra[], const float * val,
    const int * indx, const int * pntr, const float * b, int & ldb,
    const float & beta, float * c, const int & ldc,
    float * work, const int & lwork)
    {
    if( descra[0] != 0 )
      {
      std::cerr << "Must have general matrix" << "\n";
      exit(1);
      }
  
    if( beta != 0.0 )
      {
      std::cerr << "beta != 0 not implemented" << "\n";
      exit(1);
      }
  
    // To make the compiler happy
    if( work && lwork )
      {
      ;
      }
  
    if( transa == 0 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompRow_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompRow_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      }
    else if( transa == 1 || transa == 2 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompCol_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompCol_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                                       b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else
        {
        std::cerr << "Bad descra[1]" << "\n";
        exit(1);
        }
      }
    else
      {
      std::cerr << "Bad transa" << "\n";
      exit(1);
      }
    }
  
  void F77NAME(dcscsm)
    (const int & transa, const int & m, const int & n,
    const int & unitd, const double * dv, const double & alpha,
    const int descra[], const double * val,
    const int * indx, const int * pntr, const double * b, int & ldb,
    const double & beta, double * c, const int & ldc,
    double * work, const int & lwork)
    {
    if( descra[0] != 0 )
      {
      std::cerr << "Must have general matrix" << "\n";
      exit(1);
      }
  
    if( beta != 0.0 )
      {
      std::cerr << "beta != 0 not implemented" << "\n";
      exit(1);
      }
  
    // To make the compiler happy
    if( work && lwork )
      {
      ;
      }
  
    if( transa == 0 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompCol_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompCol_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      }
    else if( transa == 1 || transa == 2 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompRow_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompRow_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else
        {
        std::cerr << "Bad descra[1]" << "\n";
        exit(1);
        }
      }
    else
      {
      std::cerr << "Bad transa" << "\n";
      exit(1);
      }
    }
  
  void F77NAME(dcsrsm)
    (const int & transa, const int & m, const int & n,
    const int & unitd, const double * dv, const double & alpha,
    const int descra[], const double * val,
    const int * indx, const int * pntr, const double * b, int & ldb,
    const double & beta, double * c, const int & ldc,
    double * work, const int & lwork)
    {
    if( descra[0] != 0 )
      {
      std::cerr << "Must have general matrix" << "\n";
      exit(1);
      }
  
    if( beta != 0.0 )
      {
      std::cerr << "beta != 0 not implemented" << "\n";
      exit(1);
      }
  
    // To make the compiler happy
    if( work && lwork )
      {
      ;
      }
  
    if( transa == 0 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompRow_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompRow_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompRow_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      }
    else if( transa == 1 || transa == 2 )
      {
      if( descra[1] == 1 )
        {
        if( descra[2] == 0 )
          {
          CompCol_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else if( descra[1] == 2 )
        {
        if( descra[2] == 0 )
          {
          CompCol_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else if( descra[2] == 1 )
          {
          CompCol_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                                        b, ldb, c, ldc);
          }
        else
          {
          std::cerr << "Bad descra[2]" << "\n";
          exit(1);
          }
        }
      else
        {
        std::cerr << "Bad descra[1]" << "\n";
        exit(1);
        }
      }
    else
      {
      std::cerr << "Bad transa" << "\n";
      exit(1);
      }
    }
  