PEXSI
 All Classes Namespaces Files Functions Variables Typedefs Pages
utility.hpp
Go to the documentation of this file.
1 /*
2  Copyright (c) 2012 The Regents of the University of California,
3  through Lawrence Berkeley National Laboratory.
4 
5  Author: Lin Lin and Mathias Jacquelin
6 
7  This file is part of PEXSI. All rights reserved.
8 
9  Redistribution and use in source and binary forms, with or without
10  modification, are permitted provided that the following conditions are met:
11 
12  (1) Redistributions of source code must retain the above copyright notice, this
13  list of conditions and the following disclaimer.
14  (2) Redistributions in binary form must reproduce the above copyright notice,
15  this list of conditions and the following disclaimer in the documentation
16  and/or other materials provided with the distribution.
17  (3) Neither the name of the University of California, Lawrence Berkeley
18  National Laboratory, U.S. Dept. of Energy nor the names of its contributors may
19  be used to endorse or promote products derived from this software without
20  specific prior written permission.
21 
22  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
23  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
24  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
25  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
26  ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
27  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
28  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
29  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 
33  You are under no obligation whatsoever to provide any bug fixes, patches, or
34  upgrades to the features, functionality or performance of the source code
35  ("Enhancements") to anyone; however, if you choose to make your Enhancements
36  available either publicly, or directly to Lawrence Berkeley National
37  Laboratory, without imposing a separate written license agreement for such
38  Enhancements, then you hereby grant the following license: a non-exclusive,
39  royalty-free perpetual license to install, use, modify, prepare derivative
40  works, incorporate into other computer software, distribute, and sublicense
41  such enhancements or derivative works thereof, in binary and source code form.
42 */
46 #ifndef _PEXSI_UTILITY_HPP_
47 #define _PEXSI_UTILITY_HPP_
48 
49 #include <stdlib.h>
50 #include "environment.hpp"
51 #include "tinyvec.hpp"
52 #include "NumVec.hpp"
53 #include "NumMat.hpp"
54 #include "NumTns.hpp"
55 #include "sparse_matrix.hpp"
56 #include "mpi_interf.hpp"
57 
58 namespace PEXSI{
59 
60  // *********************************************************************
61  // Define constants
62  // *********************************************************************
63 
64  // Write format control parameters
65  const int LENGTH_VAR_NAME = 8;
66  const int LENGTH_DBL_DATA = 16;
67  const int LENGTH_INT_DATA = 5;
68  const int LENGTH_VAR_UNIT = 6;
69  const int LENGTH_DBL_PREC = 8;
70  const int LENGTH_VAR_DATA = 16;
71 
72 
73  // standard case for most serialization/deserialization process.
74  const std::vector<Int> NO_MASK(1);
75 
76  // *********************************************************************
77  // Formatted output stream
78  // *********************************************************************
79 
80  // Bool is NOT defined due to ambiguity with Int
81 
82  inline Int PrintBlock(std::ostream &os, const std::string name){
83 
84  os << std::endl<< "*********************************************************************" << std::endl;
85  os << name << std::endl;
86  os << "*********************************************************************" << std::endl << std::endl;
87  return 0;
88  }
89 
90  // String
91  inline Int Print(std::ostream &os, const std::string name) {
92  os << std::setiosflags(std::ios::left) << name << std::endl;
93  return 0;
94  };
95 
96  inline Int Print(std::ostream &os, const char* name) {
97  os << std::setiosflags(std::ios::left) << std::string(name) << std::endl;
98  return 0;
99  };
100 
101  inline Int Print(std::ostream &os, const std::string name, std::string val) {
102  os << std::setiosflags(std::ios::left)
103  << std::setw(LENGTH_VAR_NAME) << name
104  << std::setw(LENGTH_VAR_DATA) << val
105  << std::endl;
106  return 0;
107  };
108 
109  inline Int Print(std::ostream &os, const std::string name, const char* val) {
110  os << std::setiosflags(std::ios::left)
111  << std::setw(LENGTH_VAR_NAME) << name
112  << std::setw(LENGTH_VAR_DATA) << std::string(val)
113  << std::endl;
114  return 0;
115  };
116 
117 
118  // Real
119 
120  // one real number
121 
122  inline Int Print(std::ostream &os, const std::string name, Real val) {
123  os << std::setiosflags(std::ios::left)
124  << std::setw(LENGTH_VAR_NAME) << name
125  << std::setiosflags(std::ios::scientific)
126  << std::setiosflags(std::ios::showpos)
127  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
128  << std::resetiosflags(std::ios::scientific)
129  << std::resetiosflags(std::ios::showpos)
130  << std::endl;
131  return 0;
132  };
133 
134  inline Int Print(std::ostream &os, const char* name, Real val) {
135  os << std::setiosflags(std::ios::left)
136  << std::setw(LENGTH_VAR_NAME) << std::string(name)
137  << std::setiosflags(std::ios::scientific)
138  << std::setiosflags(std::ios::showpos)
139  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
140  << std::resetiosflags(std::ios::scientific)
141  << std::resetiosflags(std::ios::showpos)
142  << std::endl;
143  return 0;
144  };
145 
146 
147  inline Int Print(std::ostream &os, const std::string name, Real val, const std::string unit) {
148  os << std::setiosflags(std::ios::left)
149  << std::setw(LENGTH_VAR_NAME) << name
150  << std::setiosflags(std::ios::scientific)
151  << std::setiosflags(std::ios::showpos)
152  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
153  << std::resetiosflags(std::ios::scientific)
154  << std::resetiosflags(std::ios::showpos)
155  << std::setw(LENGTH_VAR_UNIT) << unit
156  << std::endl;
157  return 0;
158  };
159 
160  inline Int Print(std::ostream &os, const char *name, Real val, const char *unit) {
161  os << std::setiosflags(std::ios::left)
162  << std::setw(LENGTH_VAR_NAME) << std::string(name)
163  << std::setiosflags(std::ios::scientific)
164  << std::setiosflags(std::ios::showpos)
165  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
166  << std::resetiosflags(std::ios::scientific)
167  << std::resetiosflags(std::ios::showpos)
168  << std::setw(LENGTH_VAR_UNIT) << std::string(unit)
169  << std::endl;
170  return 0;
171  };
172 
173  // Two real numbers
174  inline Int Print(std::ostream &os, const std::string name1, Real val1, const std::string unit1,
175  const std::string name2, Real val2, const std::string unit2) {
176  os << std::setiosflags(std::ios::left)
177  << std::setw(LENGTH_VAR_NAME) << name1
178  << std::setiosflags(std::ios::scientific)
179  << std::setiosflags(std::ios::showpos)
180  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val1
181  << std::setw(LENGTH_VAR_UNIT) << unit1
182  << std::setw(LENGTH_VAR_NAME) << name2
183  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
184  << std::resetiosflags(std::ios::scientific)
185  << std::resetiosflags(std::ios::showpos)
186  << std::setw(LENGTH_VAR_UNIT) << unit2
187  << std::endl;
188  return 0;
189  };
190 
191  inline Int Print(std::ostream &os, const char *name1, Real val1, const char *unit1,
192  char *name2, Real val2, char *unit2) {
193  os << std::setiosflags(std::ios::left)
194  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
195  << std::setiosflags(std::ios::scientific)
196  << std::setiosflags(std::ios::showpos)
197  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val1
198  << std::setw(LENGTH_VAR_UNIT) << std::string(unit1)
199  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
200  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
201  << std::resetiosflags(std::ios::scientific)
202  << std::resetiosflags(std::ios::showpos)
203  << std::setw(LENGTH_VAR_UNIT) << std::string(unit2)
204  << std::endl;
205  return 0;
206  };
207 
208  // Int and Real
209  inline Int Print(std::ostream &os, const std::string name1, Int val1, const std::string unit1,
210  const std::string name2, Real val2, const std::string unit2) {
211  os << std::setiosflags(std::ios::left)
212  << std::setw(LENGTH_VAR_NAME) << name1
213  << std::setw(LENGTH_INT_DATA) << val1
214  << std::setw(LENGTH_VAR_UNIT) << unit1
215  << std::setw(LENGTH_VAR_NAME) << name2
216  << std::setiosflags(std::ios::scientific)
217  << std::setiosflags(std::ios::showpos)
218  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
219  << std::resetiosflags(std::ios::scientific)
220  << std::resetiosflags(std::ios::showpos)
221  << std::setw(LENGTH_VAR_UNIT) << unit2
222  << std::endl;
223  return 0;
224  };
225 
226  inline Int Print(std::ostream &os, const char *name1, Int val1, const char *unit1,
227  char *name2, Real val2, char *unit2) {
228  os << std::setiosflags(std::ios::left)
229  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
230  << std::setw(LENGTH_INT_DATA) << val1
231  << std::setw(LENGTH_VAR_UNIT) << std::string(unit1)
232  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
233  << std::setiosflags(std::ios::scientific)
234  << std::setiosflags(std::ios::showpos)
235  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
236  << std::resetiosflags(std::ios::scientific)
237  << std::resetiosflags(std::ios::showpos)
238  << std::setw(LENGTH_VAR_UNIT) << std::string(unit2)
239  << std::endl;
240  return 0;
241  };
242 
243  inline Int Print(std::ostream &os,
244  const char *name1, Int val1,
245  const char *name2, Real val2,
246  char *name3, Real val3) {
247  os << std::setiosflags(std::ios::left)
248  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
249  << std::setw(LENGTH_INT_DATA) << val1
250  << std::setiosflags(std::ios::scientific)
251  << std::setiosflags(std::ios::showpos)
252  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
253  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
254  << std::setw(LENGTH_VAR_NAME) << std::string(name3)
255  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val3
256  << std::resetiosflags(std::ios::scientific)
257  << std::resetiosflags(std::ios::showpos)
258  << std::endl;
259  return 0;
260  };
261 
262 
263  inline Int Print(std::ostream &os,
264  const char *name1, Int val1,
265  const char *name2, Real val2,
266  const char *name3, Real val3,
267  const char *name4, Real val4 ) {
268  os << std::setiosflags(std::ios::left)
269  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
270  << std::setw(LENGTH_INT_DATA) << val1
271  << std::setiosflags(std::ios::scientific)
272  << std::setiosflags(std::ios::showpos)
273  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
274  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
275  << std::setw(LENGTH_VAR_NAME) << std::string(name3)
276  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val3
277  << std::setw(LENGTH_VAR_NAME) << std::string(name4)
278  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val4
279  << std::resetiosflags(std::ios::scientific)
280  << std::resetiosflags(std::ios::showpos)
281  << std::endl;
282  return 0;
283  };
284 
285  // Int
286 
287  // one integer number
288  inline Int Print(std::ostream &os, std::string name, Int val) {
289  os << std::setiosflags(std::ios::left)
290  << std::setw(LENGTH_VAR_NAME) << name
291  << std::setw(LENGTH_VAR_DATA) << val
292  << std::endl;
293  return 0;
294  };
295 
296 
297  inline Int Print(std::ostream &os, const char *name, Int val) {
298  os << std::setiosflags(std::ios::left)
299  << std::setw(LENGTH_VAR_NAME) << std::string(name)
300  << std::setw(LENGTH_VAR_DATA) << val
301  << std::endl;
302  return 0;
303  };
304 
305 
306  inline Int Print(std::ostream &os, const std::string name, Int val, const std::string unit) {
307  os << std::setiosflags(std::ios::left)
308  << std::setw(LENGTH_VAR_NAME) << name
309  << std::setw(LENGTH_VAR_DATA) << val
310  << std::setw(LENGTH_VAR_UNIT) << unit
311  << std::endl;
312  return 0;
313  };
314 
315 
316  inline Int Print(std::ostream &os, const char* name, Int val, const std::string unit) {
317  os << std::setiosflags(std::ios::left)
318  << std::setw(LENGTH_VAR_NAME) << std::string(name)
319  << std::setw(LENGTH_VAR_DATA) << val
320  << std::setw(LENGTH_VAR_UNIT) << unit
321  << std::endl;
322  return 0;
323  };
324 
325 
326 
327  // two integer numbers
328  inline Int Print(std::ostream &os, const std::string name1, Int val1, const std::string unit1,
329  const std::string name2, Int val2, const std::string unit2) {
330  os << std::setiosflags(std::ios::left)
331  << std::setw(LENGTH_VAR_NAME) << name1
332  << std::setw(LENGTH_VAR_DATA) << val1
333  << std::setw(LENGTH_VAR_UNIT) << unit1
334  << std::setw(LENGTH_VAR_NAME) << name2
335  << std::setw(LENGTH_VAR_DATA) << val2
336  << std::setw(LENGTH_VAR_UNIT) << unit2
337  << std::endl;
338  return 0;
339  };
340 
341  inline Int Print(std::ostream &os, const char *name1, Int val1, const char *unit1,
342  char *name2, Int val2, char *unit2) {
343  os << std::setiosflags(std::ios::left)
344  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
345  << std::setw(LENGTH_VAR_DATA) << val1
346  << std::setw(LENGTH_VAR_UNIT) << std::string(unit1)
347  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
348  << std::setw(LENGTH_VAR_DATA) << val2
349  << std::setw(LENGTH_VAR_UNIT) << std::string(unit2)
350  << std::endl;
351  return 0;
352  };
353 
354  // Bool
355 
356  // one boolean number
357  inline Int Print(std::ostream &os, const std::string name, bool val) {
358  os << std::setiosflags(std::ios::left)
359  << std::setw(LENGTH_VAR_NAME) << name;
360  if( val == true )
361  os << std::setw(LENGTH_VAR_NAME) << "true" << std::endl;
362  else
363  os << std::setw(LENGTH_VAR_NAME) << "false" << std::endl;
364  return 0;
365  };
366 
367 
368  inline Int Print(std::ostream &os, const char* name, bool val) {
369  os << std::setiosflags(std::ios::left)
370  << std::setw(LENGTH_VAR_NAME) << std::string(name);
371  if( val == true )
372  os << std::setw(LENGTH_VAR_NAME) << "true" << std::endl;
373  else
374  os << std::setw(LENGTH_VAR_NAME) << "false" << std::endl;
375  return 0;
376  };
377 
378 
379  // Index 3 and Point 3
380  inline Int Print(std::ostream &os,
381  const char *name1, Index3 val ) {
382  os << std::setiosflags(std::ios::left)
383  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
384  << std::setw(LENGTH_VAR_DATA) << val[0]
385  << std::setw(LENGTH_VAR_DATA) << val[1]
386  << std::setw(LENGTH_VAR_DATA) << val[2]
387  << std::endl;
388  return 0;
389  };
390 
391  inline Int Print(std::ostream &os,
392  const char *name1, Point3 val ) {
393  os << std::setiosflags(std::ios::left)
394  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
395  << std::setiosflags(std::ios::scientific)
396  << std::setiosflags(std::ios::showpos)
397  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) << val[0]
398  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) << val[1]
399  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) << val[2]
400  << std::resetiosflags(std::ios::scientific)
401  << std::resetiosflags(std::ios::showpos)
402  << std::endl;
403  return 0;
404  };
405 
406  inline Int Print(std::ostream &os,
407  const char *name1, Int val1,
408  const char *name2, Point3 val ) {
409  os << std::setiosflags(std::ios::left)
410  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
411  << std::setw(LENGTH_INT_DATA) << val1
412  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
413  << std::setiosflags(std::ios::scientific)
414  << std::setiosflags(std::ios::showpos)
415  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) <<val[0]
416  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) <<val[1]
417  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) <<val[2]
418  << std::resetiosflags(std::ios::scientific)
419  << std::resetiosflags(std::ios::showpos)
420  << std::endl;
421  return 0;
422  };
423 
424  // *********************************************************************
425  // Overload << and >> operators for basic data types
426  // *********************************************************************
427 
428  // std::vector
429  template <class F> inline std::ostream& operator<<( std::ostream& os, const std::vector<F>& vec)
430  {
431  os<<vec.size()<<std::endl;
432  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
433  for(Int i=0; i<vec.size(); i++)
434  os<<" "<<vec[i];
435  os<<std::endl;
436  return os;
437  }
438 
439  // NumVec
440  template <class F> inline std::ostream& operator<<( std::ostream& os, const NumVec<F>& vec)
441  {
442  os<<vec.m()<<std::endl;
443  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
444  for(Int i=0; i<vec.m(); i++)
445  os<<" "<<vec(i);
446  os<<std::endl;
447  return os;
448  }
449 
450  //template <class F> inline std::istream& operator>>( std::istream& is, NumVec<F>& vec)
451  //{
452  // Int m; is>>m; vec.resize(m);
453  // for(Int i=0; i<vec.m(); i++)
454  // is >> vec(i);
455  // return is;
456  //}
457 
458  // NumMat
459  template <class F> inline std::ostream& operator<<( std::ostream& os, const NumMat<F>& mat)
460  {
461  os<<mat.m()<<" "<<mat.n()<<std::endl;
462  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
463  os << std::setprecision(16);
464  for(Int i=0; i<mat.m(); i++) {
465  for(Int j=0; j<mat.n(); j++)
466  os<<" "<<mat(i,j);
467  os<<std::endl;
468  }
469  return os;
470  }
471 
472  // NumTns
473  template <class F> inline std::ostream& operator<<( std::ostream& os, const NumTns<F>& tns)
474  {
475  os<<tns.m()<<" "<<tns.n()<<" "<<tns.p()<<std::endl;
476  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
477  for(Int i=0; i<tns.m(); i++) {
478  for(Int j=0; j<tns.n(); j++) {
479  for(Int k=0; k<tns.p(); k++) {
480  os<<" "<<tns(i,j,k);
481  }
482  os<<std::endl;
483  }
484  os<<std::endl;
485  }
486  return os;
487  }
488 
489  // *********************************************************************
490  // serialize/deserialize for basic types
491  // More specific serialize/deserialize will be defined in individual
492  // class files
493  // *********************************************************************
494 
495 
496 
497 
498  //template
499  template <typename T>
500  inline Int serialize(const T& val, std::ostream& os, const std::vector<Int>& mask)
501  {
502  os.write((char*)&val, sizeof(T));
503  return 0;
504  }
505 
506  template <typename T>
507  inline Int deserialize(T& val, std::istream& is, const std::vector<Int>& mask)
508  {
509  is.read((char*)&val, sizeof(T));
510  return 0;
511  }
512 
513 
514 
515  //bool
516  inline Int serialize(const bool& val, std::ostream& os, const std::vector<Int>& mask)
517  {
518  os.write((char*)&val, sizeof(bool));
519  return 0;
520  }
521 
522  inline Int deserialize(bool& val, std::istream& is, const std::vector<Int>& mask)
523  {
524  is.read((char*)&val, sizeof(bool));
525  return 0;
526  }
527 
528  //char
529  inline Int serialize(const char& val, std::ostream& os, const std::vector<Int>& mask)
530  {
531  os.write((char*)&val, sizeof(char));
532  return 0;
533  }
534 
535  inline Int deserialize(char& val, std::istream& is, const std::vector<Int>& mask)
536  {
537  is.read((char*)&val, sizeof(char));
538  return 0;
539  }
540 
541  inline Int combine(char& val, char& ext)
542  {
543  throw std::logic_error( "Combine operation not implemented." );
544  }
545 
546  //-------------------
547  //Int
548  inline Int serialize(const Int& val, std::ostream& os, const std::vector<Int>& mask)
549  {
550  os.write((char*)&val, sizeof(Int));
551  return 0;
552  }
553 
554  inline Int deserialize(Int& val, std::istream& is, const std::vector<Int>& mask)
555  {
556  is.read((char*)&val, sizeof(Int));
557  return 0;
558  }
559 
560  inline Int combine(Int& val, Int& ext)
561  {
562  val += ext;
563  return 0;
564  }
565 
566  //-------------------
567  //LongInt
568  inline Int serialize(const LongInt& val, std::ostream& os, const std::vector<Int>& mask)
569  {
570  os.write((char*)&val, sizeof(LongInt));
571  return 0;
572  }
573 
574  inline Int deserialize(LongInt& val, std::istream& is, const std::vector<Int>& mask)
575  {
576  is.read((char*)&val, sizeof(LongInt));
577  return 0;
578  }
579 
580  inline Int combine(LongInt& val, LongInt& ext)
581  {
582  val += ext;
583  return 0;
584  }
585 
586 
587  //-------------------
588  //Real
589  inline Int serialize(const Real& val, std::ostream& os, const std::vector<Int>& mask)
590  {
591  os.write((char*)&val, sizeof(Real));
592  return 0;
593  }
594 
595  inline Int deserialize(Real& val, std::istream& is, const std::vector<Int>& mask)
596  {
597  is.read((char*)&val, sizeof(Real));
598  return 0;
599  }
600 
601  inline Int combine(Real& val, Real& ext)
602  {
603  val += ext;
604  return 0;
605  }
606 
607  //-------------------
608  //Complex
609  inline Int serialize(const Complex& val, std::ostream& os, const std::vector<Int>& mask)
610  {
611  os.write((char*)&val, sizeof(Complex));
612  return 0;
613  }
614 
615  inline Int deserialize(Complex& val, std::istream& is, const std::vector<Int>& mask)
616  {
617  is.read((char*)&val, sizeof(Complex));
618  return 0;
619  }
620 
621  inline Int combine(Complex& val, Complex& ext)
622  {
623  val += ext;
624  return 0;
625  }
626 
627  //-------------------
628  //Index2 TODO
629  //inline Int serialize(const Index2& val, std::ostream& os, const std::vector<Int>& mask)
630  //{
631  // os.write((char*)&(val[0]), 2*sizeof(Int));
632  // return 0;
633  //}
634  //
635  //inline Int deserialize(Index2& val, std::istream& is, const std::vector<Int>& mask)
636  //{
637  // is.read((char*)&(val[0]), 2*sizeof(Int));
638  // return 0;
639  //}
640  //
641  //inline Int combine(Index2& val, Index2& ext)
642  //{
643  // throw std::logic_error( "Combine operation not implemented." );
644  // return 0;
645  //}
646 
647  //-------------------
648  //Point2 TODO
649  //inline Int serialize(const Point2& val, std::ostream& os, const std::vector<Int>& mask)
650  //{
651  // os.write((char*)&(val[0]), 2*sizeof(Real));
652  // return 0;
653  //}
654  //
655  //inline Int deserialize(Point2& val, std::istream& is, const std::vector<Int>& mask)
656  //{
657  // is.read((char*)&(val[0]), 2*sizeof(Real));
658  // return 0;
659  //}
660  //
661  //inline Int combine(Point2& val, Point2& ext)
662  //{
663  // throw std::logic_error( "Combine operation not implemented." );
664  // return 0;
665  //}
666 
667  //-------------------
668  //Index3
669  inline Int serialize(const Index3& val, std::ostream& os, const std::vector<Int>& mask)
670  {
671  os.write((char*)&(val[0]), 3*sizeof(Int));
672  return 0;
673  }
674 
675  inline Int deserialize(Index3& val, std::istream& is, const std::vector<Int>& mask)
676  {
677  is.read((char*)&(val[0]), 3*sizeof(Int));
678  return 0;
679  }
680 
681  inline Int combine(Index3& val, Index3& ext)
682  {
683  throw std::logic_error( "Combine operation not implemented." );
684  }
685 
686  //-------------------
687  //Point3
688  inline Int serialize(const Point3& val, std::ostream& os, const std::vector<Int>& mask)
689  {
690  os.write((char*)&(val[0]), 3*sizeof(Real));
691  return 0;
692  }
693 
694  inline Int deserialize(Point3& val, std::istream& is, const std::vector<Int>& mask)
695  {
696  is.read((char*)&(val[0]), 3*sizeof(Real));
697  return 0;
698  }
699 
700  inline Int combine(Point3& val, Point3& ext)
701  {
702  throw std::logic_error( "Combine operation not implemented." );
703  }
704 
705  //-------------------
706  //std::vector
707  template<class T>
708  Int serialize(const std::vector<T>& val, std::ostream& os, const std::vector<Int>& mask)
709  {
710  Int sz = val.size();
711  os.write((char*)&sz, sizeof(Int));
712  for(Int k=0; k<sz; k++)
713  serialize(val[k], os, mask);
714  return 0;
715  }
716 
717  template<class T>
718  Int deserialize(std::vector<T>& val, std::istream& is, const std::vector<Int>& mask)
719  {
720  Int sz;
721  is.read((char*)&sz, sizeof(Int));
722  val.resize(sz);
723  for(Int k=0; k<sz; k++)
724  deserialize(val[k], is, mask);
725  return 0;
726  }
727 
728  template<class T>
729  Int combine(std::vector<T>& val, std::vector<T>& ext)
730  {
731  throw std::logic_error( "Combine operation not implemented." );
732  return 0;
733  }
734 
735  //-------------------
736  //std::set
737  template<class T>
738  Int serialize(const std::set<T>& val, std::ostream& os, const std::vector<Int>& mask)
739  {
740  Int sz = val.size();
741  os.write((char*)&sz, sizeof(Int));
742  for(typename std::set<T>::const_iterator mi=val.begin(); mi!=val.end(); mi++)
743  serialize((*mi), os, mask);
744  return 0;
745  }
746 
747  template<class T>
748  Int deserialize(std::set<T>& val, std::istream& is, const std::vector<Int>& mask)
749  {
750  val.clear();
751  Int sz;
752  is.read((char*)&sz, sizeof(Int));
753  for(Int k=0; k<sz; k++) {
754  T t; deserialize(t, is, mask);
755  val.insert(t);
756  }
757  return 0;
758  }
759 
760  template<class T>
761  Int combine(std::set<T>& val, std::set<T>& ext)
762  {
763  throw std::logic_error( "Combine operation not implemented." );
764  return 0;
765  }
766 
767  //-------------------
768  //std::map
769  template<class T, class S>
770  Int serialize(const std::map<T,S>& val, std::ostream& os, const std::vector<Int>& mask)
771  {
772  Int sz = val.size();
773  os.write((char*)&sz, sizeof(Int));
774  for(typename std::map<T,S>::const_iterator mi=val.begin(); mi!=val.end(); mi++) {
775  serialize((*mi).first, os, mask);
776  serialize((*mi).second, os, mask);
777  }
778  return 0;
779  }
780 
781  template<class T, class S>
782  Int deserialize(std::map<T,S>& val, std::istream& is, const std::vector<Int>& mask)
783  {
784  val.clear();
785  Int sz;
786  is.read((char*)&sz, sizeof(Int));
787  for(Int k=0; k<sz; k++) {
788  T t; deserialize(t, is, mask);
789  S s; deserialize(s, is, mask);
790  val[t] = s;
791  }
792  return 0;
793  }
794 
795  template<class T, class S>
796  Int combine(std::map<T,S>& val, std::map<T,S>& ext)
797  {
798  throw std::logic_error( "Combine operation not implemented." );
799  return 0;
800  }
801 
802  //-------------------
803  //std::pair
804  template<class T, class S>
805  Int serialize(const std::pair<T,S>& val, std::ostream& os, const std::vector<Int>& mask)
806  {
807  serialize(val.first, os, mask);
808  serialize(val.second, os, mask);
809  return 0;
810  }
811 
812  template<class T, class S>
813  Int deserialize(std::pair<T,S>& val, std::istream& is, const std::vector<Int>& mask)
814  {
815  deserialize(val.first, is, mask);
816  deserialize(val.second, is, mask);
817  return 0;
818  }
819 
820  template<class T, class S>
821  Int combine(std::pair<T,S>& val, std::pair<T,S>& ext)
822  {
823  throw std::logic_error( "Combine operation not implemented." );
824  return 0;
825  }
826 
827  /*
828  //-------------------
829  //BolNumVec
830  inline Int serialize(const BolNumVec& val, std::ostream& os, const std::vector<Int>& mask)
831  {
832  Int m = val.m();
833  os.write((char*)&m, sizeof(Int));
834  os.write((char*)(val.Data()), m*sizeof(bool));
835  return 0;
836  }
837 
838  inline Int deserialize(BolNumVec& val, std::istream& is, const std::vector<Int>& mask)
839  {
840  Int m;
841  is.read((char*)&m, sizeof(Int));
842  val.Resize(m);
843  is.read((char*)(val.Data()), m*sizeof(bool));
844  return 0;
845  }
846 
847  //-------------------
848  //BolNumMat
849  inline Int serialize(const BolNumMat& val, std::ostream& os, const std::vector<Int>& mask)
850  {
851  Int m = val.m();
852  Int n = val.n();
853  os.write((char*)&m, sizeof(Int));
854  os.write((char*)&n, sizeof(Int));
855  os.write((char*)(val.Data()), m*n*sizeof(bool));
856  return 0;
857  }
858 
859  inline Int deserialize(BolNumMat& val, std::istream& is, const std::vector<Int>& mask)
860  {
861  Int m;
862  Int n;
863  is.read((char*)&m, sizeof(Int));
864  is.read((char*)&n, sizeof(Int));
865  val.Resize(m,n);
866  is.read((char*)(val.Data()), m*n*sizeof(bool));
867  return 0;
868  }
869 
870  //-------------------
871  //BolNumTns
872  inline Int serialize(const BolNumTns& val, std::ostream& os, const std::vector<Int>& mask)
873  {
874  Int m = val.m(); Int n = val.n(); Int p = val.p();
875  os.write((char*)&m, sizeof(Int));
876  os.write((char*)&n, sizeof(Int));
877  os.write((char*)&p, sizeof(Int));
878  os.write((char*)(val.Data()), m*n*p*sizeof(bool));
879  return 0;
880  }
881 
882  inline Int deserialize(BolNumTns& val, std::istream& is, const std::vector<Int>& mask)
883  {
884  Int m,n,p;
885  is.read((char*)&m, sizeof(Int));
886  is.read((char*)&n, sizeof(Int));
887  is.read((char*)&p, sizeof(Int));
888  val.Resize(m,n,p);
889  is.read((char*)(val.Data()), m*n*p*sizeof(bool));
890  return 0;
891  }
892  */
893 
894  //-------------------
895  //IntNumVec
896  inline Int serialize(const IntNumVec& val, std::ostream& os, const std::vector<Int>& mask)
897  {
898  Int m = val.m();
899  os.write((char*)&m, sizeof(Int));
900  os.write((char*)(val.Data()), m*sizeof(Int));
901  return 0;
902  }
903 
904  inline Int deserialize(IntNumVec& val, std::istream& is, const std::vector<Int>& mask)
905  {
906  Int m;
907  is.read((char*)&m, sizeof(Int));
908  val.Resize(m);
909  is.read((char*)(val.Data()), m*sizeof(Int));
910  return 0;
911  }
912 
913  inline Int combine(IntNumVec& val, IntNumVec& ext)
914  {
915  //val.resize(ext.m());
916  assert(val.m()==ext.m());
917  for(Int i=0; i<val.m(); i++) val(i) += ext(i);
918  return 0;
919  }
920 
921 
922  //-------------------
923  //IntNumMat
924  inline Int serialize(const IntNumMat& val, std::ostream& os, const std::vector<Int>& mask)
925  {
926  Int m = val.m();
927  Int n = val.n();
928  os.write((char*)&m, sizeof(Int));
929  os.write((char*)&n, sizeof(Int));
930  os.write((char*)(val.Data()), m*n*sizeof(Int));
931  return 0;
932  }
933 
934  inline Int deserialize(IntNumMat& val, std::istream& is, const std::vector<Int>& mask)
935  {
936  Int m;
937  Int n;
938  is.read((char*)&m, sizeof(Int));
939  is.read((char*)&n, sizeof(Int));
940  val.Resize(m,n);
941  is.read((char*)(val.Data()), m*n*sizeof(Int));
942  return 0;
943  }
944 
945  inline Int combine(IntNumMat& val, IntNumMat& ext)
946  {
947  //val.resize(ext.m(),ext.n());
948  assert(val.m()==ext.m() && val.n()==ext.n());
949  for(Int i=0; i<val.m(); i++)
950  for(Int j=0; j<val.n(); j++)
951  val(i,j) += ext(i,j);
952  return 0;
953  }
954 
955  //-------------------
956  //IntNumTns
957  inline Int serialize(const IntNumTns& val, std::ostream& os, const std::vector<Int>& mask)
958  {
959  Int m = val.m(); Int n = val.n(); Int p = val.p();
960  os.write((char*)&m, sizeof(Int));
961  os.write((char*)&n, sizeof(Int));
962  os.write((char*)&p, sizeof(Int));
963  os.write((char*)(val.Data()), m*n*p*sizeof(Int));
964  return 0;
965  }
966 
967  inline Int deserialize(IntNumTns& val, std::istream& is, const std::vector<Int>& mask)
968  {
969  Int m,n,p;
970  is.read((char*)&m, sizeof(Int));
971  is.read((char*)&n, sizeof(Int));
972  is.read((char*)&p, sizeof(Int));
973  val.Resize(m,n,p);
974  is.read((char*)(val.Data()), m*n*p*sizeof(Int));
975  return 0;
976  }
977 
978  inline Int combine(IntNumTns& val, IntNumTns& ext)
979  {
980  //val.resize(ext.m(),ext.n(),ext.p());
981  assert(val.m()==ext.m() && val.n()==ext.n() && val.p()==ext.p());
982  for(Int i=0; i<val.m(); i++)
983  for(Int j=0; j<val.n(); j++)
984  for(Int k=0; k<val.p(); k++)
985  val(i,j,k) += ext(i,j,k);
986  return 0;
987  }
988 
989  //-------------------
990  //DblNumVec
991  inline Int serialize(const DblNumVec& val, std::ostream& os, const std::vector<Int>& mask)
992  {
993  Int m = val.m();
994  os.write((char*)&m, sizeof(Int));
995  os.write((char*)(val.Data()), m*sizeof(Real));
996  return 0;
997  }
998 
999  inline Int deserialize(DblNumVec& val, std::istream& is, const std::vector<Int>& mask)
1000  {
1001  Int m;
1002  is.read((char*)&m, sizeof(Int));
1003  val.Resize(m);
1004  is.read((char*)(val.Data()), m*sizeof(Real));
1005  return 0;
1006  }
1007 
1008  inline Int combine(DblNumVec& val, DblNumVec& ext)
1009  {
1010  //val.resize(ext.m());
1011  assert(val.m()==ext.m());
1012  for(Int i=0; i<val.m(); i++) val(i) += ext(i);
1013  return 0;
1014  }
1015 
1016  //-------------------
1017  //DblNumMat
1018  inline Int serialize(const DblNumMat& val, std::ostream& os, const std::vector<Int>& mask)
1019  {
1020  Int m = val.m();
1021  Int n = val.n();
1022  os.write((char*)&m, sizeof(Int));
1023  os.write((char*)&n, sizeof(Int));
1024  os.write((char*)(val.Data()), m*n*sizeof(Real));
1025  return 0;
1026  }
1027 
1028  inline Int deserialize(DblNumMat& val, std::istream& is, const std::vector<Int>& mask)
1029  {
1030  Int m;
1031  Int n;
1032  is.read((char*)&m, sizeof(Int));
1033  is.read((char*)&n, sizeof(Int));
1034  val.Resize(m,n);
1035  is.read((char*)(val.Data()), m*n*sizeof(Real));
1036  return 0;
1037  }
1038 
1039  inline Int combine(DblNumMat& val, DblNumMat& ext)
1040  {
1041  //val.resize(ext.m(),ext.n());
1042  assert(val.m()==ext.m() && val.n()==ext.n());
1043  for(Int i=0; i<val.m(); i++)
1044  for(Int j=0; j<val.n(); j++)
1045  val(i,j) += ext(i,j);
1046  return 0;
1047  }
1048 
1049 
1050  //-------------------
1051  //DblNumTns
1052  inline Int serialize(const DblNumTns& val, std::ostream& os, const std::vector<Int>& mask)
1053  {
1054  Int m = val.m(); Int n = val.n(); Int p = val.p();
1055  os.write((char*)&m, sizeof(Int));
1056  os.write((char*)&n, sizeof(Int));
1057  os.write((char*)&p, sizeof(Int));
1058  os.write((char*)(val.Data()), m*n*p*sizeof(Real));
1059  return 0;
1060  }
1061 
1062  inline Int deserialize(DblNumTns& val, std::istream& is, const std::vector<Int>& mask)
1063  {
1064  Int m,n,p;
1065  is.read((char*)&m, sizeof(Int));
1066  is.read((char*)&n, sizeof(Int));
1067  is.read((char*)&p, sizeof(Int));
1068  val.Resize(m,n,p);
1069  is.read((char*)(val.Data()), m*n*p*sizeof(Real));
1070  return 0;
1071  }
1072 
1073  inline Int combine(DblNumTns& val, DblNumTns& ext)
1074  {
1075  //val.resize(ext.m(),ext.n(),ext.p());
1076  assert(val.m()==ext.m() && val.n()==ext.n() && val.p()==ext.p());
1077  for(Int i=0; i<val.m(); i++)
1078  for(Int j=0; j<val.n(); j++)
1079  for(Int k=0; k<val.p(); k++)
1080  val(i,j,k) += ext(i,j,k);
1081  return 0;
1082  }
1083 
1084  //-------------------
1085  //CpxNumVec
1086  inline Int serialize(const CpxNumVec& val, std::ostream& os, const std::vector<Int>& mask)
1087  {
1088  Int m = val.m();
1089  os.write((char*)&m, sizeof(Int));
1090  os.write((char*)(val.Data()), m*sizeof(Complex));
1091  return 0;
1092  }
1093 
1094  inline Int deserialize(CpxNumVec& val, std::istream& is, const std::vector<Int>& mask)
1095  {
1096  Int m;
1097  is.read((char*)&m, sizeof(Int));
1098  val.Resize(m);
1099  is.read((char*)(val.Data()), m*sizeof(Complex));
1100  return 0;
1101  }
1102 
1103  inline Int combine(CpxNumVec& val, CpxNumVec& ext)
1104  {
1105  //val.resize(ext.m());
1106  assert(val.m()==ext.m());
1107  for(Int i=0; i<val.m(); i++) val(i) += ext(i);
1108  return 0;
1109  }
1110 
1111  //-------------------
1112  //CpxNumMat
1113  inline Int serialize(const CpxNumMat& val, std::ostream& os, const std::vector<Int>& mask)
1114  {
1115  Int m = val.m();
1116  Int n = val.n();
1117  os.write((char*)&m, sizeof(Int));
1118  os.write((char*)&n, sizeof(Int));
1119  os.write((char*)(val.Data()), m*n*sizeof(Complex));
1120  return 0;
1121  }
1122 
1123  inline Int deserialize(CpxNumMat& val, std::istream& is, const std::vector<Int>& mask)
1124  {
1125  Int m;
1126  Int n;
1127  is.read((char*)&m, sizeof(Int));
1128  is.read((char*)&n, sizeof(Int));
1129  val.Resize(m,n);
1130  is.read((char*)(val.Data()), m*n*sizeof(Complex));
1131  return 0;
1132  }
1133 
1134  inline Int combine(CpxNumMat& val, CpxNumMat& ext)
1135  {
1136  //val.resize(ext.m(),ext.n());
1137  assert(val.m()==ext.m() && val.n()==ext.n());
1138  for(Int i=0; i<val.m(); i++)
1139  for(Int j=0; j<val.n(); j++)
1140  val(i,j) += ext(i,j);
1141  return 0;
1142  }
1143 
1144  //-------------------
1145  //CpxNumTns
1146  inline Int serialize(const CpxNumTns& val, std::ostream& os, const std::vector<Int>& mask)
1147  {
1148  Int m = val.m(); Int n = val.n(); Int p = val.p();
1149  os.write((char*)&m, sizeof(Int));
1150  os.write((char*)&n, sizeof(Int));
1151  os.write((char*)&p, sizeof(Int));
1152  os.write((char*)(val.Data()), m*n*p*sizeof(Complex));
1153  return 0;
1154  }
1155 
1156  inline Int deserialize(CpxNumTns& val, std::istream& is, const std::vector<Int>& mask)
1157  {
1158  Int m,n,p;
1159  is.read((char*)&m, sizeof(Int));
1160  is.read((char*)&n, sizeof(Int));
1161  is.read((char*)&p, sizeof(Int));
1162  val.Resize(m,n,p);
1163  is.read((char*)(val.Data()), m*n*p*sizeof(Complex));
1164  return 0;
1165  }
1166 
1167  inline Int combine(CpxNumTns& val, CpxNumTns& ext)
1168  {
1169  //val.resize(ext.m(),ext.n(),ext.p());
1170  assert(val.m()==ext.m() && val.n()==ext.n() && val.p()==ext.p());
1171  for(Int i=0; i<val.m(); i++)
1172  for(Int j=0; j<val.n(); j++)
1173  for(Int k=0; k<val.p(); k++)
1174  val(i,j,k) += ext(i,j,k);
1175  return 0;
1176  }
1177 
1178  //-------------------
1179  //NumVec
1180  template<class T>
1181  Int inline serialize(const NumVec<T>& val, std::ostream& os, const std::vector<Int>& mask)
1182  {
1183  Int m = val.m();
1184  os.write((char*)&m, sizeof(Int));
1185  for(Int i=0; i<m; i++)
1186  serialize(val(i), os, mask);
1187  return 0;
1188  }
1189 
1190  template<class T>
1191  Int inline deserialize(NumVec<T>& val, std::istream& is, const std::vector<Int>& mask)
1192  {
1193  Int m;
1194  is.read((char*)&m, sizeof(Int));
1195  val.Resize(m);
1196  for(Int i=0; i<m; i++)
1197  deserialize(val(i), is, mask);
1198  return 0;
1199  }
1200 
1201  template<class T>
1202  Int inline combine(NumVec<T>& val, NumVec<T>& ext)
1203  {
1204  throw std::logic_error( "Combine operation not implemented." );
1205  return 0;
1206  }
1207 
1208  //-------------------
1209  //NumMat
1210  template<class T>
1211  Int inline serialize(const NumMat<T>& val, std::ostream& os, const std::vector<Int>& mask)
1212  {
1213  Int m = val.m();
1214  Int n = val.n();
1215  os.write((char*)&m, sizeof(Int));
1216  os.write((char*)&n, sizeof(Int));
1217  for(Int j=0; j<n; j++)
1218  for(Int i=0; i<m; i++)
1219  serialize(val(i,j), os, mask);
1220  return 0;
1221  }
1222  template<class T>
1223  Int inline deserialize(NumMat<T>& val, std::istream& is, const std::vector<Int>& mask)
1224  {
1225  Int m;
1226  Int n;
1227  is.read((char*)&m, sizeof(Int));
1228  is.read((char*)&n, sizeof(Int));
1229  val.Resize(m,n);
1230  for(Int j=0; j<n; j++)
1231  for(Int i=0; i<m; i++)
1232  deserialize(val(i,j), is, mask);
1233  return 0;
1234  }
1235 
1236  template<class T>
1237  Int inline combine(NumMat<T>& val, NumMat<T>& ext)
1238  {
1239  throw std::logic_error( "Combine operation not implemented." );
1240  return 0;
1241  }
1242 
1243 
1244  //-------------------
1245  //NumTns
1246  template<class T>
1247  Int inline serialize(const NumTns<T>& val, std::ostream& os, const std::vector<Int>& mask)
1248  {
1249  Int m = val.m();
1250  Int n = val.n();
1251  Int p = val.p();
1252  os.write((char*)&m, sizeof(Int));
1253  os.write((char*)&n, sizeof(Int));
1254  os.write((char*)&p, sizeof(Int));
1255  for(Int k=0; k<p; k++)
1256  for(Int j=0; j<n; j++)
1257  for(Int i=0; i<m; i++)
1258  serialize(val(i,j,k), os, mask);
1259  return 0;
1260  }
1261 
1262  template<class T>
1263  Int inline deserialize(NumTns<T>& val, std::istream& is, const std::vector<Int>& mask)
1264  {
1265  Int m;
1266  Int n;
1267  Int p;
1268  is.read((char*)&m, sizeof(Int));
1269  is.read((char*)&n, sizeof(Int));
1270  is.read((char*)&p, sizeof(Int));
1271  val.Resize(m,n,p);
1272  for(Int k=0; k<p; k++)
1273  for(Int j=0; j<n; j++)
1274  for(Int i=0; i<m; i++)
1275  deserialize(val(i,j,k), is, mask);
1276  return 0;
1277  }
1278 
1279  template<class T>
1280  Int inline combine(NumTns<T>& val, NumTns<T>& ext)
1281  {
1282  throw std::logic_error( "Combine operation not implemented." );
1283  return 0;
1284  }
1285 
1286  //-------------------
1287  //DistSparseMatrix
1288  template<class T>
1289  Int inline serialize(const DistSparseMatrix<T>& val, std::ostream& os, const std::vector<Int>& mask)
1290  {
1291  serialize( val.size, os, mask );
1292  serialize( val.nnz, os, mask );
1293  serialize( val.nnzLocal, os, mask );
1294  serialize( val.colptrLocal, os, mask );
1295  serialize( val.rowindLocal, os, mask );
1296  serialize( val.nzvalLocal, os, mask );
1297  // No need to serialize the communicator
1298  return 0;
1299  }
1300 
1301  template<class T>
1302  Int inline deserialize(DistSparseMatrix<T>& val, std::istream& is, const std::vector<Int>& mask)
1303  {
1304  deserialize( val.size, is, mask );
1305  deserialize( val.nnz, is, mask );
1306  deserialize( val.nnzLocal, is, mask );
1307  deserialize( val.colptrLocal, is, mask );
1308  deserialize( val.rowindLocal, is, mask );
1309  deserialize( val.nzvalLocal, is, mask );
1310  // No need to deserialize the communicator
1311  return 0;
1312  }
1313 
1314  template<class T>
1315  Int inline combine(DistSparseMatrix<T>& val, DistSparseMatrix<T>& ext)
1316  {
1317  throw std::logic_error( "Combine operation not implemented." );
1318  return 0;
1319  }
1320 
1321 
1322  // *********************************************************************
1323  // Parallel IO functions
1324  // *********************************************************************
1325 
1326  Int SeparateRead(std::string name, std::istringstream& is);
1327 
1328  Int SeparateWrite(std::string name, std::ostringstream& os);
1329 
1330  Int SeparateWriteAscii(std::string name, std::ostringstream& os);
1331 
1332  Int SharedRead(std::string name, std::istringstream& is);
1333 
1334  Int SharedWrite(std::string name, std::ostringstream& os);
1335 
1336  // *********************************************************************
1337  // Ej
1338  // *********************************************************************
1339  inline void IdentityCol( Int col, NumVec<Real>& vec )
1340  {
1341  for(Int i=0; i<std::min(col,vec.m()); i++)
1342  vec(i) = 0.0;
1343 
1344  if(col<vec.m())
1345  vec(col) = 1.0;
1346 
1347 
1348  for(Int i=col+1; i<vec.m(); i++)
1349  vec(i) = 0.0;
1350  }
1351 
1352  inline void IdentityCol( Int col, NumVec<Complex>& vec )
1353  {
1354  for(Int i=0; i<std::min(col,vec.m()); i++)
1355  vec(i) = Complex(0.0,0.0);
1356 
1357  if(col<vec.m())
1358  vec(col) = Complex(1.0,0.0);
1359 
1360 
1361  for(Int i=col+1; i<vec.m(); i++)
1362  vec(i) = Complex(0.0,0.0);
1363  }
1364 
1365 
1366 
1367  inline void IdentityCol( IntNumVec & cols, NumMat<Real>& mat )
1368  {
1369  for(Int j=0;j<cols.m();j++){
1370  Int col= cols(j);
1371  for(Int i=0; i<std::min(col,mat.m()); i++)
1372  mat(i,j) = 0.0;
1373 
1374  if(col<mat.m())
1375  mat(col,j) = 1.0;
1376 
1377 
1378  for(Int i=col+1; i<mat.m(); i++)
1379  mat(i,j) = 0.0;
1380  }
1381  }
1382 
1383 
1384  inline void IdentityCol( IntNumVec & cols, NumMat<Complex>& mat )
1385  {
1386  for(Int j=0;j<cols.m();j++){
1387  Int col= cols(j);
1388  for(Int i=0; i<std::min(col,mat.m()); i++)
1389  mat(i,j) = Complex(0.0,0.0);
1390 
1391  if(col<mat.m())
1392  mat(col,j) = Complex(1.0,0.0);
1393 
1394 
1395  for(Int i=col+1; i<mat.m(); i++)
1396  mat(i,j) = Complex(0.0,0.0);
1397  }
1398  }
1399 
1400 
1401 
1402 
1403  // *********************************************************************
1404  // Random numbers
1405  // *********************************************************************
1406  inline void SetRandomSeed(long int seed){
1407  srand48(seed);
1408  }
1409 
1410  inline Real UniformRandom(){
1411  return (Real)drand48();
1412  }
1413 
1414  inline void UniformRandom( NumVec<Real>& vec )
1415  {
1416  for(Int i=0; i<vec.m(); i++)
1417  vec(i) = UniformRandom();
1418  }
1419 
1420  inline void UniformRandom( NumVec<Complex>& vec )
1421  {
1422  for(Int i=0; i<vec.m(); i++)
1423  vec(i) = Complex(UniformRandom(), UniformRandom());
1424  }
1425 
1426  inline void UniformRandom( NumMat<Real>& M )
1427  {
1428  Real *ptr = M.Data();
1429  for(Int i=0; i < M.m() * M.n(); i++)
1430  *(ptr++) = UniformRandom();
1431  }
1432 
1433  inline void UniformRandom( NumMat<Complex>& M )
1434  {
1435  Complex *ptr = M.Data();
1436  for(Int i=0; i < M.m() * M.n(); i++)
1437  *(ptr++) = Complex(UniformRandom(), UniformRandom());
1438  }
1439 
1440 
1441  inline void UniformRandom( NumTns<Real>& T )
1442  {
1443  Real *ptr = T.Data();
1444  for(Int i=0; i < T.m() * T.n() * T.p(); i++)
1445  *(ptr++) = UniformRandom();
1446  }
1447 
1448  inline void UniformRandom( NumTns<Complex>& T )
1449  {
1450  Complex *ptr = T.Data();
1451  for(Int i=0; i < T.m() * T.n() * T.p(); i++)
1452  *(ptr++) = Complex(UniformRandom(), UniformRandom());
1453  }
1454 
1455  // *********************************************************************
1456  // Timing
1457  // *********************************************************************
1458  inline void GetTime(Real& t){
1459  t = MPI_Wtime();
1460  }
1461 
1462  // *********************************************************************
1463  // Comparator
1464  // *********************************************************************
1465 
1466  // Real
1467  inline bool PairLtComparator( const std::pair<Real, Int>& l,
1468  const std::pair<Real, Int>& r ){
1469  return l.first < r.first;
1470  }
1471 
1472  inline bool PairGtComparator( const std::pair<Real, Int>& l,
1473  const std::pair<Real, Int>& r ){
1474  return l.first > r.first;
1475  }
1476 
1477  // For sorting with indices
1478  // Example usage:
1479  // std::sort(val.begin(), val.end(), IndexComp<std::vector<int>&>(indices));
1480  template<class T>
1481  struct IndexComp {
1482  private:
1483  const T indices_;
1484  public:
1485  IndexComp (const T indices) : indices_(indices) {}
1486  bool operator()(const size_t a, const size_t b) const
1487  { return indices_[a] < indices_[b]; }
1488  };
1489 
1490 
1491 
1492 
1493  // *********************************************************************
1494  // Sparse Matrix
1495  // *********************************************************************
1496 
1497  // TODO Complex format
1498  void ReadSparseMatrix ( const char* filename, SparseMatrix<Real>& spmat );
1499 
1500  void ReadDistSparseMatrix( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1501 
1502  void ParaReadDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1503 
1504  void ParaWriteDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1505 
1506  void ReadDistSparseMatrixFormatted( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1507 
1508  template <class F1, class F2>
1509  void
1510  CopyPattern ( const SparseMatrix<F1>& A, SparseMatrix<F2>& B )
1511  {
1512 #ifndef _RELEASE_
1513  PushCallStack("CopyPattern");
1514 #endif
1515  B.size = A.size;
1516  B.nnz = A.nnz;
1517  B.colptr = A.colptr;
1518  B.rowind = A.rowind;
1519  B.nzval.Resize( A.nnz );
1520 #ifndef _RELEASE_
1521  PopCallStack();
1522 #endif
1523  return ;
1524  } // ----- end of template function CopyPattern -----
1525 
1526 
1527  // TODO Real format
1528  template<typename T>
1529  void
1530  GetDiagonal ( const DistSparseMatrix<T>& A,
1531  NumVec<T>& diag );
1532 
1533 
1534  // Functions for DistSparseMatrix
1535 
1536  template <class F1, class F2>
1537  void
1538  CopyPattern ( const DistSparseMatrix<F1>& A, DistSparseMatrix<F2>& B )
1539  {
1540 #ifndef _RELEASE_
1541  PushCallStack("CopyPattern");
1542 #endif
1543  B.size = A.size;
1544  B.nnz = A.nnz;
1545  B.nnzLocal = A.nnzLocal;
1546  B.colptrLocal = A.colptrLocal;
1547  B.rowindLocal = A.rowindLocal;
1548  B.nzvalLocal.Resize( A.nnzLocal );
1549  B.comm = A.comm;
1550 #ifndef _RELEASE_
1551  PopCallStack();
1552 #endif
1553  return ;
1554  } // ----- end of template function CopyPattern -----
1555 
1560 
1573  template <class F>
1574  void
1575  DistSparseMatMultGlobalVec(
1576  const F alpha,
1577  const DistSparseMatrix<F>& AMat,
1578  const NumVec<F>& B,
1579  const F beta,
1580  NumVec<F>& C )
1581  {
1582 #ifndef _RELEASE_
1583  PushCallStack("DistSparseMatMultGlobalVec");
1584 #endif
1585  if( B.m() != AMat.size ){
1586  std::ostringstream msg;
1587  msg << std::endl
1588  << "The size of the matrix A is " << AMat.size << std::endl
1589  << "The size of the vector B is " << B.m() << std::endl
1590  << "and they must agree with each other.";
1591  throw std::runtime_error( msg.str().c_str() );
1592  }
1593  if( C.m() != AMat.size ){
1594  std::ostringstream msg;
1595  msg << std::endl
1596  << "The size of the matrix A is " << AMat.size << std::endl
1597  << "The size of the vector C is " << C.m() << std::endl
1598  << "and they must agree with each other.";
1599  throw std::runtime_error( msg.str().c_str() );
1600  }
1601 
1602  MPI_Comm comm = AMat.comm;
1603 
1604  Int mpirank, mpisize;
1605  MPI_Comm_rank( comm, &mpirank );
1606  MPI_Comm_size( comm, &mpisize );
1607 
1608  // Initiate a new vector for saving DeltaC = alpha * A * B
1609  // This is somewhat inefficient implementation but this should not
1610  // matter at this moment, since the global distribution of a dense
1611  // vector is not important anyway.
1612  NumVec<F> DeltaClocal( AMat.size );
1613  NumVec<F> DeltaC( AMat.size );
1614  SetValue( DeltaClocal, 0.0 );
1615  SetValue( DeltaC, 0.0 );
1616 
1617  Int numColLocal = AMat.colptrLocal.m() - 1;
1618  Int numColLocalFirst = AMat.size / mpisize;
1619  Int firstCol = mpirank * numColLocalFirst;
1620 
1621  // NOTE: DistSparseMatrix use 1-based indicies, while access to
1622  // vectors uses 0-based indicies.
1623  for( Int j = 0; j < numColLocal; j++ ){
1624  Int jcol = firstCol + j;
1625  for( Int i = AMat.colptrLocal(j)-1;
1626  i < AMat.colptrLocal(j+1)-1; i++ ){
1627  Int irow = AMat.rowindLocal(i) - 1;
1628  DeltaClocal( irow ) += AMat.nzvalLocal(i) * B(jcol);
1629  }
1630  } // for (j)
1631 
1632  mpi::Allreduce( DeltaClocal.Data(), DeltaC.Data(),
1633  AMat.size, MPI_SUM, comm );
1634 
1635  for( Int i = 0; i < AMat.size; i++ ){
1636  C(i) = beta * C(i) + alpha * DeltaC(i);
1637  }
1638 
1639 #ifndef _RELEASE_
1640  PopCallStack();
1641 #endif
1642  return ;
1643  } // ----- end of template function DistSparseMatMultGlobalVec -----
1644 
1645  // *********************************************************************
1646  // Other numerical routines
1647  // *********************************************************************
1648 
1649  // Interpolation
1658  void
1659  LinearInterpolation (
1660  const std::vector<Real>& x,
1661  const std::vector<Real>& y,
1662  const std::vector<Real>& xx,
1663  std::vector<Real>& yy );
1664 
1665  // Root finding
1666 
1680  inline Real
1681  MonotoneRootFinding (
1682  const std::vector<Real>& x,
1683  const std::vector<Real>& y,
1684  Real val )
1685  {
1686 #ifndef _RELEASE_
1687  PushCallStack("MonotoneRootFinding");
1688 #endif
1689  Int numX = x.size();
1690 
1691  if( val <= y[0] || val >= y[numX-1] ){
1692  std::ostringstream msg;
1693  msg
1694  << "The root finding procedure cannot find the solution for y(x)="
1695  << val << std::endl << "here [min(y) max(y)] = [" << y[0] << ", "
1696  << y[numX-1] << "]" << std::endl;
1697  throw std::runtime_error( msg.str().c_str() );
1698  }
1699 
1700  std::vector<Real>::const_iterator vi = std::lower_bound(
1701  y.begin(), y.end(), val );
1702 
1703  Int idx = vi - y.begin();
1704 
1705  Real root;
1706  if( y[idx] == y[idx-1] ){
1707  root = x[idx-1];
1708  }
1709  else{
1710  // Linear interpolation
1711  root = x[idx-1] + ( x[idx] - x[idx-1] ) / ( y[idx] - y[idx-1] )
1712  * ( val - y[idx-1] );
1713  }
1714 
1715 #ifndef _RELEASE_
1716  PopCallStack();
1717 #endif
1718 
1719  return root;
1720  } // ----- end of function MonotoneRootFinding -----
1721 
1722 
1723 
1724 } // namespace PEXSI
1725 
1726 #include "utility_impl.hpp"
1727 
1728 #endif // _PEXSI_UTILITY_HPP_
PEXSI::NumVec
Numerical vector.
Definition: NumVec.hpp:60
environment.hpp
Environmental variables.
PEXSI::DistSparseMatrix::colptrLocal
IntNumVec colptrLocal
Dimension numColLocal + 1, storing the pointers to the nonzero row indices and nonzero values in rowp...
Definition: sparse_matrix.hpp:108
PEXSI::SetValue
void SetValue(NumMat< F > &M, F val)
SetValue sets a numerical matrix to a constant val.
Definition: NumMat_impl.hpp:137
utility_impl.hpp
Implementation of utility subroutines.
PEXSI::DistSparseMatrix::nzvalLocal
NumVec< F > nzvalLocal
Dimension nnzLocal, storing the nonzero values.
Definition: sparse_matrix.hpp:116
PEXSI::SparseMatrix
SparseMatrix describes a sequential sparse matrix saved in compressed sparse column format...
Definition: sparse_matrix.hpp:66
PEXSI::DistSparseMatrix::rowindLocal
IntNumVec rowindLocal
Dimension nnzLocal, storing the nonzero row indices. The indices are 1-based (FORTRAN-convention), i.e. the first row index is 1.
Definition: sparse_matrix.hpp:113
sparse_matrix.hpp
Sparse matrix and Distributed sparse matrix in compressed column format.
NumMat.hpp
Numerical matrix.
mpi_interf.hpp
Interface with MPI to facilitate communication.
PEXSI::DistSparseMatrix::size
Int size
Matrix dimension.
Definition: sparse_matrix.hpp:93
PEXSI::DistSparseMatMultGlobalVec
void DistSparseMatMultGlobalVec(const F alpha, const DistSparseMatrix< F > &AMat, const NumVec< F > &B, const F beta, NumVec< F > &C)
Multiply a DistSparseMatrix with a vector that is distributed across all processors participating the...
Definition: utility.hpp:1575
NumTns.hpp
Numerical tensor.
PEXSI::DistSparseMatrix::comm
MPI_Comm comm
MPI communicator.
Definition: sparse_matrix.hpp:119
tinyvec.hpp
Tiny vectors of dimension 3.
NumVec.hpp
Numerical vector.
PEXSI::MonotoneRootFinding
Real MonotoneRootFinding(const std::vector< Real > &x, const std::vector< Real > &y, Real val)
Root finding for monotonic non-decreasing functions.
Definition: utility.hpp:1681
PEXSI::IndexComp
Definition: utility.hpp:1481
PEXSI::LinearInterpolation
void LinearInterpolation(const std::vector< Real > &x, const std::vector< Real > &y, const std::vector< Real > &xx, std::vector< Real > &yy)
Linear interpolates from (x,y) to (xx,yy)
Definition: utility.cpp:209
Generated on Sat May 16 2015 11:38:16 for PEXSI by   doxygen 1.8.6