blas_interface.cc

MTL C++ Numeric Library

    swap(x, strided(y,  iy));
  else if (iy == 1)
    swap(strided(x, ix), y);
  else
    swap(strided(x), strided(y,  iy));
}

//------------------------------------------------------------
//  NRM2 for floats (the 2 norm)
//------------------------------------------------------------
float
MTL_FCALL(snrm2)(int *n, float *sx, int *incx)
{
  int N = *n;
  int ix = *incx;

  float sw = 0.0;

  if (N <= 0) return sw;

  svec x(sx, N * ix);
  
  if (ix == 1)
    sw = two_norm(x);
  else
    sw = two_norm(strided(x, ix));

  return sw;

}

//------------------------------------------------------------
//  NRM2 for doubles (the 2 norm)
//------------------------------------------------------------
double
MTL_FCALL(dnrm2)(int *n, double *dx, int *incx)
{
  int N = *n;
  int ix = *incx;

  double dw = 0.0;

  if (N <= 0) return dw;

  dvec x(dx, N * ix);
  
  if (ix == 1)
    dw = two_norm(x);
  else
    dw = two_norm(strided(x, ix));

  return dw;
  
}

//------------------------------------------------------------
//  SUM of abs. values for floats
//------------------------------------------------------------

float
MTL_FCALL(sasum)(int *n, float *sx, int *incx)
{
  int N = *n;
  int ix = *incx;
  float sw = 0.0;

  if (N <= 0) return sw;

  svec x(sx, N * ix);

  if (ix == 1)
    sw = one_norm(x);
  else
    sw = one_norm(strided(x, ix));    

  return sw;
}

//------------------------------------------------------------
//  SUM of abs. values for doubles
//------------------------------------------------------------

double
MTL_FCALL(dasum)(int *n, double *dx, int *incx)
{
  int N = *n;
  int ix = *incx;

  double dw = 0.0;

  if (N <= 0) return dw;

  dvec x(dx, N * ix);

  if (ix == 1)
    dw = one_norm(x);
  else
    dw = one_norm(strided(x, ix));

  return dw;
}

//------------------------------------------------------------
//  Scale for floats
p//------------------------------------------------------------
void
MTL_FCALL(sscal)(int *n, float *sa, float *sx, int *incx)
{
  int N = *n;
  float a = *sa;
  int ix = *incx;

  if (N <= 0 || a == 1) return;
  
  svec x(sx, N * ix);

  if (ix == 1)
    scaled(x, a);
  else {
    scaled(strided(x, ix), a);
  }

}

//------------------------------------------------------------
//  Scale for doubles
//------------------------------------------------------------
void
MTL_FCALL(dscal)(int *n, double *da, double *dx, int *incx)
{
  int N = *n;
  double a = *da;
  int ix = *incx;

  if (N <= 0 || a == 1) return;
  
  dvec x(dx, N * ix);

  if (ix == 1)
    scaled(x, a);
  else {
    strided<dvec> sx(x, ix);
    scaled(sx, a);
  }

}
  
//------------------------------------------------------------
//  Largest component for floats (abs. value
//------------------------------------------------------------
int
MTL_FCALL(isamax)(int *n, float *sx, int *incx)
{
  int N = *n;
  int ix = *incx;
  int imax = 0;

  if (N <= 0) return imax;

  svec x(sx, N * ix);

  if (ix == 1)
    imax = max_index(x);
  else
    imax = max_index(strided(x, ix));

  return imax + 1;
}
 
//------------------------------------------------------------
//  Largest component for doubles (abs. value
//------------------------------------------------------------
int
MTL_FCALL(idamax)(int *n, double *dx, int *incx)
{
  int N = *n;
  int ix = *incx;
  int imax = 0;

  if (N <= 0) return imax;

  dvec x(dx, N * ix);

  if (ix == 1)
    imax = max_index(x);
  else
    imax = max_index(strided(x, ix));

  return imax + 1;
}

