tminres
|
00001 // tminres is free software; you can redistribute it and/or modify it under the 00002 // terms of the GNU Lesser General Public License (as published by the Free 00003 // Software Foundation) version 2.1 dated February 1999. 00004 // 00005 // Authors: 00006 // - Umberto Villa, Emory University - uvilla@emory.edu 00007 // - Michael Saunders, Stanford University 00008 // - Santiago Akle, Stanford University 00009 00016 #include "SimpleVector.hpp" 00017 #include <cassert> 00018 #include <algorithm> 00019 #include <cmath> 00020 #include <cstdlib> 00021 00022 SimpleVector::SimpleVector(int size_): 00023 size(size_) 00024 { 00025 assert(size > 0); 00026 vals = new double[size]; 00027 } 00028 00029 SimpleVector::~SimpleVector() 00030 { 00031 delete[] vals; 00032 } 00033 00034 SimpleVector & SimpleVector::operator=(const double & val) 00035 { 00036 std::fill(vals, vals+size, val); 00037 return *this; 00038 } 00039 // Set the entry of the Vector equal to the entries in RHS 00040 SimpleVector & SimpleVector::operator=(const SimpleVector & RHS) 00041 { 00042 std::copy(RHS.vals, RHS.vals + RHS.size, vals); 00043 return *this; 00044 } 00045 00046 // multiply THIS by a scalar value 00047 void SimpleVector::Scale(const double & val) 00048 { 00049 for( double * it(vals); it != vals + size; ++it) 00050 (*it) *= val; 00051 } 00052 // Create a new vector with the same structure of THIS. Values are not initialized. 00053 SimpleVector * SimpleVector::Clone() 00054 { 00055 return new SimpleVector(size); 00056 } 00057 00058 double & SimpleVector::operator[](const int i) 00059 { 00060 assert( i < size); 00061 return vals[i]; 00062 } 00063 const double & SimpleVector::operator[](const int i) const 00064 { 00065 assert( i < size ); 00066 return vals[i]; 00067 } 00068 00069 const double SimpleVector::at(const int i) const 00070 { 00071 00072 if (i<0 || i > size-1) 00073 return 0.0; 00074 00075 return vals[i]; 00076 } 00077 00078 void SimpleVector::Randomize(int seed) 00079 { 00080 srand(seed); 00081 for( double * it(vals); it != vals + size; ++it) 00082 (*it) = 2.*static_cast<double>(rand())/static_cast<double>(RAND_MAX) - 1.; 00083 00084 double norm2( InnerProduct(*this, *this) ); 00085 Scale(1./sqrt(norm2)); 00086 00087 } 00088 00089 void SimpleVector::Print(std::ostream & os) 00090 { 00091 for(double *it(vals); it != vals+size; ++it) 00092 os << *it << "\t "; 00093 00094 os << "\n"; 00095 } 00096 00097 // result = v1 + c2*v2 00098 void add(const SimpleVector & v1, const double & c2, const SimpleVector & v2, SimpleVector & result) 00099 { 00100 int size(result.size); 00101 double * rr(result.vals); 00102 double * vv1(v1.vals); 00103 double * vv2(v2.vals); 00104 00105 assert( v1.size == v2.size ); 00106 assert( size == v1.size ); 00107 00108 double * end(rr + size); 00109 00110 for( ; rr != end; ++rr, ++vv1, ++vv2) 00111 (*rr) = (*vv1) + c2*(*vv2); 00112 } 00113 // result = c1*v1 + c2*v2 00114 void add(const double & c1, const SimpleVector & v1, const double & c2, const SimpleVector & v2, SimpleVector & result) 00115 { 00116 int size(result.size); 00117 double * rr(result.vals); 00118 double * vv1(v1.vals); 00119 double * vv2(v2.vals); 00120 00121 assert( v1.size == v2.size ); 00122 assert( size == v1.size ); 00123 00124 double * end(rr + size); 00125 00126 for( ; rr != end; ++rr, ++vv1, ++vv2) 00127 (*rr) = c1*(*vv1) + c2*(*vv2); 00128 } 00129 00130 // result = alpha(v1 + v2) 00131 void add(const double & alpha, const SimpleVector & v1, const SimpleVector & v2, SimpleVector & result) 00132 { 00133 int size(result.size); 00134 double * rr(result.vals); 00135 double * vv1(v1.vals); 00136 double * vv2(v2.vals); 00137 00138 assert( v1.size == v2.size ); 00139 assert( size == v1.size ); 00140 00141 double * end(rr + size); 00142 00143 for( ; rr != end; ++rr, ++vv1, ++vv2) 00144 (*rr) = alpha*(*vv1 +*vv2); 00145 00146 } 00147 00148 // result = v1 + v2 + v3 00149 void add(const SimpleVector & v1, const SimpleVector & v2, const SimpleVector & v3, SimpleVector & result) 00150 { 00151 int size(result.size); 00152 double * rr(result.vals); 00153 double * vv1(v1.vals); 00154 double * vv2(v2.vals); 00155 double * vv3(v3.vals); 00156 00157 assert( v1.size == v2.size ); 00158 assert( v2.size == v3.size ); 00159 assert( size == v1.size ); 00160 00161 double * end(rr + size); 00162 00163 for( ; rr != end; ++rr, ++vv1, ++vv2, ++vv3) 00164 (*rr) = (*vv1) + (*vv2) + (*vv3); 00165 00166 } 00167 // result = v1 - v2 00168 void subtract(const SimpleVector & v1, const SimpleVector & v2, SimpleVector & result) 00169 { 00170 int size(result.size); 00171 double * rr(result.vals); 00172 double * vv1(v1.vals); 00173 double * vv2(v2.vals); 00174 00175 assert( v1.size == v2.size ); 00176 assert( size == v1.size ); 00177 00178 double * end(rr + size); 00179 00180 for( ; rr != end; ++rr, ++vv1, ++vv2) 00181 (*rr) = *vv1 - *vv2; 00182 00183 00184 } 00185 00186 // return the inner product of v1 and v2 00187 double InnerProduct(const SimpleVector & v1, const SimpleVector & v2) 00188 { 00189 double result(0); 00190 int size(v1.size); 00191 double * vv1(v1.vals); 00192 double * vv2(v2.vals); 00193 00194 assert( v1.size == v2.size ); 00195 00196 double * end(vv1 + size); 00197 00198 for( ; vv1 != end; ++vv1, ++vv2) 00199 result += (*vv1) * (*vv2); 00200 00201 return result; 00202 00203 } 00204