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 "EpetraVectorAdapter.hpp" 00017 00018 EpetraVectorAdapter::EpetraVectorAdapter(Epetra_MultiVector & v_): 00019 v(&v_), 00020 ownEpetraVector(false) 00021 { 00022 assert( v->NumVectors() == 1); 00023 v->ExtractView(&vals, &localSize); 00024 } 00025 00026 EpetraVectorAdapter::EpetraVectorAdapter(Epetra_MultiVector * v_, bool owned): 00027 v(v_), 00028 ownEpetraVector(owned) 00029 { 00030 assert( v->NumVectors() == 1); 00031 v->ExtractView(&vals, &localSize); 00032 } 00033 00034 EpetraVectorAdapter::~EpetraVectorAdapter() 00035 { 00036 if(ownEpetraVector) 00037 delete v; 00038 } 00039 00041 EpetraVectorAdapter & EpetraVectorAdapter::operator=(const double & val) 00042 { 00043 std::fill(vals, vals+localSize, val); 00044 return *this; 00045 } 00046 00048 EpetraVectorAdapter & EpetraVectorAdapter::operator=(const EpetraVectorAdapter & RHS) 00049 { 00050 assert( localSize == RHS.localSize); 00051 std::copy(RHS.vals, RHS.vals+ RHS.localSize, vals); 00052 return *this; 00053 } 00054 00056 void EpetraVectorAdapter::Scale(const double & val) 00057 { 00058 for( double * it(vals); it != vals + localSize; ++it) 00059 (*it) *= val; 00060 } 00061 00062 void EpetraVectorAdapter::Randomize(int seed) 00063 { 00064 v->SetSeed( seed ); 00065 v->Random(); 00066 double norm2( InnerProduct(*this, *this) ); 00067 Scale(1./sqrt(norm2)); 00068 00069 } 00070 00071 00073 EpetraVectorAdapter * EpetraVectorAdapter::Clone() 00074 { 00075 Epetra_MultiVector * cv(new Epetra_MultiVector(v->Map(), v->NumVectors()) ); 00076 return new EpetraVectorAdapter(cv); 00077 } 00078 00079 Epetra_MultiVector & EpetraVectorAdapter::EpetraVector() 00080 { 00081 return *v; 00082 } 00083 00084 const Epetra_MultiVector & EpetraVectorAdapter::EpetraVector() const 00085 { 00086 return *v; 00087 } 00088 00089 void EpetraVectorAdapter::Print(std::ostream & os) const 00090 { 00091 double * end(vals+localSize); 00092 for(double * it(vals); it != end; ++it) 00093 os << *it << "\t"; 00094 os << "\n"; 00095 } 00096 00098 void add(const EpetraVectorAdapter & v1, const double & c2, const EpetraVectorAdapter & v2, EpetraVectorAdapter & result) 00099 { 00100 int size(result.localSize); 00101 double * rr(result.vals); 00102 double * vv1(v1.vals); 00103 double * vv2(v2.vals); 00104 00105 assert( v1.localSize == v2.localSize ); 00106 assert( size == v1.localSize ); 00107 00108 double * end(rr + size); 00109 00110 for( ; rr != end; ++rr, ++vv1, ++vv2) 00111 (*rr) = (*vv1) + c2*(*vv2); 00112 00113 } 00115 void add(const double & c1, const EpetraVectorAdapter & v1, const double & c2, const EpetraVectorAdapter & v2, EpetraVectorAdapter & result) 00116 { 00117 int size(result.localSize); 00118 double * rr(result.vals); 00119 double * vv1(v1.vals); 00120 double * vv2(v2.vals); 00121 00122 assert( v1.localSize == v2.localSize ); 00123 assert( size == v1.localSize ); 00124 00125 double * end(rr + size); 00126 00127 for( ; rr != end; ++rr, ++vv1, ++vv2) 00128 (*rr) = c1*(*vv1) + c2*(*vv2); 00129 00130 00131 } 00133 void add(const double & alpha, const EpetraVectorAdapter & v1, const EpetraVectorAdapter & v2, EpetraVectorAdapter & result) 00134 { 00135 int size(result.localSize); 00136 double * rr(result.vals); 00137 double * vv1(v1.vals); 00138 double * vv2(v2.vals); 00139 00140 assert( v1.localSize == v2.localSize ); 00141 assert( size == v1.localSize ); 00142 00143 double * end(rr + size); 00144 00145 for( ; rr != end; ++rr, ++vv1, ++vv2) 00146 (*rr) = alpha*(*vv1 + *vv2); 00147 00148 } 00150 void add(const EpetraVectorAdapter & v1, const EpetraVectorAdapter & v2, const EpetraVectorAdapter & v3, EpetraVectorAdapter & result) 00151 { 00152 int size(result.localSize); 00153 double * rr(result.vals); 00154 double * vv1(v1.vals); 00155 double * vv2(v2.vals); 00156 double * vv3(v3.vals); 00157 00158 assert( v1.localSize == v2.localSize ); 00159 assert( v1.localSize == v3.localSize ); 00160 assert( size == v1.localSize ); 00161 00162 double * end(rr + size); 00163 00164 for( ; rr != end; ++rr, ++vv1, ++vv2, ++vv3) 00165 (*rr) = (*vv1) + (*vv2) + (*vv3); 00166 00167 00168 } 00170 void subtract(const EpetraVectorAdapter & v1, const EpetraVectorAdapter & v2, EpetraVectorAdapter & result) 00171 { 00172 int size(result.localSize); 00173 double * rr(result.vals); 00174 double * vv1(v1.vals); 00175 double * vv2(v2.vals); 00176 00177 assert( v1.localSize == v2.localSize ); 00178 assert( size == v1.localSize ); 00179 00180 double * end(rr + size); 00181 00182 for( ; rr != end; ++rr, ++vv1, ++vv2) 00183 (*rr) = (*vv1 - *vv2); 00184 00185 } 00187 double InnerProduct(const EpetraVectorAdapter & v1, const EpetraVectorAdapter & v2) 00188 { 00189 assert(v1.v->NumVectors() == 1); 00190 assert(v2.v->NumVectors() == 1); 00191 00192 double result; 00193 v1.v->Dot(*v2.v, &result); 00194 00195 return result; 00196 }