00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00032
00033
00034
00035
00036
00037
00038
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);
00044
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
00055
00056
00057
00058
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;
00069
00070 nelems = nrows;
00071 assert( !ref_counter.q_engaged() );
00072
00073
00074
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
00080
00081 else if( old_nrows - nrows > (old_nrows>>3) )
00082 elements = (REAL *)realloc(elements,nelems*sizeof(REAL));
00083
00084
00085
00086 assert( elements != 0 );
00087 }
00088
00089
00090
00091
00092
00093
00094
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
00109
00110
00111
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;
00124 row_lwb = A.row_lwb;
00125 assert( (nrows = A.nrows) > 0 );
00126
00127 nelems = nrows;
00128 assert( (elements = (REAL *)malloc(nelems*sizeof(REAL))) != 0 );
00129
00130 register REAL * tp = elements;
00131 register REAL * mp = A.elements;
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;
00139 }
00140 assert( mp == A.elements + A.nrows );
00141
00142 free((REAL *)old_vector);
00143 return *this;
00144 }
00145
00146 }