//// Copyright 1997, 1998, 1999 University of Notre Dame.// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee//// This file is part of the Matrix Template Library//// You should have received a copy of the License Agreement for the// Matrix Template Library along with the software;  see the// file LICENSE.  If not, contact Office of Research, University of Notre// Dame, Notre Dame, IN  46556.//// Permission to modify the code and to distribute modified code is// granted, provided the text of this NOTICE is retained, a notice that// the code was modified is included with the above COPYRIGHT NOTICE and// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE// file is distributed with the modified code.//// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.// By way of example, but not limitation, Licensor MAKES NO// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS// OR OTHER RIGHTS.////------------------------------------------------------------// Basic Linear Algebra Subprograms for C/C++// Version 1.0// Matthew E. Gaston// May 6, 1998// ------------------------------------------------------------#include "blas_interface.h"#include "mtl/linalg_vec.h"#include "mtl/matrix.h"#include "mtl/mtl.h"#include <math.h>// TO DO: check whether int args should be const// and then change the prototypes and func defsusing mtl;typedef external_vec<float> svec;typedef external_vec<double> dvec;typedef matrix<double, rectangle<>, dense<>, column_major>::type col_matrix_d;typedef matrix<double, rectangle<>, dense<>, column_major>::type col_matrix_s;#define MTL_FCALL(X) X##_using namespace mtl;//------------------------------------------------------------//  Dot product of floats return a float//------------------------------------------------------------floatMTL_FCALL(sdot)(int *n, float *sx, int *incx, float *sy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;    float sw = 0.0;  if (N <= 0) return sw;  svec x(sx, N * ix);  svec y(sy, N * iy);  if (ix == 1 && iy == 1)    sw = dot(x, y, sw);  else if (ix == 1)    sw = dot(x, strided(y, iy), sw);  else if (iy == 1)    sw = dot(strided(x, ix), y, sw);  else    sw = dot(strided(x, ix), strided(y, iy), sw);  return sw;}    //------------------------------------------------------------//  Dot product of floats return a double//------------------------------------------------------------doubleMTL_FCALL(dsdot)(int *n, float *sx, int *incx, float *sy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;    double dw = 0.0;  if (N <= 0) return dw;  svec x(sx, N * ix);  svec y(sy, N * iy);  if (ix == 1 && iy == 1)    dw = dot(x, y, dw);  else if (ix == 1)    dw = dot(x, strided(y, iy), dw);  else if (iy == 1)    dw = dot(strided(x, ix), y, dw);  else    dw = dot(strided(x, ix), strided(y, iy), dw);  return dw;}//------------------------------------------------------------//  Dot product of floats plus a float return a float//------------------------------------------------------------floatMTL_FCALL(sdsdot)(int *n, float *sb, float *sx, int *incx, float *sy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;  float b = *sb;    float sw = 0.0;  if (N <= 0) return sw;  svec x(sx, N * ix);  svec y(sy, N * iy);  if (ix == 1 && iy == 1)    sw = dot(x, y, sw);  else if (ix == 1)    sw = dot(x, strided(y, iy), sw);  else if (iy == 1)    sw = dot(strided(x, ix), y, sw);  else    sw = dot(strided(x, ix), strided(y, iy), sw);  sw += b;  return sw;}//------------------------------------------------------------//  Dot product of doubles return a double//------------------------------------------------------------doubleMTL_FCALL(ddot)(int *n, double *dx, int *incx, double *dy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;    double dw = 0.0;  if (N <= 0) return dw;  dvec x(dx, N * ix);  dvec y(dy, N * iy);  if (ix == 1 && iy == 1)    dw = dot(x, y, dw);  else if (ix == 1)    dw = dot(x, strided(y, iy), dw);  else if (iy == 1)    dw = dot(strided(x, ix), y, dw);  else    dw = dot(strided(x, ix), strided(y, iy), dw);  return dw;}//------------------------------------------------------------//  AXPY for floats//------------------------------------------------------------voidMTL_FCALL(saxpy)(int *n, float *sa, float *sx, int *incx, float *sy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;  float a = *sa;  if (a == 0 || N <= 0) return;  svec x(sx, N * ix);  svec y(sy, N * iy);  // no scale to x -- i.e. a = 1;  if (a == 1) {    if (ix == 1 && iy == 1)      add(x, y, y);    else if (iy == 1)      add(strided(x, ix), y, y);    else if (ix == 1)      add(x, strided(y, iy), strided(y, iy));    else      add(strided(x, ix), strided(y, iy), strided(y,iy));  }  // must scale x  else {    if (ix == 1 && iy == 1)      add(scaled(x, a), y, y);    else if (iy == 1)      add(scaled(strided(x, ix), a), y, y);    else if (ix == 1)      add(scaled(x, a), strided(y, iy), strided(y, iy));    else      add(scaled(strided(x, ix), a), strided(y, iy), strided(y, iy));  }}//------------------------------------------------------------//  AXPY for doubles//------------------------------------------------------------voidMTL_FCALL(daxpy)(int *n, double *da, double *dx, int *incx,       double *dy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;  double a = *da;  if (a == 0 || N <= 0) return;  dvec x(dx, N * ix);  dvec y(dy, N * iy);  // no scale to x -- i.e. a = 1;  if (a == 1) {    if (ix == 1 && iy == 1)      add(x, y, y);    else if (iy == 1)      add(strided(x, ix), y, y);    else if (ix == 1)      add(x, strided(y, iy), strided(y,iy));     else       add(strided(x, ix), strided(y,iy), strided(y,iy));  }  // must scale x  else {    if (ix == 1 && iy == 1)      add(scaled(x, a), y, y);    else if (iy == 1)      add(scaled(strided(x, ix), a), y, y);    else if (ix == 1)      add(scaled(x, a), strided(y,iy), strided(y,iy));    else      add(scaled(strided(x, ix), a), strided(y,  iy), strided(y,  iy));  }}//------------------------------------------------------------//  COPY for floats//------------------------------------------------------------voidMTL_FCALL(scopy)(int *n, float *sx, int *incx, float *sy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;  if (N <= 0) return;  svec x(sx, N * ix);  svec y(sy, N * iy);  if (ix == 1 && iy == 1)    copy(x, y);  else if (ix == 1)    copy(x, strided(y, iy));  else if (iy == 1)    copy(strided(x, ix), y);  else    copy(strided(x, ix), strided(y, iy));}//------------------------------------------------------------//  COPY for doubles//------------------------------------------------------------voidMTL_FCALL(dcopy)(int *n, double *dx, int *incx, double *dy, int *incy){  int N = *n;  int ix = *incx;  int iy = *incy;  if (N <= 0) return;  dvec x(dx, N * ix);  dvec y(dy, N * iy);  if (ix == 1 && iy == 1)    copy(x, y);  else if (ix == 1)    copy(x, strided(y,  iy));  else if (iy == 1)    copy(strided(x, ix), y);  else    copy(strided(x, ix), strided(y,  iy));}//------------------------------------------------------------//  SWAP for floats//------------------------------------------------------------voidMTL_FCALL(sswap)(int *n, float *sx, int *incx, float *sy, int *incy){  int N = *n;   int ix = *incx;  int iy = *incy;  if (N <= 0) return;  svec x(sx, N * ix);  svec y(sy, N * iy);  if (ix == 1 && iy == 1)    swap(x, y);  else if (ix == 1)    swap(x, strided(y,  iy));  else if (iy == 1)    swap(strided(x, ix), y);   else    swap(strided(x, ix), strided(y,  iy));}//------------------------------------------------------------//  SWAP for doubles//------------------------------------------------------------voidMTL_FCALL(dswap)(int *n, double *dx, int *incx, double *dy, int *incy){  int N = *n;   int ix = *incx;  int iy = *incy;  if (N <= 0) return;  dvec x(dx, N * ix);  dvec y(dy, N * iy);  if (ix == 1 && iy == 1)    swap(x, y);  else if (ix == 1)