//**************************************************
// This stuff is the givens transformatioin
// stuff that is not implemented yet.  That
// is why it is #ifdef'ed out.
//**************************************************

#ifdef GIVENS

//****************************************
//  Modified rotations are omitted
//****************************************

//------------------------------------------------------------
//  Construct a Givens Plane Rotation for floats
//------------------------------------------------------------
void
srotg(float *sa, float *sb, float *sc, float *ss)
{
  float a = *sa;
  float b = *sb;

  givens_rotation<float> g(a, b);

  *sa = g.a();
  *sb = g.b();
  *sc = g.c();
  *ss = g.s();

#if 0
  float c = *sc;
  float s = *ss;
  float r = 0.0;
  float s = 0.0;
  int sigma = 1;  // should always equal 1 or -1;
  if (std::abs(a) > std::abs(b))
    sigma = std::abs(a)/a;
  else if (std::abs(b) >= std::abs(a))
    sigma = std::abs(b)/b;
  else 
    cout << "Error in SROTG caculating sigma" << endl;

  r = sigma*(sqrt(a*a) + (b*b));
  
  if (r != 0){
    c = a/r;
    s = b/r;
  }
  else {
    c = 1;
    s = 0;
  }

  float z = 0.0;

  if (c != 0) {
    if (std::abs(a) > std::abs(b))
      z = s;
    else
      z = 1/c;
  } else
    z = 1;

  *sa = r;
  *sb = z;
  *sc = c;
  *ss = s;
#endif
}


//------------------------------------------------------------
//  Construct a Givens Plane Rotation for doubles
//------------------------------------------------------------
void
MTL_FCALL(drotg)(double *da, double *db, double *dc, double *ds)
{
  double a = *da;
  double b = *db;

  givens_rotation<double> g(a, b);

  *da = g.a();
  *db = g.b();
  *dc = g.c();
  *ds = g.s();
}

//------------------------------------------------------------
//  Applying a Givens Plane Rotation for floats
//------------------------------------------------------------
void
MTL_FCALL(srot)(int *n, float *sx, int *incx, float *sy, 
      int *incy, float *sc, float *ss)
{
  int N = *n;
  int ix = *incx;
  int iy = *incy;
  float c = *sc;
  float s = *ss;

  if (N <= 0) return;

  s_vec x(sx, N), y(sy, N);
  givens_rotation<float> g;
  g.set_cs(c, s);
  g.apply(strided(x, ix), strided(y, iy));
}

//------------------------------------------------------------
//  Applying a Givens Plane Rotation for doubles
//------------------------------------------------------------
void
MTL_FCALL(drot)(int *n, double *dx, int *incx, double *dy, 
		int *incy, double *dc, double *ds)
{
  int N = *n;
  int ix = *incx;
  int iy = *incy;
  double c = *dc;
  double s = *ds;

  if (N <= 0) return;

  d_vec x(dx, N), y(dy, N);
  givens_rotation<double> g;
  g.set_cs(c, s);
  g.apply(strided(x, ix), strided(y, iy));
}



void MTL_FCALL(dgemv)(char* transp, int* m, int* n, double* alphap,
		      double* a, int* ldap, double* xp, int* incxp,
		      double* betap, double* yp, int* incyp)
{
  char trans = *transp;
  int M = *m;
  int N = *n;
  double alpha = *alphap;
  double beta = *betap;
  int lda = *ldap;
  int incx = *incxp;
  int incy = *incyp;

  dvec x(xp, N), y(yp, M);
  col_matrix_d A(a, M, N, lda);

  if (trans == 'N' || trans == 'n')
     mult(A, strided(scaled(x,alpha),incx), 
	  strided(scaled(s, beta), iy));
  else
    mult(trans(A), strided(scaled(x,alpha),incx), 
	 strided(scaled(y, beta), iy));
}



void MTL_FCALL(dger)(int*, int*, double*, double*, int*, double*,
		     int*, double*, int*)
{

}


#endif