BCLS: Bound Constrained Least Squares

Version 0.1

bccblas.c

Go to the documentation of this file.

/* bccblas.c
   $Revision: 231 $ $Date: 2006-04-15 18:47:05 -0700 (Sat, 15 Apr 2006) $

   ----------------------------------------------------------------------
   This file is part of BCLS (Bound-Constrained Least Squares).

   Copyright (C) 2006 Michael P. Friedlander, Department of Computer
   Science, University of British Columbia, Canada. All rights
   reserved. E-mail: <mpf@cs.ubc.ca>.
   
   BCLS is free software; you can redistribute it and/or modify it
   under the terms of the GNU Lesser General Public License as
   published by the Free Software Foundation; either version 2.1 of the
   License, or (at your option) any later version.
   
   BCLS is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
   Public License for more details.
   
   You should have received a copy of the GNU Lesser General Public
   License along with BCLS; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
   USA
   ----------------------------------------------------------------------
*/
#include "cblas.h"

void
cblas_daxpy( const int N, const double alpha, const double *X,
             const int incX, double *Y, const int incY)
{
  int i;

  if (N     <= 0  ) return;
  if (alpha == 0.0) return;

  if (incX == 1 && incY == 1) {
      const int m = N % 4;

      for (i = 0; i < m; i++)
          Y[i] += alpha * X[i];
      
      for (i = m; i + 3 < N; i += 4) {
          Y[i    ] += alpha * X[i    ];
          Y[i + 1] += alpha * X[i + 1];
          Y[i + 2] += alpha * X[i + 2];
          Y[i + 3] += alpha * X[i + 3];
      }
  } else {
      int ix = OFFSET(N, incX);
      int iy = OFFSET(N, incY);

      for (i = 0; i < N; i++) {
          Y[iy] += alpha * X[ix];
          ix    += incX;
          iy    += incY;
      }
  }
}

void
cblas_dcopy( const int N, const double *X,
             const int incX, double *Y, const int incY)
{
  int i;
  int ix = OFFSET(N, incX);
  int iy = OFFSET(N, incY);

  for (i = 0; i < N; i++) {
      Y[iy]  = X[ix];
      ix    += incX;
      iy    += incY;
  }
}

double
cblas_ddot( const int N, const double *X,
            const int incX, const double *Y, const int incY)
{
  double r  = 0.0;
  int    i;
  int    ix = OFFSET(N, incX);
  int    iy = OFFSET(N, incY);

  for (i = 0; i < N; i++) {
      r  += X[ix] * Y[iy];
      ix += incX;
      iy += incY;
  }
  
  return r;
}

double
cblas_dnrm2( const int N, const double *X, const int incX) 
{
  double
      scale = 0.0,
      ssq   = 1.0;
  int
      i,
      ix    = 0;

  if (N <= 0 || incX <= 0) return 0;
  else if (N == 1)         return fabs(X[0]);

  for (i = 0; i < N; i++) {
      const double x = X[ix];

      if (x != 0.0) {
          const double ax = fabs(x);

          if (scale < ax) {
              ssq   = 1.0 + ssq * (scale / ax) * (scale / ax);
              scale = ax;
          } else {
              ssq += (ax / scale) * (ax / scale);
          }
      }

      ix += incX;
  }
  
  return scale * sqrt(ssq);
}

void
cblas_dscal(const int N, const double alpha, double *X, const int incX)
{
    int i, ix;

    if (incX <= 0) return;

    ix = OFFSET(N, incX);
    
    for (i = 0; i < N; i++) {
        X[ix] *= alpha;
        ix    += incX;
    }
}

Generated on Sun Mar 4 22:50:03 2007 by Doxygen 1.5.1