PEXSI
 All Classes Namespaces Files Functions Variables Typedefs Pages
utility_impl.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_IMPL_HPP_
47 #define _PEXSI_UTILITY_IMPL_HPP_
48 
49 using namespace std;
50 using std::ifstream;
51 using std::ofstream;
52 using std::vector;
53 using std::cerr;
54 
55 namespace PEXSI{
56 
57 // *********************************************************************
58 // Sparse Matrix
59 // *********************************************************************
60 
61 //---------------------------------------------------------
62 inline void ReadSparseMatrix ( const char* filename, SparseMatrix<Real>& spmat )
63 {
64 #ifndef _RELEASE_
65  PushCallStack("ReadSparseMatrix");
66 #endif
67 
68  // FIXME
69  // Binary format
70  if( 1 ){
71  std::istringstream iss;
72  SharedRead( filename, iss );
73  deserialize( spmat.size, iss, NO_MASK );
74  deserialize( spmat.nnz, iss, NO_MASK );
75  deserialize( spmat.colptr, iss, NO_MASK );
76  deserialize( spmat.rowind, iss, NO_MASK );
77  deserialize( spmat.nzval, iss, NO_MASK );
78  }
79 
80  // Ascii format
81  if( 0 ) {
82  ifstream fin(filename);
83  fin >> spmat.size >> spmat.nnz;
84 
85  spmat.colptr.Resize( spmat.size+1 );
86  spmat.rowind.Resize( spmat.nnz );
87  spmat.nzval.Resize ( spmat.nnz );
88 
89  for( Int i = 0; i < spmat.size + 1; i++ ){
90  fin >> spmat.colptr(i);
91  }
92 
93  for( Int i = 0; i < spmat.nnz; i++ ){
94  fin >> spmat.rowind(i);
95  }
96 
97  for( Int i = 0; i < spmat.nnz; i++ ){
98  fin >> spmat.nzval(i);
99  }
100 
101  fin.close();
102  }
103 #ifndef _RELEASE_
104  PopCallStack();
105 #endif
106 
107  return ;
108 } // ----- end of function ReadSparseMatrix -----
109 
110 //---------------------------------------------------------
111 inline void ReadDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
112 {
113 #ifndef _RELEASE_
114  PushCallStack("ReadDistSparseMatrix");
115 #endif
116  // Get the processor information within the current communicator
117  MPI_Barrier( comm );
118  Int mpirank; MPI_Comm_rank(comm, &mpirank);
119  Int mpisize; MPI_Comm_size(comm, &mpisize);
120  MPI_Status mpistat;
121  std::ifstream fin;
122 
123  // FIXME Maybe change to MPI_Comm_Dup.
124  pspmat.comm = comm;
125 
126  // Read basic information
127  if( mpirank == 0 ){
128  fin.open(filename);
129  if( !fin.good() ){
130  throw std::logic_error( "File cannot be opened!" );
131  }
132  fin.read((char*)&pspmat.size, sizeof(Int));
133  fin.read((char*)&pspmat.nnz, sizeof(Int));
134  }
135 
136  // FIXME Maybe need LongInt format to read the number of nonzeros later
137 
138  MPI_Bcast(&pspmat.size, 1, MPI_INT, 0, comm);
139  MPI_Bcast(&pspmat.nnz, 1, MPI_INT, 0, comm);
140 
141  // Read colptr
142  IntNumVec colptr(pspmat.size+1);
143  if( mpirank == 0 ){
144  Int tmp;
145  fin.read((char*)&tmp, sizeof(Int));
146 
147  if( tmp != pspmat.size+1 ){
148  throw std::logic_error( "colptr is not of the right size." );
149  }
150 
151  fin.read((char*)colptr.Data(), sizeof(Int)*tmp);
152 
153  }
154 
155  MPI_Bcast(colptr.Data(), pspmat.size+1, MPI_INT, 0, comm);
156 
157  // Compute the number of columns on each processor
158  IntNumVec numColLocalVec(mpisize);
159  Int numColLocal, numColFirst;
160  numColFirst = pspmat.size / mpisize;
161  SetValue( numColLocalVec, numColFirst );
162  numColLocalVec[mpisize-1] = pspmat.size - numColFirst * (mpisize-1); // Modify the last entry
163  numColLocal = numColLocalVec[mpirank];
164 
165  pspmat.colptrLocal.Resize( numColLocal + 1 );
166  for( Int i = 0; i < numColLocal + 1; i++ ){
167  pspmat.colptrLocal[i] = colptr[mpirank * numColFirst+i] - colptr[mpirank * numColFirst] + 1;
168  }
169 
170  // Calculate nnz_loc on each processor
171  pspmat.nnzLocal = pspmat.colptrLocal[numColLocal] - pspmat.colptrLocal[0];
172 
173  pspmat.rowindLocal.Resize( pspmat.nnzLocal );
174  pspmat.nzvalLocal.Resize ( pspmat.nnzLocal );
175 
176  // Read and distribute the row indices
177  if( mpirank == 0 ){
178  Int tmp;
179  fin.read((char*)&tmp, sizeof(Int));
180 
181  if( tmp != pspmat.nnz ){
182  std::ostringstream msg;
183  msg
184  << "The number of nonzeros in row indices do not match." << std::endl
185  << "nnz = " << pspmat.nnz << std::endl
186  << "size of row indices = " << tmp << std::endl;
187  throw std::logic_error( msg.str().c_str() );
188  }
189  IntNumVec buf;
190  Int numRead;
191  for( Int ip = 0; ip < mpisize; ip++ ){
192  numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
193  colptr[ip*numColFirst];
194  buf.Resize(numRead);
195  fin.read( (char*)buf.Data(), numRead*sizeof(Int) );
196 
197  if( ip > 0 ){
198  MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
199  MPI_Send(buf.Data(), numRead, MPI_INT, ip, 1, comm);
200  }
201  else{
202  pspmat.rowindLocal = buf;
203  }
204  }
205  }
206  else{
207  Int numRead;
208  MPI_Recv(&numRead, 1, MPI_INT, 0, 0, comm, &mpistat);
209  if( numRead != pspmat.nnzLocal ){
210  std::ostringstream msg;
211  msg << "The number of columns in row indices do not match." << std::endl
212  << "numRead = " << numRead << std::endl
213  << "nnzLocal = " << pspmat.nnzLocal << std::endl;
214  throw std::logic_error( msg.str().c_str() );
215  }
216 
217  pspmat.rowindLocal.Resize( numRead );
218  MPI_Recv( pspmat.rowindLocal.Data(), numRead, MPI_INT, 0, 1, comm, &mpistat );
219  }
220 
221  // Read and distribute the nonzero values
222  if( mpirank == 0 ){
223  Int tmp;
224  fin.read((char*)&tmp, sizeof(Int));
225 
226  if( tmp != pspmat.nnz ){
227  std::ostringstream msg;
228  msg
229  << "The number of nonzeros in values do not match." << std::endl
230  << "nnz = " << pspmat.nnz << std::endl
231  << "size of values = " << tmp << std::endl;
232  throw std::logic_error( msg.str().c_str() );
233  }
234  NumVec<Real> buf;
235  Int numRead;
236  for( Int ip = 0; ip < mpisize; ip++ ){
237  numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
238  colptr[ip*numColFirst];
239  buf.Resize(numRead);
240  fin.read( (char*)buf.Data(), numRead*sizeof(Real) );
241 
242  if( ip > 0 ){
243  MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
244  MPI_Send(buf.Data(), numRead, MPI_DOUBLE, ip, 1, comm);
245  }
246  else{
247  pspmat.nzvalLocal = buf;
248  }
249  }
250  }
251  else{
252  Int numRead;
253  MPI_Recv(&numRead, 1, MPI_INT, 0, 0, comm, &mpistat);
254  if( numRead != pspmat.nnzLocal ){
255  std::ostringstream msg;
256  msg << "The number of columns in values do not match." << std::endl
257  << "numRead = " << numRead << std::endl
258  << "nnzLocal = " << pspmat.nnzLocal << std::endl;
259  throw std::logic_error( msg.str().c_str() );
260  }
261 
262  pspmat.nzvalLocal.Resize( numRead );
263  MPI_Recv( pspmat.nzvalLocal.Data(), numRead, MPI_DOUBLE, 0, 1, comm, &mpistat );
264  }
265 
266  // Close the file
267  if( mpirank == 0 ){
268  fin.close();
269  }
270 
271 
272 
273  MPI_Barrier( comm );
274 
275 #ifndef _RELEASE_
276  PopCallStack();
277 #endif
278 
279  return ;
280 } // ----- end of function ReadDistSparseMatrix -----
281 
282 inline void ParaWriteDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
283 {
284 #ifndef _RELEASE_
285  PushCallStack("ParaWriteDistSparseMatrix");
286 #endif
287  // Get the processor information within the current communicator
288  MPI_Barrier( comm );
289  Int mpirank; MPI_Comm_rank(comm, &mpirank);
290  Int mpisize; MPI_Comm_size(comm, &mpisize);
291  MPI_Status mpistat;
292  Int err = 0;
293 
294 
295 
296  int filemode = MPI_MODE_WRONLY | MPI_MODE_CREATE | MPI_MODE_UNIQUE_OPEN;
297 
298  MPI_File fout;
299  MPI_Status status;
300 
301 
302 
303  err = MPI_File_open(comm,(char*) filename, filemode, MPI_INFO_NULL, &fout);
304 
305  if (err != MPI_SUCCESS) {
306  throw std::logic_error( "File cannot be opened!" );
307  }
308 
309  // FIXME Note that nnz uses the Int data type for consistency of writing / reading
310  // Write header
311  if( mpirank == 0 ){
312  err = MPI_File_write_at(fout, 0,(char*)&pspmat.size, 1, MPI_INT, &status);
313  err = MPI_File_write_at(fout, sizeof(Int),(char*)&pspmat.nnz, 1, MPI_INT, &status);
314  }
315 
316 
317  // Compute the number of columns on each processor
318  Int numColLocal = pspmat.colptrLocal.m()-1;
319  Int numColFirst = pspmat.size / mpisize;
320  IntNumVec colptrChunk(numColLocal+1);
321 
322  Int prev_nz = 0;
323  MPI_Exscan(&pspmat.nnzLocal, &prev_nz, 1, MPI_INT, MPI_SUM, comm);
324 
325  for( Int i = 0; i < numColLocal + 1; i++ ){
326  colptrChunk[i] = pspmat.colptrLocal[i] + prev_nz;
327  }
328 
329 
330  MPI_Datatype memtype, filetype;
331  MPI_Aint disps[6];
332  int blklens[6];
333  MPI_Datatype types[6] = {MPI_INT,MPI_INT, MPI_INT,MPI_INT, MPI_INT,MPI_DOUBLE};
334 
335  /* set block lengths (same for both types) */
336  blklens[0] = (mpirank==0)?1:0;
337  blklens[1] = numColLocal+1;
338  blklens[2] = (mpirank==0)?1:0;
339  blklens[3] = pspmat.nnzLocal;
340  blklens[4] = (mpirank==0)?1:0;
341  blklens[5] = pspmat.nnzLocal;
342 
343 
344 
345 
346  //Calculate offsets
347  MPI_Offset myColPtrOffset, myRowIdxOffset, myNzValOffset;
348  myColPtrOffset = 3*sizeof(int) + (mpirank*numColFirst)*sizeof(Int);
349  myRowIdxOffset = 3*sizeof(int) + (pspmat.size +1 + prev_nz)*sizeof(Int);
350  myNzValOffset = 4*sizeof(int) + (pspmat.size +1 + pspmat.nnz)*sizeof(Int)+ prev_nz*sizeof(Real);
351  disps[0] = 2*sizeof(int);
352  disps[1] = myColPtrOffset;
353  disps[2] = myRowIdxOffset;
354  disps[3] = sizeof(int)+myRowIdxOffset;
355  disps[4] = myNzValOffset;
356  disps[5] = sizeof(int)+myNzValOffset;
357 
358 
359 
360 #if ( _DEBUGlevel_ >= 1 )
361  char msg[200];
362  char * tmp = msg;
363  tmp += sprintf(tmp,"P%d ",mpirank);
364  for(int i = 0; i<6; ++i){
365  if(i==5)
366  tmp += sprintf(tmp, "%d [%d - %d] | ",i,disps[i],disps[i]+blklens[i]*sizeof(double));
367  else
368  tmp += sprintf(tmp, "%d [%d - %d] | ",i,disps[i],disps[i]+blklens[i]*sizeof(int));
369  }
370  tmp += sprintf(tmp,"\n");
371  printf("%s",msg);
372 #endif
373 
374 
375 
376 
377  MPI_Type_create_struct(6, blklens, disps, types, &filetype);
378  MPI_Type_commit(&filetype);
379 
380  /* create memory type */
381  Int np1 = pspmat.size+1;
382  MPI_Address( (void *)&np1, &disps[0]);
383  MPI_Address(colptrChunk.Data(), &disps[1]);
384  MPI_Address( (void *)&pspmat.nnz, &disps[2]);
385  MPI_Address((void *)pspmat.rowindLocal.Data(), &disps[3]);
386  MPI_Address( (void *)&pspmat.nnz, &disps[4]);
387  MPI_Address((void *)pspmat.nzvalLocal.Data(), &disps[5]);
388 
389  MPI_Type_create_struct(6, blklens, disps, types, &memtype);
390  MPI_Type_commit(&memtype);
391 
392 
393 
394  /* set file view */
395  err = MPI_File_set_view(fout, 0, MPI_BYTE, filetype, "native",MPI_INFO_NULL);
396 
397  /* everyone writes their own row offsets, columns, and
398  * data with one big noncontiguous write (in memory and
399  * file)
400  */
401  err = MPI_File_write_all(fout, MPI_BOTTOM, 1, memtype, &status);
402 
403  MPI_Type_free(&filetype);
404  MPI_Type_free(&memtype);
405 
406 
407 
408 
409 
410  MPI_Barrier( comm );
411 
412  MPI_File_close(&fout);
413 #ifndef _RELEASE_
414  PopCallStack();
415 #endif
416 
417  return ;
418 } // ----- end of function ParaWriteDistSparseMatrix -----
419 
420 inline void ParaReadDistSparseMatrix ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
421 {
422 #ifndef _RELEASE_
423  PushCallStack("ParaReadDistSparseMatrix");
424 #endif
425  // Get the processor information within the current communicator
426  MPI_Barrier( comm );
427  Int mpirank; MPI_Comm_rank(comm, &mpirank);
428  Int mpisize; MPI_Comm_size(comm, &mpisize);
429  MPI_Status mpistat;
430  MPI_Datatype type;
431  int lens[3];
432  MPI_Aint disps[3];
433  MPI_Datatype types[3];
434  Int err = 0;
435 
436  int filemode = MPI_MODE_RDONLY | MPI_MODE_UNIQUE_OPEN;
437 
438  MPI_File fin;
439  MPI_Status status;
440 
441  // FIXME Maybe change to MPI_Comm_Dup.
442  pspmat.comm = comm;
443 
444  err = MPI_File_open(comm,(char*) filename, filemode, MPI_INFO_NULL, &fin);
445 
446  if (err != MPI_SUCCESS) {
447  throw std::logic_error( "File cannot be opened!" );
448  }
449 
450  // FIXME Note that nnz uses the Int data type for consistency of writing / reading
451  // Read header
452  if( mpirank == 0 ){
453  err = MPI_File_read_at(fin, 0,(char*)&pspmat.size, 1, MPI_INT, &status);
454  err = MPI_File_read_at(fin, sizeof(Int),(char*)&pspmat.nnz, 1, MPI_INT, &status);
455  }
456 
457 
458  /* define a struct that describes all our data */
459  lens[0] = 1;
460  lens[1] = 1;
461  MPI_Address(&pspmat.size, &disps[0]);
462  MPI_Address(&pspmat.nnz, &disps[1]);
463  types[0] = MPI_INT;
464  types[1] = MPI_INT;
465  MPI_Type_struct(2, lens, disps, types, &type);
466  MPI_Type_commit(&type);
467 
468 
469  /* broadcast the header data to everyone */
470  MPI_Bcast(MPI_BOTTOM, 1, type, 0, comm);
471 
472  MPI_Type_free(&type);
473 
474  // Compute the number of columns on each processor
475  IntNumVec numColLocalVec(mpisize);
476  Int numColLocal, numColFirst;
477  numColFirst = pspmat.size / mpisize;
478  SetValue( numColLocalVec, numColFirst );
479  numColLocalVec[mpisize-1] = pspmat.size - numColFirst * (mpisize-1); // Modify the last entry
480  numColLocal = numColLocalVec[mpirank];
481  pspmat.colptrLocal.Resize( numColLocal + 1 );
482 
483 
484 
485  MPI_Offset myColPtrOffset = (2 + ((mpirank==0)?0:1) )*sizeof(int) + (mpirank*numColFirst)*sizeof(Int);
486 
487  Int np1 = 0;
488  lens[0] = (mpirank==0)?1:0;
489  lens[1] = numColLocal + 1;
490 
491  MPI_Address(&np1, &disps[0]);
492  MPI_Address(pspmat.colptrLocal.Data(), &disps[1]);
493 
494  MPI_Type_hindexed(2, lens, disps, MPI_INT, &type);
495  MPI_Type_commit(&type);
496 
497  err= MPI_File_read_at_all(fin, myColPtrOffset, MPI_BOTTOM, 1, type, &status);
498 
499  if (err != MPI_SUCCESS) {
500  throw std::logic_error( "error reading colptr" );
501  }
502  MPI_Type_free(&type);
503 
504  // Calculate nnz_loc on each processor
505  pspmat.nnzLocal = pspmat.colptrLocal[numColLocal] - pspmat.colptrLocal[0];
506 
507 
508  pspmat.rowindLocal.Resize( pspmat.nnzLocal );
509  pspmat.nzvalLocal.Resize ( pspmat.nnzLocal );
510 
511  //read rowIdx
512  MPI_Offset myRowIdxOffset = (3 + ((mpirank==0)?-1:0) )*sizeof(int) + (pspmat.size+1 + pspmat.colptrLocal[0])*sizeof(Int);
513 
514  lens[0] = (mpirank==0)?1:0;
515  lens[1] = pspmat.nnzLocal;
516 
517  MPI_Address(&np1, &disps[0]);
518  MPI_Address(pspmat.rowindLocal.Data(), &disps[1]);
519 
520  MPI_Type_hindexed(2, lens, disps, MPI_INT, &type);
521  MPI_Type_commit(&type);
522 
523  err= MPI_File_read_at_all(fin, myRowIdxOffset, MPI_BOTTOM, 1, type,&status);
524 
525  if (err != MPI_SUCCESS) {
526  throw std::logic_error( "error reading rowind" );
527  }
528  MPI_Type_free(&type);
529 
530 
531  //read nzval
532  MPI_Offset myNzValOffset = (3 + ((mpirank==0)?-1:0) )*sizeof(int) + (pspmat.size+1 + pspmat.nnz)*sizeof(Int) + pspmat.colptrLocal[0]*sizeof(double);
533 
534  lens[0] = (mpirank==0)?1:0;
535  lens[1] = pspmat.nnzLocal;
536 
537  MPI_Address(&np1, &disps[0]);
538  MPI_Address(pspmat.nzvalLocal.Data(), &disps[1]);
539 
540  types[0] = MPI_INT;
541  types[1] = MPI_DOUBLE;
542 
543  MPI_Type_create_struct(2, lens, disps, types, &type);
544  MPI_Type_commit(&type);
545 
546  err = MPI_File_read_at_all(fin, myNzValOffset, MPI_BOTTOM, 1, type,&status);
547 
548  if (err != MPI_SUCCESS) {
549  throw std::logic_error( "error reading nzval" );
550  }
551 
552  MPI_Type_free(&type);
553 
554 
555  //convert to local references
556  for( Int i = 1; i < numColLocal + 1; i++ ){
557  pspmat.colptrLocal[i] = pspmat.colptrLocal[i] - pspmat.colptrLocal[0] + 1;
558  }
559  pspmat.colptrLocal[0]=1;
560 
561  MPI_Barrier( comm );
562 
563  MPI_File_close(&fin);
564 #ifndef _RELEASE_
565  PopCallStack();
566 #endif
567 
568  return ;
569 } // ----- end of function ParaReadDistSparseMatrix -----
570 
571 inline void ReadDistSparseMatrixFormatted ( const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
572 {
573 #ifndef _RELEASE_
574  PushCallStack("ReadDistSparseMatrixFormatted");
575 #endif
576  // Get the processor information within the current communicator
577  MPI_Barrier( comm );
578  Int mpirank; MPI_Comm_rank(comm, &mpirank);
579  Int mpisize; MPI_Comm_size(comm, &mpisize);
580  MPI_Status mpistat;
581  std::ifstream fin;
582 
583  // FIXME Maybe change to MPI_Comm_Dup.
584  pspmat.comm = comm;
585 
586  // FIXME Maybe need LongInt format to read the number of nonzeros later
587 
588  // Read basic information
589  if( mpirank == 0 ){
590  fin.open(filename);
591  if( !fin.good() ){
592  throw std::logic_error( "File cannot be opened!" );
593  }
594  Int dummy;
595 
596  fin >> pspmat.size >> dummy;
597  fin >> pspmat.nnz;
598  // FIXME this is temporary and only applies to 4*4 matrix.
599 // fin >> dummy;
600  }
601 
602  MPI_Bcast(&pspmat.size, 1, MPI_INT, 0, comm);
603  MPI_Bcast(&pspmat.nnz, 1, MPI_INT, 0, comm);
604 
605  // Read colptr
606 
607  IntNumVec colptr(pspmat.size+1);
608  if( mpirank == 0 ){
609  Int* ptr = colptr.Data();
610  for( Int i = 0; i < pspmat.size+1; i++ )
611  fin >> *(ptr++);
612  }
613 
614  MPI_Bcast(colptr.Data(), pspmat.size+1, MPI_INT, 0, comm);
615 
616  // Compute the number of columns on each processor
617  IntNumVec numColLocalVec(mpisize);
618  Int numColLocal, numColFirst;
619  numColFirst = pspmat.size / mpisize;
620  SetValue( numColLocalVec, numColFirst );
621  numColLocalVec[mpisize-1] = pspmat.size - numColFirst * (mpisize-1); // Modify the last entry
622  numColLocal = numColLocalVec[mpirank];
623 
624  pspmat.colptrLocal.Resize( numColLocal + 1 );
625  for( Int i = 0; i < numColLocal + 1; i++ ){
626  pspmat.colptrLocal[i] = colptr[mpirank * numColFirst+i] - colptr[mpirank * numColFirst] + 1;
627  }
628 
629  // Calculate nnz_loc on each processor
630  pspmat.nnzLocal = pspmat.colptrLocal[numColLocal] - pspmat.colptrLocal[0];
631 
632  pspmat.rowindLocal.Resize( pspmat.nnzLocal );
633  pspmat.nzvalLocal.Resize ( pspmat.nnzLocal );
634 
635  // Read and distribute the row indices
636  if( mpirank == 0 ){
637  Int tmp;
638  IntNumVec buf;
639  Int numRead;
640  for( Int ip = 0; ip < mpisize; ip++ ){
641  numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
642  colptr[ip*numColFirst];
643  buf.Resize(numRead);
644  Int *ptr = buf.Data();
645  for( Int i = 0; i < numRead; i++ ){
646  fin >> *(ptr++);
647  }
648  if( ip > 0 ){
649  MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
650  MPI_Send(buf.Data(), numRead, MPI_INT, ip, 1, comm);
651  }
652  else{
653  pspmat.rowindLocal = buf;
654  }
655  }
656  }
657  else{
658  Int numRead;
659  MPI_Recv(&numRead, 1, MPI_INT, 0, 0, comm, &mpistat);
660  if( numRead != pspmat.nnzLocal ){
661  std::ostringstream msg;
662  msg << "The number of columns in row indices do not match." << std::endl
663  << "numRead = " << numRead << std::endl
664  << "nnzLocal = " << pspmat.nnzLocal << std::endl;
665  throw std::logic_error( msg.str().c_str() );
666  }
667 
668  pspmat.rowindLocal.Resize( numRead );
669  MPI_Recv( pspmat.rowindLocal.Data(), numRead, MPI_INT, 0, 1, comm, &mpistat );
670  }
671 
672 // std::cout << "Proc " << mpirank << " outputs rowindLocal.size() = "
673 // << pspmat.rowindLocal.m() << endl;
674 
675 
676  // Read and distribute the nonzero values
677  if( mpirank == 0 ){
678  Int tmp;
679  NumVec<Real> buf;
680  Int numRead;
681  for( Int ip = 0; ip < mpisize; ip++ ){
682  numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
683  colptr[ip*numColFirst];
684  buf.Resize(numRead);
685  Real *ptr = buf.Data();
686  for( Int i = 0; i < numRead; i++ ){
687  fin >> *(ptr++);
688  }
689  if( ip > 0 ){
690  MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
691  MPI_Send(buf.Data(), numRead, MPI_DOUBLE, ip, 1, comm);
692  }
693  else{
694  pspmat.nzvalLocal = buf;
695  }
696  }
697  }
698  else{
699  Int numRead;
700  MPI_Recv(&numRead, 1, MPI_INT, 0, 0, comm, &mpistat);
701  if( numRead != pspmat.nnzLocal ){
702  std::ostringstream msg;
703  msg << "The number of columns in values do not match." << std::endl
704  << "numRead = " << numRead << std::endl
705  << "nnzLocal = " << pspmat.nnzLocal << std::endl;
706  throw std::logic_error( msg.str().c_str() );
707  }
708 
709  pspmat.nzvalLocal.Resize( numRead );
710  MPI_Recv( pspmat.nzvalLocal.Data(), numRead, MPI_DOUBLE, 0, 1, comm, &mpistat );
711  }
712 
713  // Close the file
714  if( mpirank == 0 ){
715  fin.close();
716  }
717 
718 
719 
720  MPI_Barrier( comm );
721 
722 #ifndef _RELEASE_
723  PopCallStack();
724 #endif
725 
726  return ;
727 } // ----- end of function ReadDistSparseMatrixFormatted -----
728 
729  template<typename T> inline void GetDiagonal ( const DistSparseMatrix<T>& A, NumVec<T>& diag )
730 {
731 #ifndef _RELEASE_
732  PushCallStack("GetDiagonal");
733 #endif
734  Int mpirank, mpisize;
735  MPI_Comm_rank( A.comm, &mpirank );
736  MPI_Comm_size( A.comm, &mpisize );
737 
738  NumVec<T> diagLocal( A.size );
739  SetValue( diagLocal, ZERO<T>() );
740  diag.Resize( A.size );
741  SetValue( diag, ZERO<T>() );
742 
743  Int numColFirst = A.size / mpisize;
744  Int firstCol = mpirank * numColFirst;
745  Int numColLocal = A.colptrLocal.m() - 1;
746 
747 #if ( _DEBUGlevel_ >= 1 )
748  statusOFS << "numColFirst = " << numColFirst << std::endl;
749  statusOFS << "A.nzvalLocal.size = " << A.nzvalLocal.m() << std::endl;
750  statusOFS << "A.nnzLocal = " << A.nnzLocal << std::endl;
751 #endif
752 
753  // Note that the indices in DistSparseMatrix follows the FORTRAN convention
754  for( Int j = 0; j < numColLocal; j++ ){
755  Int jcol = j + firstCol + 1;
756  Int numRow = A.colptrLocal(j+1) - A.colptrLocal(j);
757  const Int* rowPtr = &A.rowindLocal( A.colptrLocal(j) - 1 );
758  // NOTE: The rows in DistSparseMatrix are not necessarily ordered.
759  // So lower_bound cannot be used here for fast searching. find has to be used.
760  const Int* ptr = find( rowPtr, rowPtr + numRow, jcol );
761  if( ptr == rowPtr + numRow ){
762  std::ostringstream msg;
763  msg << "Serious problem. Did not find the row corresponding to the column." << std::endl
764  << "This happens when j = " << j << ", jcol = " << jcol << ", and the row indices are " << std::endl
765  << IntNumVec( numRow, false, const_cast<Int*>(rowPtr) ) << std::endl;
766  throw std::logic_error( msg.str().c_str() );
767  }
768  Int diagIdx = ptr - A.rowindLocal.Data();
769  diagLocal( jcol - 1 ) = A.nzvalLocal( diagIdx );
770  }
771 
772  mpi::Allreduce( &diagLocal[0], &diag[0], A.size, MPI_SUM, A.comm );
773 
774 #ifndef _RELEASE_
775  PopCallStack();
776 #endif
777 
778  return ;
779 } // ----- end of function GetDiagonal -----
780 
781 
782 
783 
784 }
785 #endif //_PEXSI_UTILITY_IMPL_HPP_
void SetValue(NumMat< F > &M, F val)
SetValue sets a numerical matrix to a constant val.
Definition: NumMat_impl.hpp:137