46 #ifndef _PEXSI_UTILITY_IMPL_HPP_
47 #define _PEXSI_UTILITY_IMPL_HPP_
62 inline void ReadSparseMatrix (
const char* filename, SparseMatrix<Real>& spmat )
65 PushCallStack(
"ReadSparseMatrix");
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 );
82 ifstream fin(filename);
83 fin >> spmat.size >> spmat.nnz;
85 spmat.colptr.Resize( spmat.size+1 );
86 spmat.rowind.Resize( spmat.nnz );
87 spmat.nzval.Resize ( spmat.nnz );
89 for( Int i = 0; i < spmat.size + 1; i++ ){
90 fin >> spmat.colptr(i);
93 for( Int i = 0; i < spmat.nnz; i++ ){
94 fin >> spmat.rowind(i);
97 for( Int i = 0; i < spmat.nnz; i++ ){
98 fin >> spmat.nzval(i);
111 inline void ReadDistSparseMatrix (
const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
114 PushCallStack(
"ReadDistSparseMatrix");
118 Int mpirank; MPI_Comm_rank(comm, &mpirank);
119 Int mpisize; MPI_Comm_size(comm, &mpisize);
130 throw std::logic_error(
"File cannot be opened!" );
132 fin.read((
char*)&pspmat.size,
sizeof(Int));
133 fin.read((
char*)&pspmat.nnz,
sizeof(Int));
138 MPI_Bcast(&pspmat.size, 1, MPI_INT, 0, comm);
139 MPI_Bcast(&pspmat.nnz, 1, MPI_INT, 0, comm);
142 IntNumVec colptr(pspmat.size+1);
145 fin.read((
char*)&tmp,
sizeof(Int));
147 if( tmp != pspmat.size+1 ){
148 throw std::logic_error(
"colptr is not of the right size." );
151 fin.read((
char*)colptr.Data(),
sizeof(Int)*tmp);
155 MPI_Bcast(colptr.Data(), pspmat.size+1, MPI_INT, 0, comm);
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);
163 numColLocal = numColLocalVec[mpirank];
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;
171 pspmat.nnzLocal = pspmat.colptrLocal[numColLocal] - pspmat.colptrLocal[0];
173 pspmat.rowindLocal.Resize( pspmat.nnzLocal );
174 pspmat.nzvalLocal.Resize ( pspmat.nnzLocal );
179 fin.read((
char*)&tmp,
sizeof(Int));
181 if( tmp != pspmat.nnz ){
182 std::ostringstream 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() );
191 for( Int ip = 0; ip < mpisize; ip++ ){
192 numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
193 colptr[ip*numColFirst];
195 fin.read( (
char*)buf.Data(), numRead*
sizeof(Int) );
198 MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
199 MPI_Send(buf.Data(), numRead, MPI_INT, ip, 1, comm);
202 pspmat.rowindLocal = buf;
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() );
217 pspmat.rowindLocal.Resize( numRead );
218 MPI_Recv( pspmat.rowindLocal.Data(), numRead, MPI_INT, 0, 1, comm, &mpistat );
224 fin.read((
char*)&tmp,
sizeof(Int));
226 if( tmp != pspmat.nnz ){
227 std::ostringstream 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() );
236 for( Int ip = 0; ip < mpisize; ip++ ){
237 numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
238 colptr[ip*numColFirst];
240 fin.read( (
char*)buf.Data(), numRead*
sizeof(Real) );
243 MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
244 MPI_Send(buf.Data(), numRead, MPI_DOUBLE, ip, 1, comm);
247 pspmat.nzvalLocal = buf;
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() );
262 pspmat.nzvalLocal.Resize( numRead );
263 MPI_Recv( pspmat.nzvalLocal.Data(), numRead, MPI_DOUBLE, 0, 1, comm, &mpistat );
282 inline void ParaWriteDistSparseMatrix (
const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
285 PushCallStack(
"ParaWriteDistSparseMatrix");
289 Int mpirank; MPI_Comm_rank(comm, &mpirank);
290 Int mpisize; MPI_Comm_size(comm, &mpisize);
296 int filemode = MPI_MODE_WRONLY | MPI_MODE_CREATE | MPI_MODE_UNIQUE_OPEN;
303 err = MPI_File_open(comm,(
char*) filename, filemode, MPI_INFO_NULL, &fout);
305 if (err != MPI_SUCCESS) {
306 throw std::logic_error(
"File cannot be opened!" );
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);
318 Int numColLocal = pspmat.colptrLocal.m()-1;
319 Int numColFirst = pspmat.size / mpisize;
320 IntNumVec colptrChunk(numColLocal+1);
323 MPI_Exscan(&pspmat.nnzLocal, &prev_nz, 1, MPI_INT, MPI_SUM, comm);
325 for( Int i = 0; i < numColLocal + 1; i++ ){
326 colptrChunk[i] = pspmat.colptrLocal[i] + prev_nz;
330 MPI_Datatype memtype, filetype;
333 MPI_Datatype types[6] = {MPI_INT,MPI_INT, MPI_INT,MPI_INT, MPI_INT,MPI_DOUBLE};
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;
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;
360 #if ( _DEBUGlevel_ >= 1 )
363 tmp += sprintf(tmp,
"P%d ",mpirank);
364 for(
int i = 0; i<6; ++i){
366 tmp += sprintf(tmp,
"%d [%d - %d] | ",i,disps[i],disps[i]+blklens[i]*
sizeof(
double));
368 tmp += sprintf(tmp,
"%d [%d - %d] | ",i,disps[i],disps[i]+blklens[i]*
sizeof(
int));
370 tmp += sprintf(tmp,
"\n");
377 MPI_Type_create_struct(6, blklens, disps, types, &filetype);
378 MPI_Type_commit(&filetype);
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]);
389 MPI_Type_create_struct(6, blklens, disps, types, &memtype);
390 MPI_Type_commit(&memtype);
395 err = MPI_File_set_view(fout, 0, MPI_BYTE, filetype,
"native",MPI_INFO_NULL);
401 err = MPI_File_write_all(fout, MPI_BOTTOM, 1, memtype, &status);
403 MPI_Type_free(&filetype);
404 MPI_Type_free(&memtype);
412 MPI_File_close(&fout);
420 inline void ParaReadDistSparseMatrix (
const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
423 PushCallStack(
"ParaReadDistSparseMatrix");
427 Int mpirank; MPI_Comm_rank(comm, &mpirank);
428 Int mpisize; MPI_Comm_size(comm, &mpisize);
433 MPI_Datatype types[3];
436 int filemode = MPI_MODE_RDONLY | MPI_MODE_UNIQUE_OPEN;
444 err = MPI_File_open(comm,(
char*) filename, filemode, MPI_INFO_NULL, &fin);
446 if (err != MPI_SUCCESS) {
447 throw std::logic_error(
"File cannot be opened!" );
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);
461 MPI_Address(&pspmat.size, &disps[0]);
462 MPI_Address(&pspmat.nnz, &disps[1]);
465 MPI_Type_struct(2, lens, disps, types, &type);
466 MPI_Type_commit(&type);
470 MPI_Bcast(MPI_BOTTOM, 1, type, 0, comm);
472 MPI_Type_free(&type);
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);
480 numColLocal = numColLocalVec[mpirank];
481 pspmat.colptrLocal.Resize( numColLocal + 1 );
485 MPI_Offset myColPtrOffset = (2 + ((mpirank==0)?0:1) )*
sizeof(
int) + (mpirank*numColFirst)*
sizeof(Int);
488 lens[0] = (mpirank==0)?1:0;
489 lens[1] = numColLocal + 1;
491 MPI_Address(&np1, &disps[0]);
492 MPI_Address(pspmat.colptrLocal.Data(), &disps[1]);
494 MPI_Type_hindexed(2, lens, disps, MPI_INT, &type);
495 MPI_Type_commit(&type);
497 err= MPI_File_read_at_all(fin, myColPtrOffset, MPI_BOTTOM, 1, type, &status);
499 if (err != MPI_SUCCESS) {
500 throw std::logic_error(
"error reading colptr" );
502 MPI_Type_free(&type);
505 pspmat.nnzLocal = pspmat.colptrLocal[numColLocal] - pspmat.colptrLocal[0];
508 pspmat.rowindLocal.Resize( pspmat.nnzLocal );
509 pspmat.nzvalLocal.Resize ( pspmat.nnzLocal );
512 MPI_Offset myRowIdxOffset = (3 + ((mpirank==0)?-1:0) )*
sizeof(
int) + (pspmat.size+1 + pspmat.colptrLocal[0])*
sizeof(Int);
514 lens[0] = (mpirank==0)?1:0;
515 lens[1] = pspmat.nnzLocal;
517 MPI_Address(&np1, &disps[0]);
518 MPI_Address(pspmat.rowindLocal.Data(), &disps[1]);
520 MPI_Type_hindexed(2, lens, disps, MPI_INT, &type);
521 MPI_Type_commit(&type);
523 err= MPI_File_read_at_all(fin, myRowIdxOffset, MPI_BOTTOM, 1, type,&status);
525 if (err != MPI_SUCCESS) {
526 throw std::logic_error(
"error reading rowind" );
528 MPI_Type_free(&type);
532 MPI_Offset myNzValOffset = (3 + ((mpirank==0)?-1:0) )*
sizeof(
int) + (pspmat.size+1 + pspmat.nnz)*
sizeof(Int) + pspmat.colptrLocal[0]*
sizeof(double);
534 lens[0] = (mpirank==0)?1:0;
535 lens[1] = pspmat.nnzLocal;
537 MPI_Address(&np1, &disps[0]);
538 MPI_Address(pspmat.nzvalLocal.Data(), &disps[1]);
541 types[1] = MPI_DOUBLE;
543 MPI_Type_create_struct(2, lens, disps, types, &type);
544 MPI_Type_commit(&type);
546 err = MPI_File_read_at_all(fin, myNzValOffset, MPI_BOTTOM, 1, type,&status);
548 if (err != MPI_SUCCESS) {
549 throw std::logic_error(
"error reading nzval" );
552 MPI_Type_free(&type);
556 for( Int i = 1; i < numColLocal + 1; i++ ){
557 pspmat.colptrLocal[i] = pspmat.colptrLocal[i] - pspmat.colptrLocal[0] + 1;
559 pspmat.colptrLocal[0]=1;
563 MPI_File_close(&fin);
571 inline void ReadDistSparseMatrixFormatted (
const char* filename, DistSparseMatrix<Real>& pspmat, MPI_Comm comm )
574 PushCallStack(
"ReadDistSparseMatrixFormatted");
578 Int mpirank; MPI_Comm_rank(comm, &mpirank);
579 Int mpisize; MPI_Comm_size(comm, &mpisize);
592 throw std::logic_error(
"File cannot be opened!" );
596 fin >> pspmat.size >> dummy;
602 MPI_Bcast(&pspmat.size, 1, MPI_INT, 0, comm);
603 MPI_Bcast(&pspmat.nnz, 1, MPI_INT, 0, comm);
607 IntNumVec colptr(pspmat.size+1);
609 Int* ptr = colptr.Data();
610 for( Int i = 0; i < pspmat.size+1; i++ )
614 MPI_Bcast(colptr.Data(), pspmat.size+1, MPI_INT, 0, comm);
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);
622 numColLocal = numColLocalVec[mpirank];
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;
630 pspmat.nnzLocal = pspmat.colptrLocal[numColLocal] - pspmat.colptrLocal[0];
632 pspmat.rowindLocal.Resize( pspmat.nnzLocal );
633 pspmat.nzvalLocal.Resize ( pspmat.nnzLocal );
640 for( Int ip = 0; ip < mpisize; ip++ ){
641 numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
642 colptr[ip*numColFirst];
644 Int *ptr = buf.Data();
645 for( Int i = 0; i < numRead; i++ ){
649 MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
650 MPI_Send(buf.Data(), numRead, MPI_INT, ip, 1, comm);
653 pspmat.rowindLocal = buf;
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() );
668 pspmat.rowindLocal.Resize( numRead );
669 MPI_Recv( pspmat.rowindLocal.Data(), numRead, MPI_INT, 0, 1, comm, &mpistat );
681 for( Int ip = 0; ip < mpisize; ip++ ){
682 numRead = colptr[ip*numColFirst + numColLocalVec[ip]] -
683 colptr[ip*numColFirst];
685 Real *ptr = buf.Data();
686 for( Int i = 0; i < numRead; i++ ){
690 MPI_Send(&numRead, 1, MPI_INT, ip, 0, comm);
691 MPI_Send(buf.Data(), numRead, MPI_DOUBLE, ip, 1, comm);
694 pspmat.nzvalLocal = buf;
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() );
709 pspmat.nzvalLocal.Resize( numRead );
710 MPI_Recv( pspmat.nzvalLocal.Data(), numRead, MPI_DOUBLE, 0, 1, comm, &mpistat );
729 template<
typename T>
inline void GetDiagonal (
const DistSparseMatrix<T>& A, NumVec<T>& diag )
732 PushCallStack(
"GetDiagonal");
734 Int mpirank, mpisize;
735 MPI_Comm_rank( A.comm, &mpirank );
736 MPI_Comm_size( A.comm, &mpisize );
738 NumVec<T> diagLocal( A.size );
740 diag.Resize( A.size );
743 Int numColFirst = A.size / mpisize;
744 Int firstCol = mpirank * numColFirst;
745 Int numColLocal = A.colptrLocal.m() - 1;
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;
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 );
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() );
768 Int diagIdx = ptr - A.rowindLocal.Data();
769 diagLocal( jcol - 1 ) = A.nzvalLocal( diagIdx );
772 mpi::Allreduce( &diagLocal[0], &diag[0], A.size, MPI_SUM, A.comm );
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