41 template<
class Treal,
class Telement>
51 template<
class Treal,
class Telement = Treal>
90 (*this) *= (1.0 / this->
eucl());
94 inline Treal
eucl()
const {
104 static void axpy(Treal
const & alpha,
114 template<
typename TmatrixElement>
115 static void gemv(
bool const tA, Treal
const alpha,
124 template<
typename TmatrixElement>
125 static void symv(
char const uplo, Treal
const alpha,
134 template<
typename TmatrixElement>
135 static void trmv(
char const uplo,
const bool tA,
141 void assign_from_full(Treal
const *
const fullvector,
const int totn);
143 void fullvector(Treal *
const full,
const int totn)
const;
157 template<
class Treal,
class Telement>
160 addFromFull(fullVector);
163 template<
class Treal,
class Telement>
168 for (
int ind = 0; ind < this->n(); ind++)
169 (*
this)(ind).addFromFull(fullVector);
172 template<
class Treal,
class Telement>
175 if (this->is_zero()) {
176 fullVec.resize(this->rows.getNTotalScalars());
177 for (
int row = 0; row < this->nScalars(); ++row )
178 fullVec[this->rows.getOffset()+row] = 0;
181 for (
int ind = 0; ind < this->n(); ind++)
182 (*
this)(ind).fullVector(fullVec);
186 template<
class Treal,
class Telement>
188 delete[] this->elements;
192 template<
class Treal,
class Telement>
197 if (this->is_zero()) {
198 char * tmp = (
char*)&ZERO;
199 file.write(tmp,
sizeof(
int));
202 char * tmp = (
char*)&ONE;
203 file.write(tmp,
sizeof(
int));
204 for (
int i = 0; i < this->n(); i++)
205 this->elements[i].writeToFile(file);
208 template<
class Treal,
class Telement>
213 char tmp[
sizeof(int)];
214 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
222 for (
int i = 0; i < this->n(); i++)
223 this->elements[i].readFromFile(file);
226 throw Failure(
"Vector<Treal, Telement>::"
227 "readFromFile(std::ifstream & file):"
228 "File corruption int value not 0 or 1");
232 template<
class Treal,
class Telement>
238 throw Failure(
"Vector::operator=(int k) only "
239 "implemented for k = 0");
243 template<
class Treal,
class Telement>
247 for (
int ind = 0; ind < this->n(); ind++)
248 (*
this)(ind).random();
253 template<
class Treal,
class Telement>
256 if (!this->is_zero() && alpha != 1) {
257 for (
int ind = 0; ind < this->n(); ind++)
258 (*
this)(ind) *= alpha;
263 template<
class Treal,
class Telement>
267 assert(x.
n() == y.
n());
270 Treal dotProduct = 0;
271 for (
int ind = 0; ind < x.
n(); ind++)
277 template<
class Treal,
class Telement>
282 assert(x.
n() == y.
n());
288 for (
int ind = 0; ind < x.
n(); ind++)
298 template<
class Treal,
class Telement>
299 template<
typename TmatrixElement>
301 gemv(
bool const tA, Treal
const alpha,
310 if ((A.is_zero() || x.
is_zero() || alpha == 0) &&
314 Treal beta_tmp = beta;
319 if (A.is_zero() || x.
is_zero() || alpha == 0)
325 throw Failure(
"Vector<Treal, Telement>::"
326 "gemv(bool const, Treal const, "
327 "const Matrix<Treal, Telement>&, "
328 "const Vector<Treal, Telement>&, "
329 "Treal const, const Vector<Treal, "
331 "Incorrect dimensions for matrix-vector product");
334 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
336 for (
int row = 0; row < A.
nrows(); row++) {
339 for (
int col = 1; col < A.
ncols(); col++)
348 throw Failure(
"Vector<Treal, Telement>::"
349 "gemv(bool const, Treal const, "
350 "const Matrix<Treal, Telement>&, "
351 "const Vector<Treal, Telement>&, "
352 "Treal const, const Vector<Treal, "
354 "Incorrect dimensions for matrix-vector product");
357 #pragma omp parallel for if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared) schedule(dynamic)
359 for (
int col = 0; col < A.
ncols(); col++) {
362 for (
int row = 1; row < A.
nrows(); row++)
376 template<
class Treal,
class Telement>
377 template<
typename TmatrixElement>
379 symv(
char const uplo, Treal
const alpha,
389 throw Failure(
"Vector<Treal, Telement>::"
390 "symv(char const uplo, Treal const, "
391 "const Matrix<Treal, Telement>&, "
392 "const Vector<Treal, Telement>&, "
393 "Treal const, const Vector<Treal, Telement>&):"
394 "Incorrect dimensions for symmetric "
395 "matrix-vector product");
397 throw Failure(
"Vector<class Treal, class Telement>::"
398 "symv only implemented for symmetric matrices in "
399 "upper triangular storage");
400 if ((A.is_zero() || x.
is_zero() || alpha == 0) &&
404 Treal beta_tmp = beta;
409 if (A.is_zero() || x.
is_zero() || alpha == 0)
414 #pragma omp parallel if(A.level() == Params::getMatrixParallelLevel()) num_threads(Params::getNProcs()) default(shared)
419 #pragma omp for schedule(dynamic)
421 for (
int rc = 0; rc < A.
ncols(); rc++) {
428 #pragma omp for schedule(dynamic)
430 for (
int row = 0; row < A.
nrows() - 1; row++) {
432 for (
int col = row + 1; col < A.
ncols(); col++)
438 #pragma omp for schedule(dynamic)
440 for (
int row = 1; row < A.
nrows(); row++) {
442 for (
int col = 0; col < row; col++)
452 template<
class Treal,
class Telement>
453 template<
typename TmatrixElement>
455 trmv(
char const uplo,
const bool tA,
459 throw Failure(
"Vector<Treal, Telement>::"
461 "Incorrect dimensions for triangular "
462 "matrix-vector product");
464 throw Failure(
"Vector<class Treal, class Telement>::"
465 "trmv only implemented for upper triangular matrices");
466 if ( ( A.is_zero() || x.
is_zero() ) ) {
472 for (
int row = 0; row < A.
nrows(); row++) {
474 for (
int col = row + 1; col < A.
ncols(); col++)
480 for (
int col = A.
ncols() - 1; col >= 0; col--) {
482 for (
int row = 0; row < col; row++)
496 template<
class Treal>
507 for (
int ind = 0; ind < this->
n(); ind++)
534 (*this) *= 1 / this->
eucl();
550 static void axpy(Treal
const & alpha,
559 static void gemv(
bool const tA, Treal
const alpha,
568 static void symv(
char const uplo, Treal
const alpha,
578 static void trmv(
char const uplo,
const bool tA,
585 template<
class Treal>
591 template<
class Treal>
600 for (
int row = 0; row < this->
n(); ++row )
604 template<
class Treal>
609 for (
int row = 0; row < this->
nScalars(); ++row )
612 for (
int row = 0; row < this->
n(); ++row )
617 template<
class Treal>
624 template<
class Treal>
630 char * tmp = (
char*)&ZERO;
631 file.write(tmp,
sizeof(
int));
634 char * tmp = (
char*)&ONE;
635 file.write(tmp,
sizeof(
int));
636 char * tmpel = (
char*)this->
elements;
637 file.write(tmpel,
sizeof(Treal) * this->
n());
641 template<
class Treal>
646 char tmp[
sizeof(int)];
647 file.read(tmp, (std::ifstream::pos_type)
sizeof(
int));
655 file.read((
char*)this->
elements,
sizeof(Treal) * this->
n());
658 throw Failure(
"Vector<Treal>::"
659 "readFromFile(std::ifstream & file):"
660 "File corruption, int value not 0 or 1");
664 template<
class Treal>
670 throw Failure(
"Vector::operator=(int k) only implemented for k = 0");
674 template<
class Treal>
678 for (
int ind = 0; ind < this->
n(); ind++)
679 (*
this)(ind) = rand() / (Treal)RAND_MAX;
683 template<
class Treal>
686 if (!this->
is_zero() && alpha != 1) {
693 template<
class Treal>
697 assert(x.
n() == y.
n());
707 template<
class Treal>
712 assert(x.
n() == y.
n());
717 for (
int ind = 0; ind < y.
n(); ind++)
730 template<
class Treal>
732 gemv(
bool const tA, Treal
const alpha,
741 if ((A.is_zero() || x.
is_zero() || alpha == 0) &&
745 Treal beta_tmp = beta;
750 if (A.is_zero() || x.
is_zero() || alpha == 0)
756 throw Failure(
"Vector<Treal, Telement>::"
757 "gemv(bool const, Treal const, "
758 "const Matrix<Treal, Telement>&, "
759 "const Vector<Treal, Telement>&, "
760 "Treal const, const Vector<Treal, "
762 "Incorrect dimensions for matrix-vector product");
771 throw Failure(
"Vector<Treal, Telement>::"
772 "gemv(bool const, Treal const, "
773 "const Matrix<Treal, Telement>&, "
774 "const Vector<Treal, Telement>&, "
775 "Treal const, const Vector<Treal, "
777 "Incorrect dimensions for matrix-vector product");
790 template<
class Treal>
792 symv(
char const uplo, Treal
const alpha,
802 throw Failure(
"Vector<Treal>::"
803 "symv(char const uplo, Treal const, "
804 "const Matrix<Treal>&, "
805 "const Vector<Treal>&, "
806 "Treal const, const Vector<Treal>&):"
807 "Incorrect dimensions for symmetric "
808 "matrix-vector product");
809 if ((A.is_zero() || x.
is_zero() || alpha == 0) &&
813 Treal beta_tmp = beta;
818 if (A.is_zero() || x.
is_zero() || alpha == 0)
828 template<
class Treal>
830 trmv(
char const uplo,
const bool tA,
834 throw Failure(
"Vector<Treal>::"
835 "trmv(...): Incorrect dimensions for triangular "
836 "matrix-vector product");
838 throw Failure(
"Vector<class Treal>::"
839 "trmv only implemented for upper triangular matrices");
840 if ( ( A.is_zero() || x.
is_zero() ) ) {