From: Tim K. <tim...@ce...> - 2011-03-10 15:39:13
|
Dear all, I'm trying to write a sparse matrix to a file. SparseMatrix::print() does this, but it prints all the 0 entries, which results in a huge file that is practically useless (besides from exceeding my disk quota). PetscMatrix::print_matlab() is supposed to do better, but PETSc says: [2]PETSC ERROR: ASCII matrix output not allowed for matrices with more than 1024 rows, use binary format instead. PetscMatrix::print_personal() can only print to stdout, which makes it difficult for me to cut out the correct part of the output, in particular if I'm writing more than one matrix. I would like to implement a method, e.g. called SparseMatrix::print_sparse(), which works the same way as SparseMatrix::print() does, but does not print all the 0 entries. It would e.g. rather print one entry per line in the format row column value Any objections or thoughts? -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim...@ce... Universitaetsallee 29 tim...@me... D-28359 Bremen Phone +49-421-218-59246 Germany Fax +49-421-218-59277 |
From: John P. <pet...@cf...> - 2011-03-10 15:53:41
|
On Thu, Mar 10, 2011 at 8:39 AM, Tim Kroeger <tim...@ce...> wrote: > Dear all, > > I'm trying to write a sparse matrix to a file. > > SparseMatrix::print() does this, but it prints all the 0 entries, > which results in a huge file that is practically useless (besides from > exceeding my disk quota). > > PetscMatrix::print_matlab() is supposed to do better, but PETSc says: > > [2]PETSC ERROR: ASCII matrix output not allowed for matrices with more than 1024 rows, use binary format instead. Hmm, I didn't realize the Matlab printing in PETSc had this limitation. I've only used the SparseMatrix printing capability to print test matrices much smaller in size. > PetscMatrix::print_personal() can only print to stdout, which makes it > difficult for me to cut out the correct part of the output, in > particular if I'm writing more than one matrix. I think we should try to fix/update the PetscMatrix implementation of print_personal. I don't think there's any limitation of MatView which prevents its contents from being written to file instead of stdout. > I would like to implement a method, e.g. called > SparseMatrix::print_sparse(), which works the same way as > SparseMatrix::print() does, but does not print all the 0 entries. It > would e.g. rather print one entry per line in the format > > row column value > > Any objections or thoughts? I'd prefer not adding yet another print* function to SparseMatrix if possible. If it turns out that print_personal() can't be fixed up the way you want, I think it'd be easier to pass an additional flag (with default value) to print() to make it do what you want. You'd have to duplicate a lot of code from print() to implement print_sparse() anyway, wouldn't you? -- John |
From: Tim K. <tim...@ce...> - 2011-03-10 16:09:49
|
On Thu, 10 Mar 2011, John Peterson wrote: > On Thu, Mar 10, 2011 at 8:39 AM, Tim Kroeger > <tim...@ce...> wrote: > >> PetscMatrix::print_personal() can only print to stdout, which makes it >> difficult for me to cut out the correct part of the output, in >> particular if I'm writing more than one matrix. > > I think we should try to fix/update the PetscMatrix implementation of > print_personal. I don't think there's any limitation of MatView which > prevents its contents from being written to file instead of stdout. I don't know either. But as far as I understand, print_personal() does different things for different solver packages, and I think it would be a nice thing to have a method that supplies a uniform output format for all packages, and takes the sparsity pattern into account. >> I would like to implement a method, e.g. called >> SparseMatrix::print_sparse(), which works the same way as >> SparseMatrix::print() does, but does not print all the 0 entries. It >> would e.g. rather print one entry per line in the format >> >> row column value >> >> Any objections or thoughts? > > I'd prefer not adding yet another print* function to SparseMatrix if > possible. If it turns out that print_personal() can't be fixed up the > way you want, I think it'd be easier to pass an additional flag (with > default value) to print() to make it do what you want. You'd have to > duplicate a lot of code from print() to implement print_sparse() > anyway, wouldn't you? That's a good point. What about this: SparseMatrix::print(std::ostream& os, const bool sparse=false); By the way, there is another issue with that method, I've just ran into. I did: std::ofstream sMatrix("matrix"); matrix->print(sMatrix); But that opens libMesh::processor_id() instances of the same file at once and writes the matrix to only one of them. Since these files keep overwriting each other, the result is rubbish. This can of course be avoided on user side by if(libMesh::processor_id()==0) { std::ofstream sMatrix("matrix"); matrix->print(sMatrix); } else { matrix->print(std::cout); } but that's not really nice, is it? How are you guys calling print() when you print to a file? Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim...@ce... Universitaetsallee 29 tim...@me... D-28359 Bremen Phone +49-421-218-59246 Germany Fax +49-421-218-59277 |
From: Roy S. <roy...@ic...> - 2011-03-10 16:22:29
|
On Thu, 10 Mar 2011, Tim Kroeger wrote: > I would like to implement a method, e.g. called > SparseMatrix::print_sparse(), which works the same way as > SparseMatrix::print() does, but does not print all the 0 entries. It > would e.g. rather print one entry per line in the format > > row column value > > Any objections or thoughts? Looks like the Matlab sparse matrix format (or at least the representation of it they show to users) - lousy for efficiency but perfect for interoperability. This would be great, thanks! I personally like the print_sparse idea slightly more than John's additional flag idea, but he seems to have a stronger preference the other way so go with that. One question about formatting: do we omit or print entries which are present in the sparsity pattern but have exactly-zero coefficients? I'd lean toward "omit" but could make an argument either way. Nice gotcha (and ugly-but-probably-optimal workaround) you found on the parallel ofstream issue. I usually just print to cout and redirect to file for debugging, so never noticed this before. Not sure how to make your solution any prettier without introducing a new ostream subclass. --- Roy |
From: John P. <pet...@cf...> - 2011-03-10 16:28:15
|
On Thu, Mar 10, 2011 at 9:09 AM, Tim Kroeger <tim...@ce...> wrote: > On Thu, 10 Mar 2011, John Peterson wrote: > >> On Thu, Mar 10, 2011 at 8:39 AM, Tim Kroeger >> <tim...@ce...> wrote: >> >>> PetscMatrix::print_personal() can only print to stdout, which makes it >>> difficult for me to cut out the correct part of the output, in >>> particular if I'm writing more than one matrix. >> >> I think we should try to fix/update the PetscMatrix implementation of >> print_personal. I don't think there's any limitation of MatView which >> prevents its contents from being written to file instead of stdout. > > I don't know either. But as far as I understand, print_personal() does > different things for different solver packages, and I think it would be a > nice thing to have a method that supplies a uniform output format for all > packages, and takes the sparsity pattern into account. Well, print_personal() already takes an ostream& so I don't think we'd be stepping outside the current bounds of its responsibilities to provide it a reasonable execution path for when os != std::cout. For PETSc, I think PETSC_VIEWER_ASCII_COMMON, which uses "a sparse format common among all matrix types" would be the way to go in print_personal to avoid printing the zeros. >>> I would like to implement a method, e.g. called >>> SparseMatrix::print_sparse(), which works the same way as >>> SparseMatrix::print() does, but does not print all the 0 entries. It >>> would e.g. rather print one entry per line in the format >>> >>> row column value >>> >>> Any objections or thoughts? >> >> I'd prefer not adding yet another print* function to SparseMatrix if >> possible. If it turns out that print_personal() can't be fixed up the >> way you want, I think it'd be easier to pass an additional flag (with >> default value) to print() to make it do what you want. You'd have to >> duplicate a lot of code from print() to implement print_sparse() >> anyway, wouldn't you? > > That's a good point. What about this: > > SparseMatrix::print(std::ostream& os, const bool sparse=false); Yeah, that's pretty much what I had in mind. > By the way, there is another issue with that method, I've just ran into. I > did: > > std::ofstream sMatrix("matrix"); > matrix->print(sMatrix); > > But that opens libMesh::processor_id() instances of the same file at once > and writes the matrix to only one of them. Since these files keep > overwriting each other, the result is rubbish. This can of course be > avoided on user side by Are you using a complex-enabled build of LibMesh? The complex version of SparseMatrix::print(), which is defined in the header file, looks like it incorrectly writes on all processors in a parallel run, but the version of SparseMatrix::print() in the C file looks like it should do the right thing in parallel. The Complex version should be wrapped in if (libMesh::processor_id()==0), or, more ideally, rewritten in the style of the C file implementation. -- John |
From: John P. <pet...@cf...> - 2011-03-10 16:43:15
|
On Thu, Mar 10, 2011 at 9:27 AM, John Peterson <pet...@cf...> wrote: > On Thu, Mar 10, 2011 at 9:09 AM, Tim Kroeger > <tim...@ce...> wrote: > > Well, print_personal() already takes an ostream& so I don't think we'd > be stepping outside the current bounds of its responsibilities to > provide it a reasonable execution path for when os != std::cout. For > PETSc, I think PETSC_VIEWER_ASCII_COMMON, which uses "a sparse format > common among all matrix types" would be the way to go in > print_personal to avoid printing the zeros. Actually, print_personal() already handles the sparsity format pretty nicely, doesn't it (at least in PETSc-2.3.3)? I added a print_personal() call in ex4 and the lines look like: row 0: (0, 3.55556e+08) (1, -4.44444e+07) (2, -0.0222223) (3, -4.44444e+07) (4, 8.88889e+07) (5, 0.111111) ... row 1: (0, -4.44444e+07) (1, 3.55556e+08) (2, -0.0666666) (3, -0.0222223) (4, 8.88889e+07) ... ... Is this not already sufficient for your purposes? I mean, once we start sending it to a file instead of stdout of course... -- John |
From: John P. <pet...@cf...> - 2011-03-10 17:48:09
Attachments:
petsc_print_personal.patch
|
Hi Tim, The attached patch makes PetscMatrix::print_personal() print to file. Would you try it out and see if that works for you? -- John |
From: Roy S. <roy...@ic...> - 2011-03-10 18:52:26
|
On Thu, 10 Mar 2011, John Peterson wrote: > The attached patch makes PetscMatrix::print_personal() print to file. > Would you try it out and see if that works for you? Could you change this to use mkstemp to get the temporary file? Running multiple libMesh apps from the same working directory is probably pretty rare, but it happens, and when it does there'd be a race condition for "libmesh_temp_petsc_matrix.081702". That temp-file-to-ofstream copy is a wonderful, horrifying hack, by the way. --- Roy |
From: John P. <pet...@cf...> - 2011-03-10 20:23:42
|
On Thu, Mar 10, 2011 at 11:52 AM, Roy Stogner <roy...@ic...> wrote: > > On Thu, 10 Mar 2011, John Peterson wrote: > >> The attached patch makes PetscMatrix::print_personal() print to file. >> Would you try it out and see if that works for you? > > Could you change this to use mkstemp to get the temporary file? > Running multiple libMesh apps from the same working directory is > probably pretty rare, but it happens, and when it does there'd be a > race condition for "libmesh_temp_petsc_matrix.081702". Yeah, you probably saw I had the mkstemp call in there. I commented it out since I wasn't sure what would happen in parallel... But I just read in the docs of PetscViewerASCIIOpen() that "If a multiprocessor communicator is used (such as PETSC_COMM_WORLD), then only the first processor in the group opens the file. All other processors send their data to the first processor to print." so I think it's probably safe. > That temp-file-to-ofstream copy is a wonderful, horrifying hack, by > the way. I'm pretty sure it's legal C++, so please let me know if you find out otherwise... -- John |
From: Tim K. <tim...@ce...> - 2011-03-11 08:25:35
|
On Thu, 10 Mar 2011, John Peterson wrote: > The attached patch makes PetscMatrix::print_personal() print to file. > Would you try it out and see if that works for you? Your patch makes it even more difficult to write to a file in parallel, because my trick to open the file on one processor and use std::cout on the others confuses your method in that way that you create your custom PetscViewer on the first processor but call MatView(...,PETSC_VIEWER_STDOUT_SELF) on all other processors. This causes PETSc to crash. Anyway, I somehow dislike the idea of a temporary file. Also, I still think that a method doing the thing independently of the solver package would be desirable. I have meanwhile implemented my print_sparse() method, see attachment. Changing it to the other API (with bool argument) would be very easy. I agree that the complex version of print() (in the header file) is wrong. Since I have never used libMesh with complex numbers and don't plan to do this at the moment, I will leave this issue to somebody else. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim...@ce... Universitaetsallee 29 tim...@me... D-28359 Bremen Phone +49-421-218-59246 Germany Fax +49-421-218-59277 |
From: John P. <pet...@cf...> - 2011-03-11 08:49:10
|
On Fri, Mar 11, 2011 at 1:25 AM, Tim Kroeger <tim...@ce...> wrote: > On Thu, 10 Mar 2011, John Peterson wrote: > >> The attached patch makes PetscMatrix::print_personal() print to file. >> Would you try it out and see if that works for you? > > Your patch makes it even more difficult to write to a file in parallel, > because my trick to open the file on one processor and use std::cout on the > others confuses your method in that way that you create your custom > PetscViewer on the first processor but call > MatView(...,PETSC_VIEWER_STDOUT_SELF) on all other processors. This causes > PETSc to crash. If you have been using a hack (passing a file on one processor and std::cout on the others) then you should stop using that hack before you try out the new patch... -- John |
From: Tim K. <tim...@ce...> - 2011-03-11 09:04:36
|
On Fri, 11 Mar 2011, John Peterson wrote: > On Fri, Mar 11, 2011 at 1:25 AM, Tim Kroeger > <tim...@ce...> wrote: >> On Thu, 10 Mar 2011, John Peterson wrote: >> >>> The attached patch makes PetscMatrix::print_personal() print to file. >>> Would you try it out and see if that works for you? >> >> Your patch makes it even more difficult to write to a file in parallel, >> because my trick to open the file on one processor and use std::cout on the >> others confuses your method in that way that you create your custom >> PetscViewer on the first processor but call >> MatView(...,PETSC_VIEWER_STDOUT_SELF) on all other processors. This causes >> PETSc to crash. > > If you have been using a hack (passing a file on one processor and > std::cout on the others) then you should stop using that hack before > you try out the new patch... But what is the correct way to use your method in parallel? std::ofstream f("filename"); matrix.print_personal(f); can definitely not be correct, no matter what print_personal() does, since it opens the same file several times, and the files overwrite each other in a non-predictable way. The more I think about this, the more I come to the conclusion that the API of getting a stream on all processors is misleading. I mean, what if the stream is e.g. a std::stringstream? Then the user might want to get the result into the stream on *all* processors, which neiter your nor my method correctly performs. What would be a better API? Either a similar API as in PETSc, where either a file name is used (and the method is responsible for opening this file only once) or some special value meaning std::cout. (This disallows to use a std::stringstream, but disallowing something seems better to me than allowing it but behaving unexpectedly.) Or we might want to do print(std::ostream*=NULL); and explicitly allow that the argument differs between processors, where the meaning would be that the result is printed on all processors that do not use NULL. However, this requires additional communication. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim...@ce... Universitaetsallee 29 tim...@me... D-28359 Bremen Phone +49-421-218-59246 Germany Fax +49-421-218-59277 |
From: Tim K. <tim...@ce...> - 2011-03-14 10:14:17
Attachments:
patch
|
On Thu, 10 Mar 2011, John Peterson wrote: > On Thu, Mar 10, 2011 at 9:09 AM, Tim Kroeger > <tim...@ce...> wrote: > >> What about this: >> >> SparseMatrix::print(std::ostream& os, const bool sparse=false); > > Yeah, that's pretty much what I had in mind. I've implemented this now, see attachment. It works for me, and I would like to check it in, unless there are still objections. Note that for the complex case, I left the (wrong) implementation for the non-sparse case unchanged and added libmesh_not_implemented() for the sparse case. Everything remains downwards compatible that way, and who ever wants to fix the complex case can then add or not add the sparse case there. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim...@ce... Universitaetsallee 29 tim...@me... D-28359 Bremen Phone +49-421-218-59246 Germany Fax +49-421-218-59277 |
From: Tim K. <tim...@ce...> - 2011-03-16 09:25:14
|
On Mon, 14 Mar 2011, Tim Kroeger wrote: > On Thu, 10 Mar 2011, John Peterson wrote: > >> On Thu, Mar 10, 2011 at 9:09 AM, Tim Kroeger >> <tim...@ce...> wrote: >> >>> What about this: >>> >>> SparseMatrix::print(std::ostream& os, const bool sparse=false); >> >> Yeah, that's pretty much what I had in mind. > > I've implemented this now, see attachment. It works for me, and I would like > to check it in, unless there are still objections. Since there seem to have been no objections, I have commited this now. Best Regards, Tim -- Dr. Tim Kroeger CeVis -- Center of Complex Systems and Visualization University of Bremen tim...@ce... Universitaetsallee 29 tim...@me... D-28359 Bremen Phone +49-421-218-59246 Germany Fax +49-421-218-59277 |