Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

vector.cpp

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                      Linear Algebra Package
00006  *
00007  *              Basic Linear Algebra operations, level 1 & 2
00008  *                   concerning specifically vectors
00009  *
00010  * The present file is concerned with the operations which either
00011  *      - specifically defined for vectors, such as norms
00012  *      - some BLAS 1 & 2 operations that can be implemented more 
00013  *        efficiently than generic operations on n*1 matrices
00014  *
00015  * $Id: vector.cpp,v 1.1 2004/05/21 21:02:52 maxx Exp $
00016  *
00017  ************************************************************************
00018  */
00019 
00020 #include "LinAlg.h"
00021 #include <math.h>
00022 #include "builtin.h"
00023 #include <stdarg.h>
00024 
00025 namespace linalg 
00026 {
00027     using namespace linalg;
00028     
00029 /*
00030  *------------------------------------------------------------------------
00031  *                     Specific vector constructors
00032  */
00033 
00034                         // Make a vector and assign initial values
00035                         // Argument list should contain DOUBLE values
00036                         // to assign to vector elements. The list must
00037                         // be terminated by the string "END"
00038                         // Example: Vector foo(1,3,0.0,1.0,1.5,"END");
00039 Vector::Vector(const int lwb, const int upb, double iv1, ... )
00040   : Matrix(lwb,upb,1,1)
00041 {
00042   va_list args;
00043   va_start(args,iv1);                   // Init 'args' to the beginning of
00044                                         // the variable length list of args
00045   register int i;
00046   (*this)(lwb) = iv1;
00047   for(i=lwb+1; i<=upb; i++)
00048     (*this)(i) = (double)va_arg(args,double);
00049 
00050   assure( strcmp((char *)va_arg(args,char *),"END") == 0,
00051          "Vector: argument list must be terminated by \"END\" ");
00052 }
00053 
00054                                 // Resize the vector for a specified number
00055                                 // of elements, trying to keep intact as many
00056                                 // elements of the old vector as possible.
00057                                 // If the vector is expanded, the new elements
00058                                 // will be zeroes
00059 void Vector::resize_to(const int lwb, const int upb)
00060 {
00061   is_valid();
00062   const int old_nrows = nrows;
00063   assure( (nrows = upb-lwb+1) > 0,
00064          "can't resize vector to a non-positive number of elems" );
00065 
00066   row_lwb = lwb;
00067   if( old_nrows == nrows )
00068     return;                                     // The same number of elems
00069 
00070   nelems = nrows;
00071   assert( !ref_counter.q_engaged() );
00072 
00073                                 // If the vector is to grow, reallocate
00074                                 // and clear the newly added elements
00075   if( nrows > old_nrows )
00076     elements = (REAL *)realloc(elements,nelems*sizeof(REAL)),
00077     memset(elements+old_nrows,0,(nrows-old_nrows)*sizeof(REAL));
00078 
00079                                 // Vector is to shrink a lot (more than
00080                                 // 7/8 of the original size), reallocate
00081   else if( old_nrows - nrows > (old_nrows>>3) )
00082     elements = (REAL *)realloc(elements,nelems*sizeof(REAL));
00083 
00084                                 // If the vector shrinks only a little, don't
00085                                 // bother reallocating
00086   assert( elements != 0 );
00087 }
00088 
00089 /*
00090  *------------------------------------------------------------------------
00091  *              Multiplications specifically defined for vectors
00092  */
00093 
00094                                 // Compute the scalar product
00095 double operator * (const Vector& v1, const Vector& v2)
00096 {
00097   are_compatible(v1,v2);
00098   register REAL * v1p = v1.elements;
00099   register REAL * v2p = v2.elements;
00100   register double sum = 0;
00101 
00102   while( v1p < v1.elements + v1.nelems )
00103     sum += *v1p++ * *v2p++;
00104 
00105   return sum;
00106 }
00107 
00108                                         // "Inplace" multiplication
00109                                         // target = A*target
00110                                         // A needn't be a square one (the
00111                                         // target will be resized to fit)
00112 Vector& Vector::operator *= (const Matrix& A)
00113 {
00114   A.is_valid();
00115   is_valid();
00116 
00117   if( A.ncols != nrows || A.col_lwb != row_lwb )
00118     A.info(), info(),
00119     _error("matrices and vector above cannot be multiplied");
00120   assert( !ref_counter.q_engaged() );
00121 
00122   const int old_nrows = nrows;
00123   const REAL * old_vector = elements;   // Save the old vector elem
00124   row_lwb = A.row_lwb;
00125   assert( (nrows = A.nrows) > 0 );
00126 
00127   nelems = nrows;                       // Allocate new vector elements
00128   assert( (elements = (REAL *)malloc(nelems*sizeof(REAL))) != 0 );
00129 
00130   register REAL * tp = elements;        // Target vector ptr
00131   register REAL * mp = A.elements;      // Matrix row ptr
00132   while( tp < elements + nelems )
00133   {
00134     register double sum = 0;
00135     for( register const REAL * sp = old_vector; sp < old_vector + old_nrows; )
00136       sum += *sp++ * *mp, mp += A.nrows;
00137     *tp++ = sum;
00138     mp -= A.nelems -1;                  // mp points to the beg of the next row
00139   }
00140   assert( mp == A.elements + A.nrows );
00141 
00142   free((REAL *)old_vector);
00143   return *this;
00144 }
00145  
00146 }

Generated on Wed Dec 15 21:20:32 2004 for vuVolume by  doxygen 1.3.9.1