PetscVector.cxx
1 // Copyright (C) 2011-2012, INRIA
2 // Author(s): Marc Fragu
3 //
4 // This file is part of the linear-algebra library Seldon,
5 // http://seldon.sourceforge.net/.
6 //
7 // Seldon is free software; you can redistribute it and/or modify it under the
8 // terms of the GNU Lesser General Public License as published by the Free
9 // Software Foundation; either version 2.1 of the License, or (at your option)
10 // any later version.
11 //
12 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
13 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15 // more details.
16 //
17 // You should have received a copy of the GNU Lesser General Public License
18 // along with Seldon. If not, see http://www.gnu.org/licenses/.
19 
20 
21 #ifndef SELDON_FILE_VECTOR_PETSCVECTOR_CXX
22 
23 
24 #include "PetscVector.hxx"
25 
26 
27 namespace Seldon
28 {
29 
31 
36  template <class T, class Allocator>
37  void PETScVector<T, Allocator>
38  ::Copy(const Vec& petsc_vector)
39  {
40  Clear();
41 
42  int ierr;
43  ierr = VecGetSize(petsc_vector, &this->m_);
44  CHKERRABORT(mpi_communicator_, ierr);
45  PetscObjectGetComm(reinterpret_cast<PetscObject>(petsc_vector),
46  &mpi_communicator_);
47  CHKERRABORT(mpi_communicator_, ierr);
48  ierr = VecDuplicate(petsc_vector, &petsc_vector_);
49  CHKERRABORT(mpi_communicator_, ierr);
50  ierr = VecCopy(petsc_vector, petsc_vector_);
51  CHKERRABORT(mpi_communicator_, ierr);
52  petsc_vector_deallocated_ = false;
53  }
54 
55 
56  /************************
57  * CONVENIENT FUNCTIONS *
58  ************************/
59 
60 
62 
66  template <class T, class Allocator>
67  void PETScVector<T, Allocator>::Zero()
68  {
69  int ierr;
70  ierr = VecSet(petsc_vector_, 0.);
71  CHKERRABORT(mpi_communicator_, ierr);
72  }
73 
74 
76  template <class T, class Allocator>
77  void PETScVector<T, Allocator>::Fill()
78  {
79  int ierr;
80  Vector<int> index(this->m_);
81  Vector<double> data(this->m_);
82  index.Fill();
83  data.Fill();
84  ierr = VecSetValues(petsc_vector_, this->m_, index.GetData(),
85  data.GetData(), INSERT_VALUES);
86  CHKERRABORT(mpi_communicator_, ierr);
87  Flush();
88  }
89 
90 
92 
95  template <class T, class Allocator>
96  template <class T0>
97  void PETScVector<T, Allocator>::Fill(const T0& x)
98  {
99  int ierr;
100  ierr = VecSet(petsc_vector_, double(x));
101  CHKERRABORT(mpi_communicator_, ierr);
102  Flush();
103  }
104 
105 
107 
110  template <class T, class Allocator>
111  void PETScVector<T, Allocator>::FillRand()
112  {
113  int ierr;
114  Vector<int> index(this->m_);
115  Vector<double> data(this->m_);
116  index.Fill();
117  data.FillRand();
118  ierr = VecSetValues(petsc_vector_, this->m_, index.GetData(),
119  data.GetData(), INSERT_VALUES);
120  CHKERRABORT(mpi_communicator_, ierr);
121  Flush();
122  }
123 
124 
125  /*********
126  * NORMS *
127  *********/
128 
129 
131 
134  template <class T, class Allocator>
135  typename PETScVector<T, Allocator>::value_type
136  PETScVector<T, Allocator>::GetNormInf() const
137  {
138  int ierr;
139  value_type res;
140  ierr = VecNorm(petsc_vector_, NORM_INFINITY, &res);
141  CHKERRABORT(mpi_communicator_, ierr);
142  return res;
143  }
144 
145 
147 
150  template <class T, class Allocator>
151  int PETScVector<T, Allocator>::GetNormInfIndex() const
152  {
153  throw Undefined("PETScVector<T, Allocator>::GetNormInfIndex()");
154  }
155 
156 
158  // SEQUENTIAL PETSC VECTOR //
160 
161 
162  /**************************
163  * OUTPUT/INPUT FUNCTIONS *
164  **************************/
165 
166 
168 
174  template <class T, class Allocator>
175  void Vector<T, PETScSeq, Allocator>
176  ::Write(string FileName, bool with_size) const
177  {
178  throw Undefined("Vector<T, PETScSeq, Allocator>"
179  "::Write(string FileName, bool with_size) const");
180  }
181 
182 
184 
190  template <class T, class Allocator>
191  void Vector<T, PETScSeq, Allocator>
192  ::Write(ostream& FileStream, bool with_size) const
193  {
194  throw Undefined("Vector<T, PETScSeq, Allocator>"
195  "::Write(string FileName, bool with_size) const");
196  }
197 
198 
200 
205  template <class T, class Allocator>
206  void Vector<T, PETScSeq, Allocator>::WriteText(string FileName) const
207  {
208  ofstream FileStream;
209  FileStream.precision(cout.precision());
210  FileStream.flags(cout.flags());
211  FileStream.open(FileName.c_str());
212 #ifdef SELDON_CHECK_IO
213  // Checks if the file was opened.
214  if (!FileStream.is_open())
215  throw IOError("Vector<T, PETScSeq, Allocator>"
216  "::WriteText(string FileName)",
217  string("Unable to open file \"") + FileName + "\".");
218 #endif
219  this->WriteText(FileStream);
220  FileStream.close();
221  }
222 
223 
225 
230  template <class T, class Allocator>
231  void Vector<T, PETScSeq, Allocator>::WriteText(ostream& FileStream) const
232  {
233  throw Undefined("Vector<T, PETScSeq, Allocator>"
234  "::WriteText(ostream& FileStream) const");
235  }
236 
237 
239 
247  template <class T, class Allocator>
248  void Vector<T, PETScSeq, Allocator>
249  ::Read(string FileName, bool with_size)
250  {
251  ifstream FileStream;
252  FileStream.open(FileName.c_str());
253 #ifdef SELDON_CHECK_IO
254  // Checks if the file was opened.
255  if (!FileStream.is_open())
256  throw IOError("Vector<T, PETScSeq, Allocator>::Read(string FileName)",
257  string("Unable to open file \"") + FileName + "\".");
258 #endif
259  this->Read(FileStream, with_size);
260  FileStream.close();
261  }
262 
263 
265 
273  template <class T, class Allocator>
274  void Vector<T, PETScSeq, Allocator>
275  ::Read(istream& FileStream, bool with_size)
276  {
277  throw Undefined("Vector<T, PETScSeq, Allocator>"
278  "::Read(istream& FileStream, bool with_size)");
279  }
280 
281 
283 
288  template <class T, class Allocator>
289  void Vector<T, PETScSeq, Allocator>::ReadText(string FileName)
290  {
291  throw Undefined("Vector<T, PETScSeq, Allocator>"
292  "::ReadText(string FileName)");
293  }
294 
295 
297 
302  template <class T, class Allocator>
303  void Vector<T, PETScSeq, Allocator>::ReadText(istream& FileStream)
304  {
305  throw Undefined("Vector<T, PETScSeq, Allocator>"
306  "::ReadText(istream& FileStream)");
307  }
308 
309 
311 
316  template <class T, class Allocator>
317  ostream& operator << (ostream& out,
318  const Vector<T, PETScSeq, Allocator>& V)
319  {
320  out << '\t';
321  for (int i = 0; i < V.GetLength() - 1; i++)
322  out << V(i) << '\t';
323  if (V.GetLength() != 0)
324  out << V(V.GetLength() - 1);
325  return out;
326  }
327 
328 
330  // PARALLEL PETSC VECTOR //
332 
333 
334  /**************************
335  * OUTPUT/INPUT FUNCTIONS *
336  **************************/
337 
338 
340 
346  template <class T, class Allocator>
347  void Vector<T, PETScPar, Allocator>
348  ::Write(string FileName, bool with_size) const
349  {
350  throw Undefined("Vector<T, PETScPar, Allocator>"
351  "::Write(string FileName, bool with_size) const");
352  }
353 
354 
356 
362  template <class T, class Allocator>
363  void Vector<T, PETScPar, Allocator>
364  ::Write(ostream& FileStream, bool with_size) const
365  {
366  int local_n;
367  VecGetLocalSize(this->petsc_vector_, &local_n);
368 
369  Vector<int> index(local_n);
370  index.Fill();
371  int i_start, i_end;
372  this->GetProcessorRange(i_start, i_end);
373  for (int i = 0; i < local_n; i++)
374  index(i) += i_start;
375  Vector<T> data(local_n);
376  VecGetValues(this->petsc_vector_, local_n, index.GetData(),
377  data.GetData());
378  data.Write(FileStream, with_size);
379  }
380 
381 
383 
388  template <class T, class Allocator>
389  void Vector<T, PETScPar, Allocator>::WriteText(string FileName) const
390  {
391  ofstream FileStream;
392  FileStream.precision(cout.precision());
393  FileStream.flags(cout.flags());
394  FileStream.open(FileName.c_str());
395 #ifdef SELDON_CHECK_IO
396  // Checks if the file was opened.
397  if (!FileStream.is_open())
398  throw IOError("Vector<T, PETScPar, Allocator>"
399  "::WriteText(string FileName)",
400  string("Unable to open file \"") + FileName + "\".");
401 #endif
402  this->WriteText(FileStream);
403  FileStream.close();
404 
405  }
406 
407 
409 
414  template <class T, class Allocator>
415  void Vector<T, PETScPar, Allocator>::WriteText(ostream& FileStream) const
416  {
417  int local_n;
418  VecGetLocalSize(this->petsc_vector_, &local_n);
419 
420  Vector<int> index(local_n);
421  index.Fill();
422  int i_start, i_end;
423  this->GetProcessorRange(i_start, i_end);
424  for (int i = 0; i < local_n; i++)
425  index(i) += i_start;
426  Vector<T> data(local_n);
427  VecGetValues(this->petsc_vector_, local_n, index.GetData(),
428  data.GetData());
429  data.WriteText(FileStream);
430  }
431 
432 
434 
442  template <class T, class Allocator>
443  void Vector<T, PETScPar, Allocator>
444  ::Read(string FileName, bool with_size)
445  {
446  ifstream FileStream;
447  FileStream.open(FileName.c_str());
448 #ifdef SELDON_CHECK_IO
449  // Checks if the file was opened.
450  if (!FileStream.is_open())
451  throw IOError("Vector<T, PETScPar, Allocator>::Read(string FileName)",
452  string("Unable to open file \"") + FileName + "\".");
453 #endif
454  this->Read(FileStream, with_size);
455  FileStream.close();
456  }
457 
458 
460 
468  template <class T, class Allocator>
469  void Vector<T, PETScPar, Allocator>
470  ::Read(istream& FileStream, bool with_size)
471  {
472  throw Undefined("PETScVector<T, PETScPar, Allocator>"
473  "::Read(istream& FileStream, bool with_size)");
474  }
475 
476 
478 
483  template <class T, class Allocator>
484  void Vector<T, PETScPar, Allocator>::ReadText(string FileName)
485  {
486  throw Undefined("PETScVector<T, PETScPar, Allocator>"
487  "::ReadText(string FileName)");
488  }
489 
490 
492 
497  template <class T, class Allocator>
498  void Vector<T, PETScPar, Allocator>::ReadText(istream& FileStream)
499  {
500  throw Undefined("PETScVector<T, PETScPar, Allocator>"
501  "::ReadText(istream& FileStream)");
502  }
503 
504 
506 
511  template <class T, class Allocator>
512  ostream& operator << (ostream& out,
513  const Vector<T, PETScPar, Allocator>& V)
514  {
515  out << '\t';
516  for (int i = 0; i < V.GetLength() - 1; i++)
517  out << V(i) << '\t';
518  if (V.GetLength() != 0)
519  out << V(V.GetLength() - 1);
520  return out;
521  }
522 
523 
524 } // namespace Seldon.
525 
526 
527 #define SELDON_FILE_PETSCVECTOR_CXX
528 #endif
Seldon::Vector< T, PETScPar, Allocator >
Definition: PetscVector.hxx:162
Seldon::Vector< int >
Seldon::PETScVector::Copy
void Copy(const PETScVector< T, Allocator > &X)
Duplicates a vector.
Definition: PetscVectorInline.cxx:296
Seldon::PETScVector::GetNormInf
value_type GetNormInf() const
Returns the infinite norm.
Definition: PetscVector.cxx:136
Seldon::Undefined
Definition: Errors.hxx:62
Seldon::PETScVector::Zero
void Zero()
Sets all elements to zero.
Definition: PetscVector.cxx:67
Seldon::Vector_Base::GetLength
size_t GetLength() const
Returns the number of elements.
Definition: VectorInline.cxx:142
Seldon::PETScVector::FillRand
void FillRand()
Fills the vector randomly.
Definition: PetscVector.cxx:111
Seldon::PETScVector::Fill
void Fill()
Fills the vector with 0, 1, 2, ...
Definition: PetscVector.cxx:77
Seldon::Vector< T, PETScSeq, Allocator >
Definition: PetscVector.hxx:105
Seldon
Seldon namespace.
Definition: Array.cxx:24
Seldon::PETScVector::GetNormInfIndex
int GetNormInfIndex() const
Returns the index of the highest absolute value.
Definition: PetscVector.cxx:151
Seldon::Vector< T, PETScPar, Allocator >::Write
void Write(string FileName, bool with_size=true) const
Writes the vector in a file.
Definition: PetscVector.cxx:348
Seldon::operator<<
ostream & operator<<(ostream &out, const Array< T, N, Allocator > &A)
operator<< overloaded for a 3D array.
Definition: Array.cxx:1617
Seldon::IOError
Definition: Errors.hxx:150