46 #ifndef _PEXSI_UTILITY_HPP_
47 #define _PEXSI_UTILITY_HPP_
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;
74 const std::vector<Int> NO_MASK(1);
82 inline Int PrintBlock(std::ostream &os,
const std::string name){
84 os << std::endl<<
"*********************************************************************" << std::endl;
85 os << name << std::endl;
86 os <<
"*********************************************************************" << std::endl << std::endl;
91 inline Int Print(std::ostream &os,
const std::string name) {
92 os << std::setiosflags(std::ios::left) << name << std::endl;
96 inline Int Print(std::ostream &os,
const char* name) {
97 os << std::setiosflags(std::ios::left) << std::string(name) << std::endl;
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
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)
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)
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)
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
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)
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
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)
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
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)
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)
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)
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
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
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
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
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
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)
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;
361 os << std::setw(LENGTH_VAR_NAME) <<
"true" << std::endl;
363 os << std::setw(LENGTH_VAR_NAME) <<
"false" << std::endl;
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);
372 os << std::setw(LENGTH_VAR_NAME) <<
"true" << std::endl;
374 os << std::setw(LENGTH_VAR_NAME) <<
"false" << std::endl;
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]
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)
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)
429 template <
class F>
inline std::ostream& operator<<( std::ostream& os, const std::vector<F>& vec)
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++)
440 template <
class F>
inline std::ostream& operator<<( std::ostream& os, const NumVec<F>& vec)
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++)
459 template <
class F>
inline std::ostream& operator<<( std::ostream& os, const NumMat<F>& mat)
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++)
473 template <
class F>
inline std::ostream& operator<<( std::ostream& os, const NumTns<F>& tns)
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++) {
499 template <
typename T>
500 inline Int serialize(
const T& val, std::ostream& os,
const std::vector<Int>& mask)
502 os.write((
char*)&val,
sizeof(T));
506 template <
typename T>
507 inline Int deserialize(T& val, std::istream& is,
const std::vector<Int>& mask)
509 is.read((
char*)&val,
sizeof(T));
516 inline Int serialize(
const bool& val, std::ostream& os,
const std::vector<Int>& mask)
518 os.write((
char*)&val,
sizeof(
bool));
522 inline Int deserialize(
bool& val, std::istream& is,
const std::vector<Int>& mask)
524 is.read((
char*)&val,
sizeof(
bool));
529 inline Int serialize(
const char& val, std::ostream& os,
const std::vector<Int>& mask)
531 os.write((
char*)&val,
sizeof(
char));
535 inline Int deserialize(
char& val, std::istream& is,
const std::vector<Int>& mask)
537 is.read((
char*)&val,
sizeof(
char));
541 inline Int combine(
char& val,
char& ext)
543 throw std::logic_error(
"Combine operation not implemented." );
548 inline Int serialize(
const Int& val, std::ostream& os,
const std::vector<Int>& mask)
550 os.write((
char*)&val,
sizeof(Int));
554 inline Int deserialize(Int& val, std::istream& is,
const std::vector<Int>& mask)
556 is.read((
char*)&val,
sizeof(Int));
560 inline Int combine(Int& val, Int& ext)
568 inline Int serialize(
const LongInt& val, std::ostream& os,
const std::vector<Int>& mask)
570 os.write((
char*)&val,
sizeof(LongInt));
574 inline Int deserialize(LongInt& val, std::istream& is,
const std::vector<Int>& mask)
576 is.read((
char*)&val,
sizeof(LongInt));
580 inline Int combine(LongInt& val, LongInt& ext)
589 inline Int serialize(
const Real& val, std::ostream& os,
const std::vector<Int>& mask)
591 os.write((
char*)&val,
sizeof(Real));
595 inline Int deserialize(Real& val, std::istream& is,
const std::vector<Int>& mask)
597 is.read((
char*)&val,
sizeof(Real));
601 inline Int combine(Real& val, Real& ext)
609 inline Int serialize(
const Complex& val, std::ostream& os,
const std::vector<Int>& mask)
611 os.write((
char*)&val,
sizeof(Complex));
615 inline Int deserialize(Complex& val, std::istream& is,
const std::vector<Int>& mask)
617 is.read((
char*)&val,
sizeof(Complex));
621 inline Int combine(Complex& val, Complex& ext)
669 inline Int serialize(
const Index3& val, std::ostream& os,
const std::vector<Int>& mask)
671 os.write((
char*)&(val[0]), 3*
sizeof(Int));
675 inline Int deserialize(Index3& val, std::istream& is,
const std::vector<Int>& mask)
677 is.read((
char*)&(val[0]), 3*
sizeof(Int));
681 inline Int combine(Index3& val, Index3& ext)
683 throw std::logic_error(
"Combine operation not implemented." );
688 inline Int serialize(
const Point3& val, std::ostream& os,
const std::vector<Int>& mask)
690 os.write((
char*)&(val[0]), 3*
sizeof(Real));
694 inline Int deserialize(Point3& val, std::istream& is,
const std::vector<Int>& mask)
696 is.read((
char*)&(val[0]), 3*
sizeof(Real));
700 inline Int combine(Point3& val, Point3& ext)
702 throw std::logic_error(
"Combine operation not implemented." );
708 Int serialize(
const std::vector<T>& val, std::ostream& os,
const std::vector<Int>& mask)
711 os.write((
char*)&sz,
sizeof(Int));
712 for(Int k=0; k<sz; k++)
713 serialize(val[k], os, mask);
718 Int deserialize(std::vector<T>& val, std::istream& is,
const std::vector<Int>& mask)
721 is.read((
char*)&sz,
sizeof(Int));
723 for(Int k=0; k<sz; k++)
724 deserialize(val[k], is, mask);
729 Int combine(std::vector<T>& val, std::vector<T>& ext)
731 throw std::logic_error(
"Combine operation not implemented." );
738 Int serialize(
const std::set<T>& val, std::ostream& os,
const std::vector<Int>& mask)
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);
748 Int deserialize(std::set<T>& val, std::istream& is,
const std::vector<Int>& mask)
752 is.read((
char*)&sz,
sizeof(Int));
753 for(Int k=0; k<sz; k++) {
754 T t; deserialize(t, is, mask);
761 Int combine(std::set<T>& val, std::set<T>& ext)
763 throw std::logic_error(
"Combine operation not implemented." );
769 template<
class T,
class S>
770 Int serialize(
const std::map<T,S>& val, std::ostream& os,
const std::vector<Int>& mask)
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);
781 template<
class T,
class S>
782 Int deserialize(std::map<T,S>& val, std::istream& is,
const std::vector<Int>& mask)
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);
795 template<
class T,
class S>
796 Int combine(std::map<T,S>& val, std::map<T,S>& ext)
798 throw std::logic_error(
"Combine operation not implemented." );
804 template<
class T,
class S>
805 Int serialize(
const std::pair<T,S>& val, std::ostream& os,
const std::vector<Int>& mask)
807 serialize(val.first, os, mask);
808 serialize(val.second, os, mask);
812 template<
class T,
class S>
813 Int deserialize(std::pair<T,S>& val, std::istream& is,
const std::vector<Int>& mask)
815 deserialize(val.first, is, mask);
816 deserialize(val.second, is, mask);
820 template<
class T,
class S>
821 Int combine(std::pair<T,S>& val, std::pair<T,S>& ext)
823 throw std::logic_error(
"Combine operation not implemented." );
896 inline Int serialize(
const IntNumVec& val, std::ostream& os,
const std::vector<Int>& mask)
899 os.write((
char*)&m,
sizeof(Int));
900 os.write((
char*)(val.Data()), m*
sizeof(Int));
904 inline Int deserialize(IntNumVec& val, std::istream& is,
const std::vector<Int>& mask)
907 is.read((
char*)&m,
sizeof(Int));
909 is.read((
char*)(val.Data()), m*
sizeof(Int));
913 inline Int combine(IntNumVec& val, IntNumVec& ext)
916 assert(val.m()==ext.m());
917 for(Int i=0; i<val.m(); i++) val(i) += ext(i);
924 inline Int serialize(
const IntNumMat& val, std::ostream& os,
const std::vector<Int>& mask)
928 os.write((
char*)&m,
sizeof(Int));
929 os.write((
char*)&n,
sizeof(Int));
930 os.write((
char*)(val.Data()), m*n*
sizeof(Int));
934 inline Int deserialize(IntNumMat& val, std::istream& is,
const std::vector<Int>& mask)
938 is.read((
char*)&m,
sizeof(Int));
939 is.read((
char*)&n,
sizeof(Int));
941 is.read((
char*)(val.Data()), m*n*
sizeof(Int));
945 inline Int combine(IntNumMat& val, IntNumMat& ext)
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);
957 inline Int serialize(
const IntNumTns& val, std::ostream& os,
const std::vector<Int>& mask)
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));
967 inline Int deserialize(IntNumTns& val, std::istream& is,
const std::vector<Int>& mask)
970 is.read((
char*)&m,
sizeof(Int));
971 is.read((
char*)&n,
sizeof(Int));
972 is.read((
char*)&p,
sizeof(Int));
974 is.read((
char*)(val.Data()), m*n*p*
sizeof(Int));
978 inline Int combine(IntNumTns& val, IntNumTns& ext)
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);
991 inline Int serialize(
const DblNumVec& val, std::ostream& os,
const std::vector<Int>& mask)
994 os.write((
char*)&m,
sizeof(Int));
995 os.write((
char*)(val.Data()), m*
sizeof(Real));
999 inline Int deserialize(DblNumVec& val, std::istream& is,
const std::vector<Int>& mask)
1002 is.read((
char*)&m,
sizeof(Int));
1004 is.read((
char*)(val.Data()), m*
sizeof(Real));
1008 inline Int combine(DblNumVec& val, DblNumVec& ext)
1011 assert(val.m()==ext.m());
1012 for(Int i=0; i<val.m(); i++) val(i) += ext(i);
1018 inline Int serialize(
const DblNumMat& val, std::ostream& os,
const std::vector<Int>& mask)
1022 os.write((
char*)&m,
sizeof(Int));
1023 os.write((
char*)&n,
sizeof(Int));
1024 os.write((
char*)(val.Data()), m*n*
sizeof(Real));
1028 inline Int deserialize(DblNumMat& val, std::istream& is,
const std::vector<Int>& mask)
1032 is.read((
char*)&m,
sizeof(Int));
1033 is.read((
char*)&n,
sizeof(Int));
1035 is.read((
char*)(val.Data()), m*n*
sizeof(Real));
1039 inline Int combine(DblNumMat& val, DblNumMat& ext)
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);
1052 inline Int serialize(
const DblNumTns& val, std::ostream& os,
const std::vector<Int>& mask)
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));
1062 inline Int deserialize(DblNumTns& val, std::istream& is,
const std::vector<Int>& mask)
1065 is.read((
char*)&m,
sizeof(Int));
1066 is.read((
char*)&n,
sizeof(Int));
1067 is.read((
char*)&p,
sizeof(Int));
1069 is.read((
char*)(val.Data()), m*n*p*
sizeof(Real));
1073 inline Int combine(DblNumTns& val, DblNumTns& ext)
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);
1086 inline Int serialize(
const CpxNumVec& val, std::ostream& os,
const std::vector<Int>& mask)
1089 os.write((
char*)&m,
sizeof(Int));
1090 os.write((
char*)(val.Data()), m*
sizeof(Complex));
1094 inline Int deserialize(CpxNumVec& val, std::istream& is,
const std::vector<Int>& mask)
1097 is.read((
char*)&m,
sizeof(Int));
1099 is.read((
char*)(val.Data()), m*
sizeof(Complex));
1103 inline Int combine(CpxNumVec& val, CpxNumVec& ext)
1106 assert(val.m()==ext.m());
1107 for(Int i=0; i<val.m(); i++) val(i) += ext(i);
1113 inline Int serialize(
const CpxNumMat& val, std::ostream& os,
const std::vector<Int>& mask)
1117 os.write((
char*)&m,
sizeof(Int));
1118 os.write((
char*)&n,
sizeof(Int));
1119 os.write((
char*)(val.Data()), m*n*
sizeof(Complex));
1123 inline Int deserialize(CpxNumMat& val, std::istream& is,
const std::vector<Int>& mask)
1127 is.read((
char*)&m,
sizeof(Int));
1128 is.read((
char*)&n,
sizeof(Int));
1130 is.read((
char*)(val.Data()), m*n*
sizeof(Complex));
1134 inline Int combine(CpxNumMat& val, CpxNumMat& ext)
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);
1146 inline Int serialize(
const CpxNumTns& val, std::ostream& os,
const std::vector<Int>& mask)
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));
1156 inline Int deserialize(CpxNumTns& val, std::istream& is,
const std::vector<Int>& mask)
1159 is.read((
char*)&m,
sizeof(Int));
1160 is.read((
char*)&n,
sizeof(Int));
1161 is.read((
char*)&p,
sizeof(Int));
1163 is.read((
char*)(val.Data()), m*n*p*
sizeof(Complex));
1167 inline Int combine(CpxNumTns& val, CpxNumTns& ext)
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);
1181 Int
inline serialize(
const NumVec<T>& val, std::ostream& os,
const std::vector<Int>& mask)
1184 os.write((
char*)&m,
sizeof(Int));
1185 for(Int i=0; i<m; i++)
1186 serialize(val(i), os, mask);
1191 Int
inline deserialize(NumVec<T>& val, std::istream& is,
const std::vector<Int>& mask)
1194 is.read((
char*)&m,
sizeof(Int));
1196 for(Int i=0; i<m; i++)
1197 deserialize(val(i), is, mask);
1202 Int
inline combine(NumVec<T>& val, NumVec<T>& ext)
1204 throw std::logic_error(
"Combine operation not implemented." );
1211 Int
inline serialize(
const NumMat<T>& val, std::ostream& os,
const std::vector<Int>& mask)
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);
1223 Int
inline deserialize(NumMat<T>& val, std::istream& is,
const std::vector<Int>& mask)
1227 is.read((
char*)&m,
sizeof(Int));
1228 is.read((
char*)&n,
sizeof(Int));
1230 for(Int j=0; j<n; j++)
1231 for(Int i=0; i<m; i++)
1232 deserialize(val(i,j), is, mask);
1237 Int
inline combine(NumMat<T>& val, NumMat<T>& ext)
1239 throw std::logic_error(
"Combine operation not implemented." );
1247 Int
inline serialize(
const NumTns<T>& val, std::ostream& os,
const std::vector<Int>& mask)
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);
1263 Int
inline deserialize(NumTns<T>& val, std::istream& is,
const std::vector<Int>& mask)
1268 is.read((
char*)&m,
sizeof(Int));
1269 is.read((
char*)&n,
sizeof(Int));
1270 is.read((
char*)&p,
sizeof(Int));
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);
1280 Int
inline combine(NumTns<T>& val, NumTns<T>& ext)
1282 throw std::logic_error(
"Combine operation not implemented." );
1289 Int
inline serialize(
const DistSparseMatrix<T>& val, std::ostream& os,
const std::vector<Int>& mask)
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 );
1302 Int
inline deserialize(DistSparseMatrix<T>& val, std::istream& is,
const std::vector<Int>& mask)
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 );
1315 Int
inline combine(DistSparseMatrix<T>& val, DistSparseMatrix<T>& ext)
1317 throw std::logic_error(
"Combine operation not implemented." );
1326 Int SeparateRead(std::string name, std::istringstream& is);
1328 Int SeparateWrite(std::string name, std::ostringstream& os);
1330 Int SeparateWriteAscii(std::string name, std::ostringstream& os);
1332 Int SharedRead(std::string name, std::istringstream& is);
1334 Int SharedWrite(std::string name, std::ostringstream& os);
1339 inline void IdentityCol( Int col, NumVec<Real>& vec )
1341 for(Int i=0; i<std::min(col,vec.m()); i++)
1348 for(Int i=col+1; i<vec.m(); i++)
1352 inline void IdentityCol( Int col, NumVec<Complex>& vec )
1354 for(Int i=0; i<std::min(col,vec.m()); i++)
1355 vec(i) = Complex(0.0,0.0);
1358 vec(col) = Complex(1.0,0.0);
1361 for(Int i=col+1; i<vec.m(); i++)
1362 vec(i) = Complex(0.0,0.0);
1367 inline void IdentityCol( IntNumVec & cols, NumMat<Real>& mat )
1369 for(Int j=0;j<cols.m();j++){
1371 for(Int i=0; i<std::min(col,mat.m()); i++)
1378 for(Int i=col+1; i<mat.m(); i++)
1384 inline void IdentityCol( IntNumVec & cols, NumMat<Complex>& mat )
1386 for(Int j=0;j<cols.m();j++){
1388 for(Int i=0; i<std::min(col,mat.m()); i++)
1389 mat(i,j) = Complex(0.0,0.0);
1392 mat(col,j) = Complex(1.0,0.0);
1395 for(Int i=col+1; i<mat.m(); i++)
1396 mat(i,j) = Complex(0.0,0.0);
1406 inline void SetRandomSeed(
long int seed){
1410 inline Real UniformRandom(){
1411 return (Real)drand48();
1414 inline void UniformRandom( NumVec<Real>& vec )
1416 for(Int i=0; i<vec.m(); i++)
1417 vec(i) = UniformRandom();
1420 inline void UniformRandom( NumVec<Complex>& vec )
1422 for(Int i=0; i<vec.m(); i++)
1423 vec(i) = Complex(UniformRandom(), UniformRandom());
1426 inline void UniformRandom( NumMat<Real>& M )
1428 Real *ptr = M.Data();
1429 for(Int i=0; i < M.m() * M.n(); i++)
1430 *(ptr++) = UniformRandom();
1433 inline void UniformRandom( NumMat<Complex>& M )
1435 Complex *ptr = M.Data();
1436 for(Int i=0; i < M.m() * M.n(); i++)
1437 *(ptr++) = Complex(UniformRandom(), UniformRandom());
1441 inline void UniformRandom( NumTns<Real>& T )
1443 Real *ptr = T.Data();
1444 for(Int i=0; i < T.m() * T.n() * T.p(); i++)
1445 *(ptr++) = UniformRandom();
1448 inline void UniformRandom( NumTns<Complex>& T )
1450 Complex *ptr = T.Data();
1451 for(Int i=0; i < T.m() * T.n() * T.p(); i++)
1452 *(ptr++) = Complex(UniformRandom(), UniformRandom());
1458 inline void GetTime(Real& t){
1467 inline bool PairLtComparator(
const std::pair<Real, Int>& l,
1468 const std::pair<Real, Int>& r ){
1469 return l.first < r.first;
1472 inline bool PairGtComparator(
const std::pair<Real, Int>& l,
1473 const std::pair<Real, Int>& r ){
1474 return l.first > r.first;
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]; }
1508 void ReadDistSparseMatrixFormatted(
const char* filename,
DistSparseMatrix<Real>& pspmat, MPI_Comm comm );
1518 template <
class F1,
class F2>
1523 PushCallStack(
"CopyPattern");
1527 B.colptr = A.colptr;
1528 B.rowind = A.rowind;
1529 B.nzval.Resize( A.nnz );
1538 template<
typename T>
1540 GetDiagonal (
const DistSparseMatrix<T>& A,
1546 template <
class F1,
class F2>
1548 CopyPattern (
const DistSparseMatrix<F1>& A, DistSparseMatrix<F2>& B )
1551 PushCallStack(
"CopyPattern");
1555 B.nnzLocal = A.nnzLocal;
1556 B.colptrLocal = A.colptrLocal;
1557 B.rowindLocal = A.rowindLocal;
1558 B.nzvalLocal.Resize( A.nnzLocal );
1593 PushCallStack(
"DistSparseMatMultGlobalVec");
1595 if( B.m() != AMat.
size ){
1596 std::ostringstream msg;
1598 <<
"The size of the matrix A is " << AMat.
size << std::endl
1599 <<
"The size of the vector B is " << B.m() << std::endl
1600 <<
"and they must agree with each other.";
1604 throw std::runtime_error( msg.str().c_str() );
1606 if( C.m() != AMat.
size ){
1607 std::ostringstream msg;
1609 <<
"The size of the matrix A is " << AMat.
size << std::endl
1610 <<
"The size of the vector C is " << C.m() << std::endl
1611 <<
"and they must agree with each other.";
1615 throw std::runtime_error( msg.str().c_str() );
1618 MPI_Comm comm = AMat.
comm;
1620 Int mpirank, mpisize;
1621 MPI_Comm_rank( comm, &mpirank );
1622 MPI_Comm_size( comm, &mpisize );
1634 Int numColLocalFirst = AMat.
size / mpisize;
1635 Int firstCol = mpirank * numColLocalFirst;
1639 for( Int j = 0; j < numColLocal; j++ ){
1640 Int jcol = firstCol + j;
1644 DeltaClocal( irow ) += AMat.
nzvalLocal(i) * B(jcol);
1648 mpi::Allreduce( DeltaClocal.Data(), DeltaC.Data(),
1649 AMat.
size, MPI_SUM, comm );
1651 for( Int i = 0; i < AMat.
size; i++ ){
1652 C(i) = beta * C(i) + alpha * DeltaC(i);
1676 const std::vector<Real>& x,
1677 const std::vector<Real>& y,
1678 const std::vector<Real>& xx,
1679 std::vector<Real>& yy );
1698 const std::vector<Real>& x,
1699 const std::vector<Real>& y,
1703 PushCallStack(
"MonotoneRootFinding");
1705 Int numX = x.size();
1707 if( val <= y[0] || val >= y[numX-1] ){
1708 std::ostringstream msg;
1710 <<
"The root finding procedure cannot find the solution for y(x)="
1711 << val << std::endl <<
"here [min(y) max(y)] = [" << y[0] <<
", "
1712 << y[numX-1] <<
"]" << std::endl;
1716 throw std::runtime_error( msg.str().c_str() );
1719 std::vector<Real>::const_iterator vi = std::lower_bound(
1720 y.begin(), y.end(), val );
1722 Int idx = vi - y.begin();
1725 if( y[idx] == y[idx-1] ){
1730 root = x[idx-1] + ( x[idx] - x[idx-1] ) / ( y[idx] - y[idx-1] )
1731 * ( val - y[idx-1] );
1749 #endif // _PEXSI_UTILITY_HPP_
Numerical vector.
Definition: NumVec.hpp:60
IntNumVec colptrLocal
Dimension numColLocal + 1, storing the pointers to the nonzero row indices and nonzero values in rowp...
Definition: sparse_matrix.hpp:108
void SetValue(NumMat< F > &M, F val)
SetValue sets a numerical matrix to a constant val.
Definition: NumMat_impl.hpp:171
Implementation of utility subroutines.
NumVec< F > nzvalLocal
Dimension nnzLocal, storing the nonzero values.
Definition: sparse_matrix.hpp:116
SparseMatrix describes a sequential sparse matrix saved in compressed sparse column format...
Definition: sparse_matrix.hpp:66
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 and Distributed sparse matrix in compressed column format.
Interface with MPI to facilitate communication.
Int size
Matrix dimension.
Definition: sparse_matrix.hpp:93
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:1585
MPI_Comm comm
MPI communicator.
Definition: sparse_matrix.hpp:119
Tiny vectors of dimension 3.
Real MonotoneRootFinding(const std::vector< Real > &x, const std::vector< Real > &y, Real val)
Root finding for monotonic non-decreasing functions.
Definition: utility.hpp:1697
Definition: utility.hpp:1481
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:224
DistSparseMatrix describes a Sparse matrix in the compressed sparse column format (CSC) and distribut...
Definition: sparse_matrix.hpp:91