[petsc-dev] [Fwd: Re: [petsc-users] PETSc recommended visualization packages]
Ethan Coon
ecoon at lanl.gov
Tue Jul 19 12:45:57 CDT 2011
Taking this to dev.
Attached is a patch for PetscBinaryRead, which now has been tested on
Mat, IS, and Vec. Supports returning dense numpy.matrix, sparse CSR,
and scipy.sparse matrices (all have been tested).
Also cc'ed Ata directly since I'm not sure he's on dev. Ata -- could
you please test the attached PetscBinaryRead.py for me on your system to
make sure I've correctly handled your version of numpy? Note that the
return format has changed from what you have now. Thanks.
PetscBinaryWrite.py when I get to it...
Thanks,
Ethan
-------- Forwarded Message --------
From: Lisandro Dalcin <dalcinl at gmail.com>
Reply-to: PETSc users list <petsc-users at mcs.anl.gov>
To: PETSc users list <petsc-users at mcs.anl.gov>
Subject: Re: [petsc-users] PETSc recommended visualization packages
Date: Thu, 14 Jul 2011 19:34:42 -0300
On 14 July 2011 18:33, Ethan Coon <ecoon at lanl.gov> wrote:
> On Thu, 2011-07-14 at 17:47 -0300, Lisandro Dalcin wrote:
>
>>
>> 1) readMatDense() is wrong. The I,J,V arrays are the CSR structure,
>> not the COO (coordinate) format, Then you cannot simply:
>>
>
> Ok, thanks, I wasn't sure of which format you were using,
Well, PETSc's AIJ matrices are CSR in memory (but not in disk, binary
files store row counts instead of row pointers)
>>
>> 2) It would be nice to provie a 'scipy' Mat format returning a
>> 'scipy.sparse.csr_matrix' instance. A possible implementation would
>> be:
>>
>
> Agreed, I considered doing that as the default for sparse, but wasn't
> sure what others would want. Providing both makes sense.
>
Note however that this would require a SciPy installation.
PS: Sorry for not being more helpful and contributing actual code, but
I'm really busy.
--
------------------------------------
Ethan Coon
Post-Doctoral Researcher
Applied Mathematics - T-5
Los Alamos National Laboratory
505-665-8289
http://www.ldeo.columbia.edu/~ecoon/
------------------------------------
-------------- next part --------------
# HG changeset patch
# Parent 2ae64ecc368c681373625f390ef9360746a3edf6
diff -r 2ae64ecc368c bin/python/PetscBinaryRead.py
--- a/bin/python/PetscBinaryRead.py Mon Jul 18 15:11:19 2011 -0500
+++ b/bin/python/PetscBinaryRead.py Tue Jul 19 11:31:13 2011 -0600
@@ -1,3 +1,18 @@
+"""PetscBinaryRead
+===============
+
+Provides
+ 1. PETSc-named objects Vec, Mat, and IS that inherit numpy.ndarray
+ 2. A function to read these objects from PETSc binary files.
+
+The standard usage of this module should look like:
+
+ >>> import PetscBinaryRead
+ >>> objects = PetscBinaryRead.readBinaryFile('filename')
+
+See readBinaryFile.__doc__
+"""
+
import numpy as np
import types
@@ -10,6 +25,23 @@
class DoneWithFile(Exception): pass
+class Vec(np.ndarray):
+ """Vec represented as 1D numpy array"""
+ pass
+
+class MatDense(np.matrix):
+ """Mat represented as 2D numpy array"""
+ pass
+
+class MatSparse(tuple):
+ """Mat represented as CSR tuple ((M, N), (row, col, val))"""
+ def __repr__(self):
+ return 'MatSparse: %s'%super(MatSparse, self).__repr__()
+
+class IS(np.ndarray):
+ """IS represented as 1D numpy array"""
+ pass
+
def readVec(fh):
"""Reads a PETSc Vec from a binary file handle, returning just the data."""
@@ -18,7 +50,9 @@
vals = np.fromfile(fh, dtype=ScalarType, count=nz)
except MemoryError:
raise IOError('Inconsistent or invalid Vec data in file')
- return vals
+ if (len(vals) is 0):
+ raise IOError('Inconsistent or invalid Vec data in file')
+ return vals.view(Vec)
def readMatSparse(fh):
"""Reads a PETSc Mat, returning a sparse representation of the data.
@@ -42,11 +76,13 @@
assert I[-1] == nz
#
J = np.fromfile(fh, dtype=IntType, count=nz)
+ assert len(J) == nz
V = np.fromfile(fh, dtype=ScalarType, count=nz)
- except AssertionError, MemoryError:
+ assert len(V) == nz
+ except (AssertionError, MemoryError, IndexError):
raise IOError('Inconsistent or invalid Mat data in file')
#
- return (M, N), (I, J, V)
+ return MatSparse(((M, N), (I, J, V)))
def readMatDense(fh):
"""Reads a PETSc Mat, returning a dense represention of the data."""
@@ -59,15 +95,24 @@
np.cumsum(rownz, out=I[1:])
assert I[-1] == nz
#
- J = np.fromfile(fh, dtype=IntType, count=nz)
+ J = np.fromfile(fh, dtype=IntType, count=nz)
+ assert len(J) == nz
V = np.fromfile(fh, dtype=ScalarType, count=nz)
- mat = np.zeros((M,N), dtype=ScalarType)
- for i,j,v in zip(I,J,V):
- mat[i,j] = v
- except AssertionError, MemoryError:
+ assert len(V) == nz
+
+ except (AssertionError, MemoryError, IndexError):
raise IOError('Inconsistent or invalid Mat data in file')
#
- return mat
+ mat = np.zeros((M,N), dtype=ScalarType)
+ for row in range(M):
+ rstart, rend = I[row:row+2]
+ mat[row, J[rstart:rend]] = V[rstart:rend]
+ return mat.view(MatDense)
+
+def readMatSciPy(fh):
+ from scipy.sparse import csr_matrix
+ (M, N), (I, J, V) = readMatSparse(fh)
+ return csr_matrix((V, J, I), shape=(M, N))
def readMat(fh, mattype='sparse'):
"""Reads a PETSc Mat from binary file handle.
@@ -80,6 +125,8 @@
return readMatSparse(fh)
elif mattype == 'dense':
return readMatDense(fh)
+ elif mattype == 'scipy.sparse':
+ return readMatSciPy(fh)
else:
raise RuntimeError('Invalid matrix type requested: choose sparse/dense')
@@ -88,21 +135,25 @@
try:
nz = np.fromfile(fh, dtype=IntType, count=1)[0]
v = np.fromfile(fh, dtype=IntType, count=nz)
- except MemoryError:
+ assert len(v) == nz
+ except (MemoryError,IndexError):
raise IOError('Inconsistent or invalid IS data in file')
- return v
+ return v.view(IS)
def readBinaryFile(fid, mattype='sparse'):
- """Reads a PETSc binary file, returning a tuple of objects contained in the file.
+ """Reads a PETSc binary file, returning a tuple of the contained objects.
objects = readBinaryFile(fid, mattype='sparse')
Input:
- fid : either file handle to open binary file, or filename
- mattype : ['sparse'] Read matrices as 'sparse' (row, col, val) or 'dense'.
+ fid : either file handle to an open binary file, or filename.
+ mattype :
+ 'sparse': Return matrices as raw CSR: (M, N), (row, col, val).
+ 'dense': Return matrices as MxN numpy arrays.
+ 'scipy.sparse': Return matrices as scipy.sparse objects.
Output:
- objects : tuple of objects representing the data.
+ objects : tuple of objects representing the data in numpy arrays.
"""
close = False
@@ -116,19 +167,19 @@
# read header
try:
header = np.fromfile(fid, dtype=IntType, count=1)[0]
- except MemoryError:
+ except (MemoryError, IndexError):
raise DoneWithFile
try:
objecttype = _classid[header]
except KeyError:
- raise IOError('Invalid PetscObject CLASSID or CLASSID not yet implemented')
+ raise IOError('Invalid PetscObject CLASSID or object not implemented for python')
if objecttype == 'Vec':
- objects.append(('Vec', readVec(fid)))
+ objects.append(readVec(fid))
elif objecttype == 'IS':
- objects.append(('IS', readIS(fid)))
+ objects.append(readIS(fid))
elif objecttype == 'Mat':
- objects.append(('Mat', readMat(fid,mattype)))
+ objects.append(readMat(fid,mattype))
elif objecttype == 'Bag':
raise NotImplementedError('Bag Reader not yet implemented')
except DoneWithFile:
@@ -141,10 +192,9 @@
if __name__ == '__main__':
import sys
- objects = readBinaryFile(sys.argv[1])
- for otype, data in objects:
- print 'Read a', otype
- print data
- print '\n'
+ petsc_objects = readBinaryFile(sys.argv[1])
+ for petsc_obj in petsc_objects:
+ print 'Read a', petsc_obj
+ print ''
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PetscBinaryRead.py
Type: text/x-python
Size: 5951 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20110719/275a722e/attachment.py>
More information about the petsc-dev
mailing list