BCLS: Bound Constrained Least Squares

Version 0.1

bccblas.c

Go to the documentation of this file.
00001 /* bccblas.c
00002    $Revision: 231 $ $Date: 2006-04-15 18:47:05 -0700 (Sat, 15 Apr 2006) $
00003 
00004    ----------------------------------------------------------------------
00005    This file is part of BCLS (Bound-Constrained Least Squares).
00006 
00007    Copyright (C) 2006 Michael P. Friedlander, Department of Computer
00008    Science, University of British Columbia, Canada. All rights
00009    reserved. E-mail: <mpf@cs.ubc.ca>.
00010    
00011    BCLS is free software; you can redistribute it and/or modify it
00012    under the terms of the GNU Lesser General Public License as
00013    published by the Free Software Foundation; either version 2.1 of the
00014    License, or (at your option) any later version.
00015    
00016    BCLS is distributed in the hope that it will be useful, but WITHOUT
00017    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00018    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
00019    Public License for more details.
00020    
00021    You should have received a copy of the GNU Lesser General Public
00022    License along with BCLS; if not, write to the Free Software
00023    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00024    USA
00025    ----------------------------------------------------------------------
00026 */
00044 #include "cblas.h"
00045 
00054 void
00055 cblas_daxpy( const int N, const double alpha, const double *X,
00056              const int incX, double *Y, const int incY)
00057 {
00058   int i;
00059 
00060   if (N     <= 0  ) return;
00061   if (alpha == 0.0) return;
00062 
00063   if (incX == 1 && incY == 1) {
00064       const int m = N % 4;
00065 
00066       for (i = 0; i < m; i++)
00067           Y[i] += alpha * X[i];
00068       
00069       for (i = m; i + 3 < N; i += 4) {
00070           Y[i    ] += alpha * X[i    ];
00071           Y[i + 1] += alpha * X[i + 1];
00072           Y[i + 2] += alpha * X[i + 2];
00073           Y[i + 3] += alpha * X[i + 3];
00074       }
00075   } else {
00076       int ix = OFFSET(N, incX);
00077       int iy = OFFSET(N, incY);
00078 
00079       for (i = 0; i < N; i++) {
00080           Y[iy] += alpha * X[ix];
00081           ix    += incX;
00082           iy    += incY;
00083       }
00084   }
00085 }
00086 
00094 void
00095 cblas_dcopy( const int N, const double *X,
00096              const int incX, double *Y, const int incY)
00097 {
00098   int i;
00099   int ix = OFFSET(N, incX);
00100   int iy = OFFSET(N, incY);
00101 
00102   for (i = 0; i < N; i++) {
00103       Y[iy]  = X[ix];
00104       ix    += incX;
00105       iy    += incY;
00106   }
00107 }
00108 
00119 double
00120 cblas_ddot( const int N, const double *X,
00121             const int incX, const double *Y, const int incY)
00122 {
00123   double r  = 0.0;
00124   int    i;
00125   int    ix = OFFSET(N, incX);
00126   int    iy = OFFSET(N, incY);
00127 
00128   for (i = 0; i < N; i++) {
00129       r  += X[ix] * Y[iy];
00130       ix += incX;
00131       iy += incY;
00132   }
00133   
00134   return r;
00135 }
00136 
00144 double
00145 cblas_dnrm2( const int N, const double *X, const int incX) 
00146 {
00147   double
00148       scale = 0.0,
00149       ssq   = 1.0;
00150   int
00151       i,
00152       ix    = 0;
00153 
00154   if (N <= 0 || incX <= 0) return 0;
00155   else if (N == 1)         return fabs(X[0]);
00156 
00157   for (i = 0; i < N; i++) {
00158       const double x = X[ix];
00159 
00160       if (x != 0.0) {
00161           const double ax = fabs(x);
00162 
00163           if (scale < ax) {
00164               ssq   = 1.0 + ssq * (scale / ax) * (scale / ax);
00165               scale = ax;
00166           } else {
00167               ssq += (ax / scale) * (ax / scale);
00168           }
00169       }
00170 
00171       ix += incX;
00172   }
00173   
00174   return scale * sqrt(ssq);
00175 }
00176 
00183 void
00184 cblas_dscal(const int N, const double alpha, double *X, const int incX)
00185 {
00186     int i, ix;
00187 
00188     if (incX <= 0) return;
00189 
00190     ix = OFFSET(N, incX);
00191     
00192     for (i = 0; i < N; i++) {
00193         X[ix] *= alpha;
00194         ix    += incX;
00195     }
00196 }

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