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 "pexsi/environment.hpp"
51 #include "pexsi/tinyvec.hpp"
52 #include "pexsi/NumVec.hpp"
53 #include "pexsi/NumMat.hpp"
54 #include "pexsi/NumTns.hpp"
55 #include "pexsi/sparse_matrix.hpp"
56 #include "pexsi/mpi_interf.hpp"
57 
58 #ifdef WITH_SYMPACK
59 #include <sympack.hpp>
60 #endif
61 
62 namespace PEXSI{
63 
64 // *********************************************************************
65 // Define constants
66 // *********************************************************************
67 
68 // Write format control parameters
69 const int LENGTH_VAR_NAME = 8;
70 const int LENGTH_DBL_DATA = 16;
71 const int LENGTH_INT_DATA = 5;
72 const int LENGTH_VAR_UNIT = 6;
73 const int LENGTH_DBL_PREC = 8;
74 const int LENGTH_VAR_DATA = 16;
75 
76 
77 // standard case for most serialization/deserialization process.
78 const std::vector<Int> NO_MASK(1);
79 
80 // *********************************************************************
81 // Formatted output stream
82 // *********************************************************************
83 
84 // Bool is NOT defined due to ambiguity with Int
85 
86 inline Int PrintBlock(std::ostream &os, const std::string name){
87 
88  os << std::endl<< "*********************************************************************" << std::endl;
89  os << name << std::endl;
90  os << "*********************************************************************" << std::endl << std::endl;
91  return 0;
92 }
93 
94 // String
95 inline Int Print(std::ostream &os, const std::string name) {
96  os << std::setiosflags(std::ios::left) << name << std::endl;
97  return 0;
98 };
99 
100 inline Int Print(std::ostream &os, const char* name) {
101  os << std::setiosflags(std::ios::left) << std::string(name) << std::endl;
102  return 0;
103 };
104 
105 inline Int Print(std::ostream &os, const std::string name, std::string val) {
106  os << std::setiosflags(std::ios::left)
107  << std::setw(LENGTH_VAR_NAME) << name
108  << std::setw(LENGTH_VAR_DATA) << val
109  << std::endl;
110  return 0;
111 };
112 
113 inline Int Print(std::ostream &os, const std::string name, const char* val) {
114  os << std::setiosflags(std::ios::left)
115  << std::setw(LENGTH_VAR_NAME) << name
116  << std::setw(LENGTH_VAR_DATA) << std::string(val)
117  << std::endl;
118  return 0;
119 };
120 
121 
122 // Real
123 
124 // one real number
125 
126 inline Int Print(std::ostream &os, const std::string name, Real val) {
127  os << std::setiosflags(std::ios::left)
128  << std::setw(LENGTH_VAR_NAME) << name
129  << std::setiosflags(std::ios::scientific)
130  << std::setiosflags(std::ios::showpos)
131  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
132  << std::resetiosflags(std::ios::scientific)
133  << std::resetiosflags(std::ios::showpos)
134  << std::endl;
135  return 0;
136 };
137 
138 inline Int Print(std::ostream &os, const char* name, Real val) {
139  os << std::setiosflags(std::ios::left)
140  << std::setw(LENGTH_VAR_NAME) << std::string(name)
141  << std::setiosflags(std::ios::scientific)
142  << std::setiosflags(std::ios::showpos)
143  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
144  << std::resetiosflags(std::ios::scientific)
145  << std::resetiosflags(std::ios::showpos)
146  << std::endl;
147  return 0;
148 };
149 
150 
151 inline Int Print(std::ostream &os, const std::string name, Real val, const std::string unit) {
152  os << std::setiosflags(std::ios::left)
153  << std::setw(LENGTH_VAR_NAME) << name
154  << std::setiosflags(std::ios::scientific)
155  << std::setiosflags(std::ios::showpos)
156  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
157  << std::resetiosflags(std::ios::scientific)
158  << std::resetiosflags(std::ios::showpos)
159  << std::setw(LENGTH_VAR_UNIT) << unit
160  << std::endl;
161  return 0;
162 };
163 
164 inline Int Print(std::ostream &os, const char *name, Real val, const char *unit) {
165  os << std::setiosflags(std::ios::left)
166  << std::setw(LENGTH_VAR_NAME) << std::string(name)
167  << std::setiosflags(std::ios::scientific)
168  << std::setiosflags(std::ios::showpos)
169  << std::setw(LENGTH_DBL_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val
170  << std::resetiosflags(std::ios::scientific)
171  << std::resetiosflags(std::ios::showpos)
172  << std::setw(LENGTH_VAR_UNIT) << std::string(unit)
173  << std::endl;
174  return 0;
175 };
176 
177 // Two real numbers
178 inline Int Print(std::ostream &os, const std::string name1, Real val1, const std::string unit1,
179  const std::string name2, Real val2, const std::string unit2) {
180  os << std::setiosflags(std::ios::left)
181  << std::setw(LENGTH_VAR_NAME) << name1
182  << std::setiosflags(std::ios::scientific)
183  << std::setiosflags(std::ios::showpos)
184  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val1
185  << std::setw(LENGTH_VAR_UNIT) << unit1
186  << std::setw(LENGTH_VAR_NAME) << name2
187  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
188  << std::resetiosflags(std::ios::scientific)
189  << std::resetiosflags(std::ios::showpos)
190  << std::setw(LENGTH_VAR_UNIT) << unit2
191  << std::endl;
192  return 0;
193 };
194 
195 inline Int Print(std::ostream &os, const char *name1, Real val1, const char *unit1,
196  char *name2, Real val2, char *unit2) {
197  os << std::setiosflags(std::ios::left)
198  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
199  << std::setiosflags(std::ios::scientific)
200  << std::setiosflags(std::ios::showpos)
201  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val1
202  << std::setw(LENGTH_VAR_UNIT) << std::string(unit1)
203  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
204  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
205  << std::resetiosflags(std::ios::scientific)
206  << std::resetiosflags(std::ios::showpos)
207  << std::setw(LENGTH_VAR_UNIT) << std::string(unit2)
208  << std::endl;
209  return 0;
210 };
211 
212 // Int and Real
213 inline Int Print(std::ostream &os, const std::string name1, Int val1, const std::string unit1,
214  const std::string name2, Real val2, const std::string unit2) {
215  os << std::setiosflags(std::ios::left)
216  << std::setw(LENGTH_VAR_NAME) << name1
217  << std::setw(LENGTH_INT_DATA) << val1
218  << std::setw(LENGTH_VAR_UNIT) << unit1
219  << std::setw(LENGTH_VAR_NAME) << name2
220  << std::setiosflags(std::ios::scientific)
221  << std::setiosflags(std::ios::showpos)
222  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
223  << std::resetiosflags(std::ios::scientific)
224  << std::resetiosflags(std::ios::showpos)
225  << std::setw(LENGTH_VAR_UNIT) << unit2
226  << std::endl;
227  return 0;
228 };
229 
230 inline Int Print(std::ostream &os, const char *name1, Int val1, const char *unit1,
231  char *name2, Real val2, char *unit2) {
232  os << std::setiosflags(std::ios::left)
233  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
234  << std::setw(LENGTH_INT_DATA) << val1
235  << std::setw(LENGTH_VAR_UNIT) << std::string(unit1)
236  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
237  << std::setiosflags(std::ios::scientific)
238  << std::setiosflags(std::ios::showpos)
239  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
240  << std::resetiosflags(std::ios::scientific)
241  << std::resetiosflags(std::ios::showpos)
242  << std::setw(LENGTH_VAR_UNIT) << std::string(unit2)
243  << std::endl;
244  return 0;
245 };
246 
247 inline Int Print(std::ostream &os,
248  const char *name1, Int val1,
249  const char *name2, Real val2,
250  char *name3, Real val3) {
251  os << std::setiosflags(std::ios::left)
252  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
253  << std::setw(LENGTH_INT_DATA) << val1
254  << std::setiosflags(std::ios::scientific)
255  << std::setiosflags(std::ios::showpos)
256  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
257  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
258  << std::setw(LENGTH_VAR_NAME) << std::string(name3)
259  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val3
260  << std::resetiosflags(std::ios::scientific)
261  << std::resetiosflags(std::ios::showpos)
262  << std::endl;
263  return 0;
264 };
265 
266 
267 inline Int Print(std::ostream &os,
268  const char *name1, Int val1,
269  const char *name2, Real val2,
270  const char *name3, Real val3,
271  const char *name4, Real val4 ) {
272  os << std::setiosflags(std::ios::left)
273  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
274  << std::setw(LENGTH_INT_DATA) << val1
275  << std::setiosflags(std::ios::scientific)
276  << std::setiosflags(std::ios::showpos)
277  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
278  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val2
279  << std::setw(LENGTH_VAR_NAME) << std::string(name3)
280  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val3
281  << std::setw(LENGTH_VAR_NAME) << std::string(name4)
282  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC)<< val4
283  << std::resetiosflags(std::ios::scientific)
284  << std::resetiosflags(std::ios::showpos)
285  << std::endl;
286  return 0;
287 };
288 
289 // Int
290 
291 // one integer number
292 inline Int Print(std::ostream &os, std::string name, Int val) {
293  os << std::setiosflags(std::ios::left)
294  << std::setw(LENGTH_VAR_NAME) << name
295  << std::setw(LENGTH_VAR_DATA) << val
296  << std::endl;
297  return 0;
298 };
299 
300 
301 inline Int Print(std::ostream &os, const char *name, Int val) {
302  os << std::setiosflags(std::ios::left)
303  << std::setw(LENGTH_VAR_NAME) << std::string(name)
304  << std::setw(LENGTH_VAR_DATA) << val
305  << std::endl;
306  return 0;
307 };
308 
309 
310 inline Int Print(std::ostream &os, const std::string name, Int val, const std::string unit) {
311  os << std::setiosflags(std::ios::left)
312  << std::setw(LENGTH_VAR_NAME) << name
313  << std::setw(LENGTH_VAR_DATA) << val
314  << std::setw(LENGTH_VAR_UNIT) << unit
315  << std::endl;
316  return 0;
317 };
318 
319 
320 inline Int Print(std::ostream &os, const char* name, Int val, const std::string unit) {
321  os << std::setiosflags(std::ios::left)
322  << std::setw(LENGTH_VAR_NAME) << std::string(name)
323  << std::setw(LENGTH_VAR_DATA) << val
324  << std::setw(LENGTH_VAR_UNIT) << unit
325  << std::endl;
326  return 0;
327 };
328 
329 
330 
331 // two integer numbers
332 inline Int Print(std::ostream &os, const std::string name1, Int val1, const std::string unit1,
333  const std::string name2, Int val2, const std::string unit2) {
334  os << std::setiosflags(std::ios::left)
335  << std::setw(LENGTH_VAR_NAME) << name1
336  << std::setw(LENGTH_VAR_DATA) << val1
337  << std::setw(LENGTH_VAR_UNIT) << unit1
338  << std::setw(LENGTH_VAR_NAME) << name2
339  << std::setw(LENGTH_VAR_DATA) << val2
340  << std::setw(LENGTH_VAR_UNIT) << unit2
341  << std::endl;
342  return 0;
343 };
344 
345 inline Int Print(std::ostream &os, const char *name1, Int val1, const char *unit1,
346  char *name2, Int val2, char *unit2) {
347  os << std::setiosflags(std::ios::left)
348  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
349  << std::setw(LENGTH_VAR_DATA) << val1
350  << std::setw(LENGTH_VAR_UNIT) << std::string(unit1)
351  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
352  << std::setw(LENGTH_VAR_DATA) << val2
353  << std::setw(LENGTH_VAR_UNIT) << std::string(unit2)
354  << std::endl;
355  return 0;
356 };
357 
358 // Bool
359 
360 // one boolean number
361 inline Int Print(std::ostream &os, const std::string name, bool val) {
362  os << std::setiosflags(std::ios::left)
363  << std::setw(LENGTH_VAR_NAME) << name;
364  if( val == true )
365  os << std::setw(LENGTH_VAR_NAME) << "true" << std::endl;
366  else
367  os << std::setw(LENGTH_VAR_NAME) << "false" << std::endl;
368  return 0;
369 };
370 
371 
372 inline Int Print(std::ostream &os, const char* name, bool val) {
373  os << std::setiosflags(std::ios::left)
374  << std::setw(LENGTH_VAR_NAME) << std::string(name);
375  if( val == true )
376  os << std::setw(LENGTH_VAR_NAME) << "true" << std::endl;
377  else
378  os << std::setw(LENGTH_VAR_NAME) << "false" << std::endl;
379  return 0;
380 };
381 
382 
383 // Index 3 and Point 3
384 inline Int Print(std::ostream &os,
385  const char *name1, Index3 val ) {
386  os << std::setiosflags(std::ios::left)
387  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
388  << std::setw(LENGTH_VAR_DATA) << val[0]
389  << std::setw(LENGTH_VAR_DATA) << val[1]
390  << std::setw(LENGTH_VAR_DATA) << val[2]
391  << std::endl;
392  return 0;
393 };
394 
395 inline Int Print(std::ostream &os,
396  const char *name1, Point3 val ) {
397  os << std::setiosflags(std::ios::left)
398  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
399  << std::setiosflags(std::ios::scientific)
400  << std::setiosflags(std::ios::showpos)
401  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) << val[0]
402  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) << val[1]
403  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) << val[2]
404  << std::resetiosflags(std::ios::scientific)
405  << std::resetiosflags(std::ios::showpos)
406  << std::endl;
407  return 0;
408 };
409 
410 inline Int Print(std::ostream &os,
411  const char *name1, Int val1,
412  const char *name2, Point3 val ) {
413  os << std::setiosflags(std::ios::left)
414  << std::setw(LENGTH_VAR_NAME) << std::string(name1)
415  << std::setw(LENGTH_INT_DATA) << val1
416  << std::setw(LENGTH_VAR_NAME) << std::string(name2)
417  << std::setiosflags(std::ios::scientific)
418  << std::setiosflags(std::ios::showpos)
419  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) <<val[0]
420  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) <<val[1]
421  << std::setw(LENGTH_VAR_DATA) << std::setprecision(LENGTH_DBL_PREC) <<val[2]
422  << std::resetiosflags(std::ios::scientific)
423  << std::resetiosflags(std::ios::showpos)
424  << std::endl;
425  return 0;
426 };
427 
428 // *********************************************************************
429 // Overload << and >> operators for basic data types
430 // *********************************************************************
431 
432 // std::vector
433 template <class F> inline std::ostream& operator<<( std::ostream& os, const std::vector<F>& vec)
434 {
435  os<<vec.size()<<std::endl;
436  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
437  for(Int i=0; i<vec.size(); i++)
438  os<<" "<<vec[i];
439  os<<std::endl;
440  return os;
441 }
442 
443 // NumVec
444 template <class F> inline std::ostream& operator<<( std::ostream& os, const NumVec<F>& vec)
445 {
446  os<<vec.m()<<std::endl;
447  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
448  for(Int i=0; i<vec.m(); i++)
449  os<<" "<<vec(i);
450  os<<std::endl;
451  return os;
452 }
453 
454 //template <class F> inline std::istream& operator>>( std::istream& is, NumVec<F>& vec)
455 //{
456 // Int m; is>>m; vec.resize(m);
457 // for(Int i=0; i<vec.m(); i++)
458 // is >> vec(i);
459 // return is;
460 //}
461 
462 // NumMat
463 template <class F> inline std::ostream& operator<<( std::ostream& os, const NumMat<F>& mat)
464 {
465  os<<mat.m()<<" "<<mat.n()<<std::endl;
466  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
467  os << std::setprecision(16);
468  for(Int i=0; i<mat.m(); i++) {
469  for(Int j=0; j<mat.n(); j++)
470  os<<" "<<mat(i,j);
471  os<<std::endl;
472  }
473  return os;
474 }
475 
476 // NumTns
477 template <class F> inline std::ostream& operator<<( std::ostream& os, const NumTns<F>& tns)
478 {
479  os<<tns.m()<<" "<<tns.n()<<" "<<tns.p()<<std::endl;
480  os.setf(std::ios_base::scientific, std::ios_base::floatfield);
481  for(Int i=0; i<tns.m(); i++) {
482  for(Int j=0; j<tns.n(); j++) {
483  for(Int k=0; k<tns.p(); k++) {
484  os<<" "<<tns(i,j,k);
485  }
486  os<<std::endl;
487  }
488  os<<std::endl;
489  }
490  return os;
491 }
492 
493 // *********************************************************************
494 // serialize/deserialize for basic types
495 // More specific serialize/deserialize will be defined in individual
496 // class files
497 // *********************************************************************
498 
499 
500 
501 
502 //template
503 template <typename T>
504 inline Int serialize(const T& val, std::ostream& os, const std::vector<Int>& mask)
505 {
506  os.write((char*)&val, sizeof(T));
507  return 0;
508 }
509 
510 template <typename T>
511 inline Int deserialize(T& val, std::istream& is, const std::vector<Int>& mask)
512 {
513  is.read((char*)&val, sizeof(T));
514  return 0;
515 }
516 
517 
518 
519 //bool
520 inline Int serialize(const bool& val, std::ostream& os, const std::vector<Int>& mask)
521 {
522  os.write((char*)&val, sizeof(bool));
523  return 0;
524 }
525 
526 inline Int deserialize(bool& val, std::istream& is, const std::vector<Int>& mask)
527 {
528  is.read((char*)&val, sizeof(bool));
529  return 0;
530 }
531 
532 //char
533 inline Int serialize(const char& val, std::ostream& os, const std::vector<Int>& mask)
534 {
535  os.write((char*)&val, sizeof(char));
536  return 0;
537 }
538 
539 inline Int deserialize(char& val, std::istream& is, const std::vector<Int>& mask)
540 {
541  is.read((char*)&val, sizeof(char));
542  return 0;
543 }
544 
545 inline Int combine(char& val, char& ext)
546 {
547  ErrorHandling( "Combine operation not implemented." );
548 }
549 
550 //-------------------
551 //Int
552 inline Int serialize(const Int& val, std::ostream& os, const std::vector<Int>& mask)
553 {
554  os.write((char*)&val, sizeof(Int));
555  return 0;
556 }
557 
558 inline Int deserialize(Int& val, std::istream& is, const std::vector<Int>& mask)
559 {
560  is.read((char*)&val, sizeof(Int));
561  return 0;
562 }
563 
564 inline Int combine(Int& val, Int& ext)
565 {
566  val += ext;
567  return 0;
568 }
569 
570 //-------------------
571 //LongInt
572 inline Int serialize(const LongInt& val, std::ostream& os, const std::vector<Int>& mask)
573 {
574  os.write((char*)&val, sizeof(LongInt));
575  return 0;
576 }
577 
578 inline Int deserialize(LongInt& val, std::istream& is, const std::vector<Int>& mask)
579 {
580  is.read((char*)&val, sizeof(LongInt));
581  return 0;
582 }
583 
584 inline Int combine(LongInt& val, LongInt& ext)
585 {
586  val += ext;
587  return 0;
588 }
589 
590 
591 //-------------------
592 //Real
593 inline Int serialize(const Real& val, std::ostream& os, const std::vector<Int>& mask)
594 {
595  os.write((char*)&val, sizeof(Real));
596  return 0;
597 }
598 
599 inline Int deserialize(Real& val, std::istream& is, const std::vector<Int>& mask)
600 {
601  is.read((char*)&val, sizeof(Real));
602  return 0;
603 }
604 
605 inline Int combine(Real& val, Real& ext)
606 {
607  val += ext;
608  return 0;
609 }
610 
611 //-------------------
612 //Complex
613 inline Int serialize(const Complex& val, std::ostream& os, const std::vector<Int>& mask)
614 {
615  os.write((char*)&val, sizeof(Complex));
616  return 0;
617 }
618 
619 inline Int deserialize(Complex& val, std::istream& is, const std::vector<Int>& mask)
620 {
621  is.read((char*)&val, sizeof(Complex));
622  return 0;
623 }
624 
625 inline Int combine(Complex& val, Complex& ext)
626 {
627  val += ext;
628  return 0;
629 }
630 
631 //-------------------
632 //Index2 TODO
633 //inline Int serialize(const Index2& val, std::ostream& os, const std::vector<Int>& mask)
634 //{
635 // os.write((char*)&(val[0]), 2*sizeof(Int));
636 // return 0;
637 //}
638 //
639 //inline Int deserialize(Index2& val, std::istream& is, const std::vector<Int>& mask)
640 //{
641 // is.read((char*)&(val[0]), 2*sizeof(Int));
642 // return 0;
643 //}
644 //
645 //inline Int combine(Index2& val, Index2& ext)
646 //{
647 // ErrorHandling( "Combine operation not implemented." );
648 // return 0;
649 //}
650 
651 //-------------------
652 //Point2 TODO
653 //inline Int serialize(const Point2& val, std::ostream& os, const std::vector<Int>& mask)
654 //{
655 // os.write((char*)&(val[0]), 2*sizeof(Real));
656 // return 0;
657 //}
658 //
659 //inline Int deserialize(Point2& val, std::istream& is, const std::vector<Int>& mask)
660 //{
661 // is.read((char*)&(val[0]), 2*sizeof(Real));
662 // return 0;
663 //}
664 //
665 //inline Int combine(Point2& val, Point2& ext)
666 //{
667 // ErrorHandling( "Combine operation not implemented." );
668 // return 0;
669 //}
670 
671 //-------------------
672 //Index3
673 inline Int serialize(const Index3& val, std::ostream& os, const std::vector<Int>& mask)
674 {
675  os.write((char*)&(val[0]), 3*sizeof(Int));
676  return 0;
677 }
678 
679 inline Int deserialize(Index3& val, std::istream& is, const std::vector<Int>& mask)
680 {
681  is.read((char*)&(val[0]), 3*sizeof(Int));
682  return 0;
683 }
684 
685 inline Int combine(Index3& val, Index3& ext)
686 {
687  ErrorHandling( "Combine operation not implemented." );
688 }
689 
690 //-------------------
691 //Point3
692 inline Int serialize(const Point3& val, std::ostream& os, const std::vector<Int>& mask)
693 {
694  os.write((char*)&(val[0]), 3*sizeof(Real));
695  return 0;
696 }
697 
698 inline Int deserialize(Point3& val, std::istream& is, const std::vector<Int>& mask)
699 {
700  is.read((char*)&(val[0]), 3*sizeof(Real));
701  return 0;
702 }
703 
704 inline Int combine(Point3& val, Point3& ext)
705 {
706  ErrorHandling( "Combine operation not implemented." );
707 }
708 
709 //-------------------
710 //std::vector
711 template<class T>
712 Int serialize(const std::vector<T>& val, std::ostream& os, const std::vector<Int>& mask)
713 {
714  Int sz = val.size();
715  os.write((char*)&sz, sizeof(Int));
716  for(Int k=0; k<sz; k++)
717  serialize(val[k], os, mask);
718  return 0;
719 }
720 
721 template<class T>
722 Int deserialize(std::vector<T>& val, std::istream& is, const std::vector<Int>& mask)
723 {
724  Int sz;
725  is.read((char*)&sz, sizeof(Int));
726  val.resize(sz);
727  for(Int k=0; k<sz; k++)
728  deserialize(val[k], is, mask);
729  return 0;
730 }
731 
732 template<class T>
733 Int combine(std::vector<T>& val, std::vector<T>& ext)
734 {
735  ErrorHandling( "Combine operation not implemented." );
736  return 0;
737 }
738 
739 //-------------------
740 //std::set
741 template<class T>
742 Int serialize(const std::set<T>& val, std::ostream& os, const std::vector<Int>& mask)
743 {
744  Int sz = val.size();
745  os.write((char*)&sz, sizeof(Int));
746  for(typename std::set<T>::const_iterator mi=val.begin(); mi!=val.end(); mi++)
747  serialize((*mi), os, mask);
748  return 0;
749 }
750 
751 template<class T>
752 Int deserialize(std::set<T>& val, std::istream& is, const std::vector<Int>& mask)
753 {
754  val.clear();
755  Int sz;
756  is.read((char*)&sz, sizeof(Int));
757  for(Int k=0; k<sz; k++) {
758  T t; deserialize(t, is, mask);
759  val.insert(t);
760  }
761  return 0;
762 }
763 
764 template<class T>
765 Int combine(std::set<T>& val, std::set<T>& ext)
766 {
767  ErrorHandling( "Combine operation not implemented." );
768  return 0;
769 }
770 
771 //-------------------
772 //std::map
773 template<class T, class S>
774 Int serialize(const std::map<T,S>& val, std::ostream& os, const std::vector<Int>& mask)
775 {
776  Int sz = val.size();
777  os.write((char*)&sz, sizeof(Int));
778  for(typename std::map<T,S>::const_iterator mi=val.begin(); mi!=val.end(); mi++) {
779  serialize((*mi).first, os, mask);
780  serialize((*mi).second, os, mask);
781  }
782  return 0;
783 }
784 
785 template<class T, class S>
786 Int deserialize(std::map<T,S>& val, std::istream& is, const std::vector<Int>& mask)
787 {
788  val.clear();
789  Int sz;
790  is.read((char*)&sz, sizeof(Int));
791  for(Int k=0; k<sz; k++) {
792  T t; deserialize(t, is, mask);
793  S s; deserialize(s, is, mask);
794  val[t] = s;
795  }
796  return 0;
797 }
798 
799 template<class T, class S>
800 Int combine(std::map<T,S>& val, std::map<T,S>& ext)
801 {
802  ErrorHandling( "Combine operation not implemented." );
803  return 0;
804 }
805 
806 //-------------------
807 //std::pair
808 template<class T, class S>
809 Int serialize(const std::pair<T,S>& val, std::ostream& os, const std::vector<Int>& mask)
810 {
811  serialize(val.first, os, mask);
812  serialize(val.second, os, mask);
813  return 0;
814 }
815 
816 template<class T, class S>
817 Int deserialize(std::pair<T,S>& val, std::istream& is, const std::vector<Int>& mask)
818 {
819  deserialize(val.first, is, mask);
820  deserialize(val.second, is, mask);
821  return 0;
822 }
823 
824 template<class T, class S>
825 Int combine(std::pair<T,S>& val, std::pair<T,S>& ext)
826 {
827  ErrorHandling( "Combine operation not implemented." );
828  return 0;
829 }
830 
831 /*
832 //-------------------
833 //BolNumVec
834 inline Int serialize(const BolNumVec& val, std::ostream& os, const std::vector<Int>& mask)
835 {
836 Int m = val.m();
837 os.write((char*)&m, sizeof(Int));
838 os.write((char*)(val.Data()), m*sizeof(bool));
839 return 0;
840 }
841 
842 inline Int deserialize(BolNumVec& val, std::istream& is, const std::vector<Int>& mask)
843 {
844 Int m;
845 is.read((char*)&m, sizeof(Int));
846 val.Resize(m);
847 is.read((char*)(val.Data()), m*sizeof(bool));
848 return 0;
849 }
850 
851 //-------------------
852 //BolNumMat
853 inline Int serialize(const BolNumMat& val, std::ostream& os, const std::vector<Int>& mask)
854 {
855 Int m = val.m();
856 Int n = val.n();
857 os.write((char*)&m, sizeof(Int));
858 os.write((char*)&n, sizeof(Int));
859 os.write((char*)(val.Data()), m*n*sizeof(bool));
860 return 0;
861 }
862 
863 inline Int deserialize(BolNumMat& val, std::istream& is, const std::vector<Int>& mask)
864 {
865 Int m;
866 Int n;
867 is.read((char*)&m, sizeof(Int));
868 is.read((char*)&n, sizeof(Int));
869 val.Resize(m,n);
870 is.read((char*)(val.Data()), m*n*sizeof(bool));
871 return 0;
872 }
873 
874 //-------------------
875 //BolNumTns
876 inline Int serialize(const BolNumTns& val, std::ostream& os, const std::vector<Int>& mask)
877 {
878 Int m = val.m(); Int n = val.n(); Int p = val.p();
879 os.write((char*)&m, sizeof(Int));
880 os.write((char*)&n, sizeof(Int));
881 os.write((char*)&p, sizeof(Int));
882 os.write((char*)(val.Data()), m*n*p*sizeof(bool));
883 return 0;
884 }
885 
886 inline Int deserialize(BolNumTns& val, std::istream& is, const std::vector<Int>& mask)
887 {
888 Int m,n,p;
889 is.read((char*)&m, sizeof(Int));
890 is.read((char*)&n, sizeof(Int));
891 is.read((char*)&p, sizeof(Int));
892 val.Resize(m,n,p);
893 is.read((char*)(val.Data()), m*n*p*sizeof(bool));
894 return 0;
895 }
896  */
897 
898 //-------------------
899 //IntNumVec
900 inline Int serialize(const IntNumVec& val, std::ostream& os, const std::vector<Int>& mask)
901 {
902  Int m = val.m();
903  os.write((char*)&m, sizeof(Int));
904  os.write((char*)(val.Data()), m*sizeof(Int));
905  return 0;
906 }
907 
908 inline Int deserialize(IntNumVec& val, std::istream& is, const std::vector<Int>& mask)
909 {
910  Int m;
911  is.read((char*)&m, sizeof(Int));
912  val.Resize(m);
913  is.read((char*)(val.Data()), m*sizeof(Int));
914  return 0;
915 }
916 
917 inline Int combine(IntNumVec& val, IntNumVec& ext)
918 {
919  //val.resize(ext.m());
920  assert(val.m()==ext.m());
921  for(Int i=0; i<val.m(); i++) val(i) += ext(i);
922  return 0;
923 }
924 
925 
926 //-------------------
927 //IntNumMat
928 inline Int serialize(const IntNumMat& val, std::ostream& os, const std::vector<Int>& mask)
929 {
930  Int m = val.m();
931  Int n = val.n();
932  os.write((char*)&m, sizeof(Int));
933  os.write((char*)&n, sizeof(Int));
934  os.write((char*)(val.Data()), m*n*sizeof(Int));
935  return 0;
936 }
937 
938 inline Int deserialize(IntNumMat& val, std::istream& is, const std::vector<Int>& mask)
939 {
940  Int m;
941  Int n;
942  is.read((char*)&m, sizeof(Int));
943  is.read((char*)&n, sizeof(Int));
944  val.Resize(m,n);
945  is.read((char*)(val.Data()), m*n*sizeof(Int));
946  return 0;
947 }
948 
949 inline Int combine(IntNumMat& val, IntNumMat& ext)
950 {
951  //val.resize(ext.m(),ext.n());
952  assert(val.m()==ext.m() && val.n()==ext.n());
953  for(Int i=0; i<val.m(); i++)
954  for(Int j=0; j<val.n(); j++)
955  val(i,j) += ext(i,j);
956  return 0;
957 }
958 
959 //-------------------
960 //IntNumTns
961 inline Int serialize(const IntNumTns& val, std::ostream& os, const std::vector<Int>& mask)
962 {
963  Int m = val.m(); Int n = val.n(); Int p = val.p();
964  os.write((char*)&m, sizeof(Int));
965  os.write((char*)&n, sizeof(Int));
966  os.write((char*)&p, sizeof(Int));
967  os.write((char*)(val.Data()), m*n*p*sizeof(Int));
968  return 0;
969 }
970 
971 inline Int deserialize(IntNumTns& val, std::istream& is, const std::vector<Int>& mask)
972 {
973  Int m,n,p;
974  is.read((char*)&m, sizeof(Int));
975  is.read((char*)&n, sizeof(Int));
976  is.read((char*)&p, sizeof(Int));
977  val.Resize(m,n,p);
978  is.read((char*)(val.Data()), m*n*p*sizeof(Int));
979  return 0;
980 }
981 
982 inline Int combine(IntNumTns& val, IntNumTns& ext)
983 {
984  //val.resize(ext.m(),ext.n(),ext.p());
985  assert(val.m()==ext.m() && val.n()==ext.n() && val.p()==ext.p());
986  for(Int i=0; i<val.m(); i++)
987  for(Int j=0; j<val.n(); j++)
988  for(Int k=0; k<val.p(); k++)
989  val(i,j,k) += ext(i,j,k);
990  return 0;
991 }
992 
993 //-------------------
994 //DblNumVec
995 inline Int serialize(const DblNumVec& val, std::ostream& os, const std::vector<Int>& mask)
996 {
997  Int m = val.m();
998  os.write((char*)&m, sizeof(Int));
999  os.write((char*)(val.Data()), m*sizeof(Real));
1000  return 0;
1001 }
1002 
1003 inline Int deserialize(DblNumVec& val, std::istream& is, const std::vector<Int>& mask)
1004 {
1005  Int m;
1006  is.read((char*)&m, sizeof(Int));
1007  val.Resize(m);
1008  is.read((char*)(val.Data()), m*sizeof(Real));
1009  return 0;
1010 }
1011 
1012 inline Int combine(DblNumVec& val, DblNumVec& ext)
1013 {
1014  //val.resize(ext.m());
1015  assert(val.m()==ext.m());
1016  for(Int i=0; i<val.m(); i++) val(i) += ext(i);
1017  return 0;
1018 }
1019 
1020 //-------------------
1021 //DblNumMat
1022 inline Int serialize(const DblNumMat& val, std::ostream& os, const std::vector<Int>& mask)
1023 {
1024  Int m = val.m();
1025  Int n = val.n();
1026  os.write((char*)&m, sizeof(Int));
1027  os.write((char*)&n, sizeof(Int));
1028  os.write((char*)(val.Data()), m*n*sizeof(Real));
1029  return 0;
1030 }
1031 
1032 inline Int deserialize(DblNumMat& val, std::istream& is, const std::vector<Int>& mask)
1033 {
1034  Int m;
1035  Int n;
1036  is.read((char*)&m, sizeof(Int));
1037  is.read((char*)&n, sizeof(Int));
1038  val.Resize(m,n);
1039  is.read((char*)(val.Data()), m*n*sizeof(Real));
1040  return 0;
1041 }
1042 
1043 inline Int combine(DblNumMat& val, DblNumMat& ext)
1044 {
1045  //val.resize(ext.m(),ext.n());
1046  assert(val.m()==ext.m() && val.n()==ext.n());
1047  for(Int i=0; i<val.m(); i++)
1048  for(Int j=0; j<val.n(); j++)
1049  val(i,j) += ext(i,j);
1050  return 0;
1051 }
1052 
1053 
1054 //-------------------
1055 //DblNumTns
1056 inline Int serialize(const DblNumTns& val, std::ostream& os, const std::vector<Int>& mask)
1057 {
1058  Int m = val.m(); Int n = val.n(); Int p = val.p();
1059  os.write((char*)&m, sizeof(Int));
1060  os.write((char*)&n, sizeof(Int));
1061  os.write((char*)&p, sizeof(Int));
1062  os.write((char*)(val.Data()), m*n*p*sizeof(Real));
1063  return 0;
1064 }
1065 
1066 inline Int deserialize(DblNumTns& val, std::istream& is, const std::vector<Int>& mask)
1067 {
1068  Int m,n,p;
1069  is.read((char*)&m, sizeof(Int));
1070  is.read((char*)&n, sizeof(Int));
1071  is.read((char*)&p, sizeof(Int));
1072  val.Resize(m,n,p);
1073  is.read((char*)(val.Data()), m*n*p*sizeof(Real));
1074  return 0;
1075 }
1076 
1077 inline Int combine(DblNumTns& val, DblNumTns& ext)
1078 {
1079  //val.resize(ext.m(),ext.n(),ext.p());
1080  assert(val.m()==ext.m() && val.n()==ext.n() && val.p()==ext.p());
1081  for(Int i=0; i<val.m(); i++)
1082  for(Int j=0; j<val.n(); j++)
1083  for(Int k=0; k<val.p(); k++)
1084  val(i,j,k) += ext(i,j,k);
1085  return 0;
1086 }
1087 
1088 //-------------------
1089 //CpxNumVec
1090 inline Int serialize(const CpxNumVec& val, std::ostream& os, const std::vector<Int>& mask)
1091 {
1092  Int m = val.m();
1093  os.write((char*)&m, sizeof(Int));
1094  os.write((char*)(val.Data()), m*sizeof(Complex));
1095  return 0;
1096 }
1097 
1098 inline Int deserialize(CpxNumVec& val, std::istream& is, const std::vector<Int>& mask)
1099 {
1100  Int m;
1101  is.read((char*)&m, sizeof(Int));
1102  val.Resize(m);
1103  is.read((char*)(val.Data()), m*sizeof(Complex));
1104  return 0;
1105 }
1106 
1107 inline Int combine(CpxNumVec& val, CpxNumVec& ext)
1108 {
1109  //val.resize(ext.m());
1110  assert(val.m()==ext.m());
1111  for(Int i=0; i<val.m(); i++) val(i) += ext(i);
1112  return 0;
1113 }
1114 
1115 //-------------------
1116 //CpxNumMat
1117 inline Int serialize(const CpxNumMat& val, std::ostream& os, const std::vector<Int>& mask)
1118 {
1119  Int m = val.m();
1120  Int n = val.n();
1121  os.write((char*)&m, sizeof(Int));
1122  os.write((char*)&n, sizeof(Int));
1123  os.write((char*)(val.Data()), m*n*sizeof(Complex));
1124  return 0;
1125 }
1126 
1127 inline Int deserialize(CpxNumMat& val, std::istream& is, const std::vector<Int>& mask)
1128 {
1129  Int m;
1130  Int n;
1131  is.read((char*)&m, sizeof(Int));
1132  is.read((char*)&n, sizeof(Int));
1133  val.Resize(m,n);
1134  is.read((char*)(val.Data()), m*n*sizeof(Complex));
1135  return 0;
1136 }
1137 
1138 inline Int combine(CpxNumMat& val, CpxNumMat& ext)
1139 {
1140  //val.resize(ext.m(),ext.n());
1141  assert(val.m()==ext.m() && val.n()==ext.n());
1142  for(Int i=0; i<val.m(); i++)
1143  for(Int j=0; j<val.n(); j++)
1144  val(i,j) += ext(i,j);
1145  return 0;
1146 }
1147 
1148 //-------------------
1149 //CpxNumTns
1150 inline Int serialize(const CpxNumTns& val, std::ostream& os, const std::vector<Int>& mask)
1151 {
1152  Int m = val.m(); Int n = val.n(); Int p = val.p();
1153  os.write((char*)&m, sizeof(Int));
1154  os.write((char*)&n, sizeof(Int));
1155  os.write((char*)&p, sizeof(Int));
1156  os.write((char*)(val.Data()), m*n*p*sizeof(Complex));
1157  return 0;
1158 }
1159 
1160 inline Int deserialize(CpxNumTns& val, std::istream& is, const std::vector<Int>& mask)
1161 {
1162  Int m,n,p;
1163  is.read((char*)&m, sizeof(Int));
1164  is.read((char*)&n, sizeof(Int));
1165  is.read((char*)&p, sizeof(Int));
1166  val.Resize(m,n,p);
1167  is.read((char*)(val.Data()), m*n*p*sizeof(Complex));
1168  return 0;
1169 }
1170 
1171 inline Int combine(CpxNumTns& val, CpxNumTns& ext)
1172 {
1173  //val.resize(ext.m(),ext.n(),ext.p());
1174  assert(val.m()==ext.m() && val.n()==ext.n() && val.p()==ext.p());
1175  for(Int i=0; i<val.m(); i++)
1176  for(Int j=0; j<val.n(); j++)
1177  for(Int k=0; k<val.p(); k++)
1178  val(i,j,k) += ext(i,j,k);
1179  return 0;
1180 }
1181 
1182 //-------------------
1183 //NumVec
1184 template<class T>
1185 Int inline serialize(const NumVec<T>& val, std::ostream& os, const std::vector<Int>& mask)
1186 {
1187  Int m = val.m();
1188  os.write((char*)&m, sizeof(Int));
1189  for(Int i=0; i<m; i++)
1190  serialize(val(i), os, mask);
1191  return 0;
1192 }
1193 
1194 template<class T>
1195 Int inline deserialize(NumVec<T>& val, std::istream& is, const std::vector<Int>& mask)
1196 {
1197  Int m;
1198  is.read((char*)&m, sizeof(Int));
1199  val.Resize(m);
1200  for(Int i=0; i<m; i++)
1201  deserialize(val(i), is, mask);
1202  return 0;
1203 }
1204 
1205 template<class T>
1206 Int inline combine(NumVec<T>& val, NumVec<T>& ext)
1207 {
1208  ErrorHandling( "Combine operation not implemented." );
1209  return 0;
1210 }
1211 
1212 //-------------------
1213 //NumMat
1214 template<class T>
1215 Int inline serialize(const NumMat<T>& val, std::ostream& os, const std::vector<Int>& mask)
1216 {
1217  Int m = val.m();
1218  Int n = val.n();
1219  os.write((char*)&m, sizeof(Int));
1220  os.write((char*)&n, sizeof(Int));
1221  for(Int j=0; j<n; j++)
1222  for(Int i=0; i<m; i++)
1223  serialize(val(i,j), os, mask);
1224  return 0;
1225 }
1226 template<class T>
1227 Int inline deserialize(NumMat<T>& val, std::istream& is, const std::vector<Int>& mask)
1228 {
1229  Int m;
1230  Int n;
1231  is.read((char*)&m, sizeof(Int));
1232  is.read((char*)&n, sizeof(Int));
1233  val.Resize(m,n);
1234  for(Int j=0; j<n; j++)
1235  for(Int i=0; i<m; i++)
1236  deserialize(val(i,j), is, mask);
1237  return 0;
1238 }
1239 
1240 template<class T>
1241 Int inline combine(NumMat<T>& val, NumMat<T>& ext)
1242 {
1243  ErrorHandling( "Combine operation not implemented." );
1244  return 0;
1245 }
1246 
1247 
1248 //-------------------
1249 //NumTns
1250 template<class T>
1251 Int inline serialize(const NumTns<T>& val, std::ostream& os, const std::vector<Int>& mask)
1252 {
1253  Int m = val.m();
1254  Int n = val.n();
1255  Int p = val.p();
1256  os.write((char*)&m, sizeof(Int));
1257  os.write((char*)&n, sizeof(Int));
1258  os.write((char*)&p, sizeof(Int));
1259  for(Int k=0; k<p; k++)
1260  for(Int j=0; j<n; j++)
1261  for(Int i=0; i<m; i++)
1262  serialize(val(i,j,k), os, mask);
1263  return 0;
1264 }
1265 
1266 template<class T>
1267 Int inline deserialize(NumTns<T>& val, std::istream& is, const std::vector<Int>& mask)
1268 {
1269  Int m;
1270  Int n;
1271  Int p;
1272  is.read((char*)&m, sizeof(Int));
1273  is.read((char*)&n, sizeof(Int));
1274  is.read((char*)&p, sizeof(Int));
1275  val.Resize(m,n,p);
1276  for(Int k=0; k<p; k++)
1277  for(Int j=0; j<n; j++)
1278  for(Int i=0; i<m; i++)
1279  deserialize(val(i,j,k), is, mask);
1280  return 0;
1281 }
1282 
1283 template<class T>
1284 Int inline combine(NumTns<T>& val, NumTns<T>& ext)
1285 {
1286  ErrorHandling( "Combine operation not implemented." );
1287  return 0;
1288 }
1289 
1290 //-------------------
1291 //DistSparseMatrix
1292 template<class T>
1293 Int inline serialize(const DistSparseMatrix<T>& val, std::ostream& os, const std::vector<Int>& mask)
1294 {
1295  serialize( val.size, os, mask );
1296  serialize( val.nnz, os, mask );
1297  serialize( val.nnzLocal, os, mask );
1298  serialize( val.colptrLocal, os, mask );
1299  serialize( val.rowindLocal, os, mask );
1300  serialize( val.nzvalLocal, os, mask );
1301  // No need to serialize the communicator
1302  return 0;
1303 }
1304 
1305 template<class T>
1306 Int inline deserialize(DistSparseMatrix<T>& val, std::istream& is, const std::vector<Int>& mask)
1307 {
1308  deserialize( val.size, is, mask );
1309  deserialize( val.nnz, is, mask );
1310  deserialize( val.nnzLocal, is, mask );
1311  deserialize( val.colptrLocal, is, mask );
1312  deserialize( val.rowindLocal, is, mask );
1313  deserialize( val.nzvalLocal, is, mask );
1314  // No need to deserialize the communicator
1315  return 0;
1316 }
1317 
1318 
1319 #ifdef WITH_SYMPACK
1320 Int inline serialize(const symPACK::DistSparseMatrixGraph & val, std::ostream& os, const std::vector<Int>& mask)
1321 {
1322  serialize( val.size, os, mask );
1323  serialize( val.nnz, os, mask );
1324  serialize( val.bIsExpanded, os, mask );
1325  serialize( val.mpirank, os, mask );
1326  serialize( val.mpisize, os, mask );
1327  serialize( val.baseval, os, mask );
1328  serialize( val.keepDiag, os, mask );
1329  serialize( val.sorted, os, mask );
1330  serialize( val.vertexDist, os, mask );
1331  serialize( val.colptr, os, mask );
1332  serialize( val.rowind, os, mask );
1333  // No need to serialize the communicator
1334  return 0;
1335 }
1336 
1337 Int inline deserialize(symPACK::DistSparseMatrixGraph& val, std::istream& is, const std::vector<Int>& mask)
1338 {
1339  deserialize( val.size, is, mask );
1340  deserialize( val.nnz, is, mask );
1341  deserialize( val.bIsExpanded, is, mask );
1342  deserialize( val.mpirank, is, mask );
1343  deserialize( val.mpisize, is, mask );
1344  deserialize( val.baseval, is, mask );
1345  deserialize( val.keepDiag, is, mask );
1346  deserialize( val.sorted, is, mask );
1347  deserialize( val.vertexDist, is, mask );
1348  deserialize( val.colptr, is, mask );
1349  deserialize( val.rowind, is, mask );
1350  // No need to deserialize the communicator
1351  return 0;
1352 }
1353 
1354 template<class T>
1355 Int inline serialize(const symPACK::DistSparseMatrix<T>& val, std::ostream& os, const std::vector<Int>& mask)
1356 {
1357  serialize( val.size, os, mask );
1358  serialize( val.nnz, os, mask );
1359  serialize( val.GetLocalGraph(), os, mask );
1360  serialize( val.nzvalLocal, os, mask );
1361  // No need to serialize the communicator
1362  return 0;
1363 }
1364 
1365 template<class T>
1366 Int inline deserialize(symPACK::DistSparseMatrix<T>& val, std::istream& is, const std::vector<Int>& mask)
1367 {
1368  deserialize( val.size, is, mask );
1369  deserialize( val.nnz, is, mask );
1370  deserialize( val.GetLocalGraph(), is, mask );
1371  deserialize( val.nzvalLocal, is, mask );
1372  // No need to deserialize the communicator
1373  return 0;
1374 }
1375 
1376 template<class T>
1377 void inline Convert(const DistSparseMatrix<T>& A, symPACK::DistSparseMatrix<T> & B)
1378 {
1379  const MPI_Comm & comm = A.comm;
1380  int mpisize,mpirank;
1381  MPI_Comm_size(comm,&mpisize);
1382  MPI_Comm_rank(comm,&mpirank);
1383 
1384  B.comm = comm;
1385  B.size = A.size;
1386  B.nnz = A.nnz;
1387 
1388  //Initialize the graph
1389  symPACK::DistSparseMatrixGraph & BGraph = B.GetLocalGraph();
1390  //Int colPerProc = A.size / mpisize;
1391  //BGraph.vertexDist.resize(mpisize+1,colPerProc);
1392  //BGraph.vertexDist[0] = 1;
1393  //std::partial_sum(BGraph.vertexDist.begin(),BGraph.vertexDist.end(),BGraph.vertexDist.begin());
1394  //BGraph.vertexDist.back() = A.size+1;
1395  BGraph.size = B.size;
1396  BGraph.nnz = B.nnz;
1397  BGraph.SetComm(B.comm);
1398  BGraph.bIsExpanded = true;
1399  BGraph.baseval = 1;
1400  BGraph.keepDiag = 1;
1401  BGraph.sorted = 1;
1402  BGraph.colptr.resize(BGraph.LocalVertexCount()+1);
1403  symPACK::bassert(A.colptrLocal.m()== BGraph.colptr.size());
1404  std::copy(A.colptrLocal.Data(),A.colptrLocal.Data()+BGraph.LocalVertexCount()+1,BGraph.colptr.data());
1405 
1406  //Copy nzval
1407  BGraph.rowind.resize(A.nnzLocal);
1408  std::copy(A.rowindLocal.Data(),A.rowindLocal.Data()+A.nnzLocal,BGraph.rowind.data());
1409 
1410  B.nzvalLocal.resize(A.nzvalLocal.m());
1411  std::copy(A.nzvalLocal.Data(),A.nzvalLocal.Data()+A.nzvalLocal.m(),B.nzvalLocal.data());
1412 }
1413 
1414 
1415 template<class T>
1416 void inline Convert(symPACK::DistSparseMatrix<T>& B, DistSparseMatrix<T> & A)
1417 {
1418 
1419  MPI_Comm & comm = B.comm;
1420  int mpisize,mpirank;
1421  MPI_Comm_size(comm,&mpisize);
1422  MPI_Comm_rank(comm,&mpirank);
1423  bool wasExpanded = B.GetLocalGraph().bIsExpanded;
1424  if(!wasExpanded){
1425  B.ExpandSymmetric();
1426  B.SortGraph();
1427  B.GetLocalGraph().SetBaseval(1);
1428  }
1429 
1430  //TODO do we have to redistribute when vertexDist is not balanced ?
1431  Int colPerProc = B.size / mpisize;
1432  if(mpirank==mpisize-1){ colPerProc = B.size - mpirank*colPerProc;}
1433  symPACK::bassert(B.GetLocalGraph().LocalVertexCount() == colPerProc);
1434 
1435  A.comm = B.comm;
1436  A.size = B.size;
1437  A.nnz = B.nnz;
1438  A.colptrLocal.Resize(B.GetLocalGraph().colptr.size());
1439  std::copy(B.GetLocalGraph().colptr.begin(), B.GetLocalGraph().colptr.end(), &A.colptrLocal[0]);
1440  A.rowindLocal.Resize(B.GetLocalGraph().rowind.size());
1441  std::copy(B.GetLocalGraph().rowind.begin(), B.GetLocalGraph().rowind.end(), &A.rowindLocal[0]);
1442  A.nnzLocal = B.GetLocalGraph().rowind.size();
1443 
1444  A.nzvalLocal.Resize(B.nzvalLocal.size());
1445  std::copy(B.nzvalLocal.begin(),B.nzvalLocal.end(),A.nzvalLocal.Data());
1446 
1447  if(!wasExpanded){
1448  B.ToLowerTriangular();
1449  }
1450 }
1451 
1452 
1453 
1454 
1455 #endif
1456 
1457 
1458 
1459 template<class T>
1460 Int inline combine(DistSparseMatrix<T>& val, DistSparseMatrix<T>& ext)
1461 {
1462  ErrorHandling( "Combine operation not implemented." );
1463  return 0;
1464 }
1465 
1466 
1467 // *********************************************************************
1468 // Parallel IO functions
1469 // *********************************************************************
1470 
1471 Int SeparateRead(std::string name, std::istringstream& is);
1472 
1473 Int SeparateWrite(std::string name, std::ostringstream& os);
1474 
1475 Int SeparateWriteAscii(std::string name, std::ostringstream& os);
1476 
1477 Int SharedRead(std::string name, std::istringstream& is);
1478 
1479 Int SharedWrite(std::string name, std::ostringstream& os);
1480 
1481 // *********************************************************************
1482 // Ej
1483 // *********************************************************************
1484 inline void IdentityCol( Int col, NumVec<Real>& vec )
1485 {
1486  for(Int i=0; i<std::min(col,vec.m()); i++)
1487  vec(i) = 0.0;
1488 
1489  if(col<vec.m())
1490  vec(col) = 1.0;
1491 
1492 
1493  for(Int i=col+1; i<vec.m(); i++)
1494  vec(i) = 0.0;
1495 }
1496 
1497 inline void IdentityCol( Int col, NumVec<Complex>& vec )
1498 {
1499  for(Int i=0; i<std::min(col,vec.m()); i++)
1500  vec(i) = Complex(0.0,0.0);
1501 
1502  if(col<vec.m())
1503  vec(col) = Complex(1.0,0.0);
1504 
1505 
1506  for(Int i=col+1; i<vec.m(); i++)
1507  vec(i) = Complex(0.0,0.0);
1508 }
1509 
1510 
1511 
1512 inline void IdentityCol( IntNumVec & cols, NumMat<Real>& mat )
1513 {
1514  for(Int j=0;j<cols.m();j++){
1515  Int col= cols(j);
1516  for(Int i=0; i<std::min(col,mat.m()); i++)
1517  mat(i,j) = 0.0;
1518 
1519  if(col<mat.m())
1520  mat(col,j) = 1.0;
1521 
1522 
1523  for(Int i=col+1; i<mat.m(); i++)
1524  mat(i,j) = 0.0;
1525  }
1526 }
1527 
1528 
1529 inline void IdentityCol( IntNumVec & cols, NumMat<Complex>& mat )
1530 {
1531  for(Int j=0;j<cols.m();j++){
1532  Int col= cols(j);
1533  for(Int i=0; i<std::min(col,mat.m()); i++)
1534  mat(i,j) = Complex(0.0,0.0);
1535 
1536  if(col<mat.m())
1537  mat(col,j) = Complex(1.0,0.0);
1538 
1539 
1540  for(Int i=col+1; i<mat.m(); i++)
1541  mat(i,j) = Complex(0.0,0.0);
1542  }
1543 }
1544 
1545 
1546 
1547 
1548 // *********************************************************************
1549 // Random numbers
1550 // *********************************************************************
1551 inline void SetRandomSeed(long int seed){
1552  srand48(seed);
1553 }
1554 
1555 inline Real UniformRandom(){
1556  return (Real)drand48();
1557 }
1558 
1559 inline void UniformRandom( NumVec<Real>& vec )
1560 {
1561  for(Int i=0; i<vec.m(); i++)
1562  vec(i) = UniformRandom();
1563 }
1564 
1565 inline void UniformRandom( NumVec<Complex>& vec )
1566 {
1567  for(Int i=0; i<vec.m(); i++)
1568  vec(i) = Complex(UniformRandom(), UniformRandom());
1569 }
1570 
1571 inline void UniformRandom( NumMat<Real>& M )
1572 {
1573  Real *ptr = M.Data();
1574  for(Int i=0; i < M.m() * M.n(); i++)
1575  *(ptr++) = UniformRandom();
1576 }
1577 
1578 inline void UniformRandom( NumMat<Complex>& M )
1579 {
1580  Complex *ptr = M.Data();
1581  for(Int i=0; i < M.m() * M.n(); i++)
1582  *(ptr++) = Complex(UniformRandom(), UniformRandom());
1583 }
1584 
1585 
1586 inline void UniformRandom( NumTns<Real>& T )
1587 {
1588  Real *ptr = T.Data();
1589  for(Int i=0; i < T.m() * T.n() * T.p(); i++)
1590  *(ptr++) = UniformRandom();
1591 }
1592 
1593 inline void UniformRandom( NumTns<Complex>& T )
1594 {
1595  Complex *ptr = T.Data();
1596  for(Int i=0; i < T.m() * T.n() * T.p(); i++)
1597  *(ptr++) = Complex(UniformRandom(), UniformRandom());
1598 }
1599 
1600 // *********************************************************************
1601 // Timing
1602 // *********************************************************************
1603 inline void GetTime(Real& t){
1604  t = MPI_Wtime();
1605 }
1606 
1607 // *********************************************************************
1608 // Comparator
1609 // *********************************************************************
1610 
1611 // Real
1612 inline bool PairLtComparator( const std::pair<Real, Int>& l,
1613  const std::pair<Real, Int>& r ){
1614  return l.first < r.first;
1615 }
1616 
1617 inline bool PairGtComparator( const std::pair<Real, Int>& l,
1618  const std::pair<Real, Int>& r ){
1619  return l.first > r.first;
1620 }
1621 
1622 // For sorting with indices
1623 // Example usage:
1624 // std::sort(val.begin(), val.end(), IndexComp<std::vector<int>&>(indices));
1625 template<class T>
1626 struct IndexComp {
1627 private:
1628  const T indices_;
1629 public:
1630  IndexComp (const T indices) : indices_(indices) {}
1631  bool operator()(const size_t a, const size_t b) const
1632  { return indices_[a] < indices_[b]; }
1633 };
1634 
1635 
1636 
1637 
1638 // *********************************************************************
1639 // Sparse Matrix
1640 // *********************************************************************
1641 
1642 template <typename T> void CSCToCSR(DistSparseMatrix<T>& sparseA, DistSparseMatrix<T> & sparseB );
1643 
1644 // Real format
1645 void ReadSparseMatrix ( const char* filename, SparseMatrix<Real>& spmat );
1646 
1647 void ReadDistSparseMatrix( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1648 
1649 void ParaReadDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1650 
1651 void ParaWriteDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1652 
1653 void ReadDistSparseMatrixFormatted( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1654 
1655 // Complex format
1656 void ParaReadDistSparseMatrix ( const char* filename, DistSparseMatrix<Complex>& pspmat, MPI_Comm comm );
1657 
1658 void ParaWriteDistSparseMatrix ( const char* filename, DistSparseMatrix<Complex>& pspmat, MPI_Comm comm );
1659 
1660 void ReadDistSparseMatrixFormatted( const char* filename, DistSparseMatrix<Complex>& pspmat, MPI_Comm comm );
1661 
1662 
1663 template <class F1, class F2>
1664 void
1665 CopyPattern ( const SparseMatrix<F1>& A, SparseMatrix<F2>& B )
1666 {
1667  B.size = A.size;
1668  B.nnz = A.nnz;
1669  B.colptr = A.colptr;
1670  B.rowind = A.rowind;
1671  B.nzval.Resize( A.nnz );
1672  return ;
1673 } // ----- end of template function CopyPattern -----
1674 
1675 
1676 // TODO Real format
1677 template<typename T>
1678 void
1679 GetDiagonal ( const DistSparseMatrix<T>& A,
1680  NumVec<T>& diag );
1681 
1682 
1683 // Functions for DistSparseMatrix
1684 
1685 template <class F1, class F2>
1686 void
1687 CopyPattern ( const DistSparseMatrix<F1>& A, DistSparseMatrix<F2>& B )
1688 {
1689  B.size = A.size;
1690  B.nnz = A.nnz;
1691  B.nnzLocal = A.nnzLocal;
1692  B.colptrLocal = A.colptrLocal;
1693  B.rowindLocal = A.rowindLocal;
1694  B.nzvalLocal.Resize( A.nnzLocal );
1695  B.comm = A.comm;
1696  return ;
1697 } // ----- end of template function CopyPattern -----
1698 
1699 #ifdef WITH_SYMPACK
1700 template <class F1, class F2>
1701 void
1702 CopyPattern ( const symPACK::DistSparseMatrix<F1>& A, symPACK::DistSparseMatrix<F2>& B )
1703 {
1704  B.size = A.size;
1705  B.nnz = A.nnz;
1706  const symPACK::DistSparseMatrixGraph & Agraph = A.GetLocalGraph();
1707  symPACK::DistSparseMatrixGraph & Bgraph = B.GetLocalGraph();
1708  //copy the graph
1709  Bgraph = Agraph;
1710  B.nzvalLocal.resize( Bgraph.LocalEdgeCount() );
1711  B.comm = A.comm;
1712  return ;
1713 } // ----- end of template function CopyPattern -----
1714 #endif
1715 
1716 
1721 
1734 template <class F>
1735 void
1736 DistSparseMatMultGlobalVec(
1737  const F alpha,
1738  const DistSparseMatrix<F>& AMat,
1739  const NumVec<F>& B,
1740  const F beta,
1741  NumVec<F>& C )
1742 {
1743  if( B.m() != AMat.size ){
1744  std::ostringstream msg;
1745  msg << std::endl
1746  << "The size of the matrix A is " << AMat.size << std::endl
1747  << "The size of the vector B is " << B.m() << std::endl
1748  << "and they must agree with each other.";
1749  ErrorHandling( msg.str().c_str() );
1750  }
1751  if( C.m() != AMat.size ){
1752  std::ostringstream msg;
1753  msg << std::endl
1754  << "The size of the matrix A is " << AMat.size << std::endl
1755  << "The size of the vector C is " << C.m() << std::endl
1756  << "and they must agree with each other.";
1757  ErrorHandling( msg.str().c_str() );
1758  }
1759 
1760  MPI_Comm comm = AMat.comm;
1761 
1762  Int mpirank, mpisize;
1763  MPI_Comm_rank( comm, &mpirank );
1764  MPI_Comm_size( comm, &mpisize );
1765 
1766  // Initiate a new vector for saving DeltaC = alpha * A * B
1767  // This is somewhat inefficient implementation but this should not
1768  // matter at this moment, since the global distribution of a dense
1769  // vector is not important anyway.
1770  NumVec<F> DeltaClocal( AMat.size );
1771  NumVec<F> DeltaC( AMat.size );
1772  SetValue( DeltaClocal, 0.0 );
1773  SetValue( DeltaC, 0.0 );
1774 
1775  Int numColLocal = AMat.colptrLocal.m() - 1;
1776  Int numColLocalFirst = AMat.size / mpisize;
1777  Int firstCol = mpirank * numColLocalFirst;
1778 
1779  // NOTE: DistSparseMatrix use 1-based indicies, while access to
1780  // vectors uses 0-based indicies.
1781  for( Int j = 0; j < numColLocal; j++ ){
1782  Int jcol = firstCol + j;
1783  for( Int i = AMat.colptrLocal(j)-1;
1784  i < AMat.colptrLocal(j+1)-1; i++ ){
1785  Int irow = AMat.rowindLocal(i) - 1;
1786  DeltaClocal( irow ) += AMat.nzvalLocal(i) * B(jcol);
1787  }
1788  } // for (j)
1789 
1790  mpi::Allreduce( DeltaClocal.Data(), DeltaC.Data(),
1791  AMat.size, MPI_SUM, comm );
1792 
1793  for( Int i = 0; i < AMat.size; i++ ){
1794  C(i) = beta * C(i) + alpha * DeltaC(i);
1795  }
1796 
1797  return ;
1798 } // ----- end of template function DistSparseMatMultGlobalVec -----
1799 
1800 // *********************************************************************
1801 // Other numerical routines
1802 // *********************************************************************
1803 
1804 // Interpolation
1813 void
1814 LinearInterpolation (
1815  const std::vector<Real>& x,
1816  const std::vector<Real>& y,
1817  const std::vector<Real>& xx,
1818  std::vector<Real>& yy );
1819 
1820 // Root finding
1821 
1835 inline Real
1836 MonotoneRootFinding (
1837  const std::vector<Real>& x,
1838  const std::vector<Real>& y,
1839  Real val )
1840 {
1841  Int numX = x.size();
1842 
1843  if( val <= y[0] || val >= y[numX-1] ){
1844  std::ostringstream msg;
1845  msg
1846  << "The root finding procedure cannot find the solution for y(x)="
1847  << val << std::endl << "here [min(y) max(y)] = [" << y[0] << ", "
1848  << y[numX-1] << "]" << std::endl;
1849  ErrorHandling( msg.str().c_str() );
1850  }
1851 
1852  std::vector<Real>::const_iterator vi = std::lower_bound(
1853  y.begin(), y.end(), val );
1854 
1855  Int idx = vi - y.begin();
1856 
1857  Real root;
1858  if( y[idx] == y[idx-1] ){
1859  root = x[idx-1];
1860  }
1861  else{
1862  // Linear interpolation
1863  root = x[idx-1] + ( x[idx] - x[idx-1] ) / ( y[idx] - y[idx-1] )
1864  * ( val - y[idx-1] );
1865  }
1866 
1867 
1868  return root;
1869 } // ----- end of function MonotoneRootFinding -----
1870 
1871 
1872 
1873 
1874 
1875 } // namespace PEXSI
1876 
1877 #include "pexsi/utility_impl.hpp"
1878 
1879 #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:153
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:1736
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:1836
PEXSI::IndexComp
Definition: utility.hpp:1626
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:180
PEXSI::DistSparseMatrix
DistSparseMatrix describes a Sparse matrix in the compressed sparse column format (CSC) and distribut...
Definition: sparse_matrix.hpp:91
Generated on Tue Nov 8 2016 21:04:58 for PEXSI by   doxygen 1.8.6