[Scipy-svn] r2496 - in trunk/Lib/sparse: . sparsetools tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Jan 6 00:34:54 CST 2007
Author: wnbell
Date: 2007-01-06 00:34:41 -0600 (Sat, 06 Jan 2007)
New Revision: 2496
Added:
trunk/Lib/sparse/sparsetools/complex_ops.h
trunk/Lib/sparse/sparsetools/numpy.i
trunk/Lib/sparse/sparsetools/sparsetools.h
trunk/Lib/sparse/sparsetools/sparsetools.i
trunk/Lib/sparse/sparsetools/sparsetools.py
trunk/Lib/sparse/sparsetools/sparsetools_wrap.cxx
Removed:
trunk/Lib/sparse/sparsetools.txt
trunk/Lib/sparse/sparsetools/sparsetools.pyf.src
trunk/Lib/sparse/sparsetools/spblas.f.src
trunk/Lib/sparse/sparsetools/spconv.f.src
Modified:
trunk/Lib/sparse/setup.py
trunk/Lib/sparse/sparse.py
trunk/Lib/sparse/tests/test_sparse.py
Log:
Rewrite of sparsetools. additional sparse unittests, and
changes to sparse.py to interface w/ the new sparsetools
Modified: trunk/Lib/sparse/setup.py
===================================================================
--- trunk/Lib/sparse/setup.py 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/setup.py 2007-01-06 06:34:41 UTC (rev 2496)
@@ -15,10 +15,16 @@
# sources = [join('sparsekit','*.f')]
# )
- sources = ['spblas.f.src','spconv.f.src','sparsetools.pyf.src']
+## sources = ['spblas.f.src','spconv.f.src','sparsetools.pyf.src']
+## sources = [join('sparsetools',x) for x in sources]
+
+## config.add_extension('sparsetools',
+## sources = sources,
+## )
+ sources = ['sparsetools_wrap.cxx','sparsetools.py']
sources = [join('sparsetools',x) for x in sources]
- config.add_extension('sparsetools',
+ config.add_extension('_sparsetools',
sources = sources,
)
Modified: trunk/Lib/sparse/sparse.py
===================================================================
--- trunk/Lib/sparse/sparse.py 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparse.py 2007-01-06 06:34:41 UTC (rev 2496)
@@ -2,6 +2,7 @@
Original code by Travis Oliphant.
Modified and extended by Ed Schofield and Robert Cimrman.
+Revision of sparsetools by Nathan Bell
"""
from numpy import zeros, isscalar, real, imag, asarray, asmatrix, matrix, \
@@ -23,14 +24,6 @@
ALLOCSIZE=1000
NZMAX = 100
-_coerce_rules = {('f', 'f'):'f', ('f', 'd'):'d', ('f', 'F'):'F',
- ('f', 'D'):'D', ('d', 'f'):'d', ('d', 'd'):'d',
- ('d', 'F'):'D', ('d', 'D'):'D', ('F', 'f'):'F',
- ('F', 'd'):'D', ('F', 'F'):'F', ('F', 'D'):'D',
- ('D', 'f'):'D', ('D', 'd'):'d', ('D', 'F'):'D',
- ('D', 'D'):'D'}
-_transtabl = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
-_itranstabl = {'s':'f', 'd':'d', 'c':'F', 'z':'D'}
# The formats that we might potentially understand.
_formats = {'csc':[0, "Compressed Sparse Column"],
@@ -55,6 +48,17 @@
'und':[19, "Undefined"]
}
+
+
+#these four don't appear to be used anymore (Jan 6th, 2007)
+_coerce_rules = {('f', 'f'):'f', ('f', 'd'):'d', ('f', 'F'):'F',
+ ('f', 'D'):'D', ('d', 'f'):'d', ('d', 'd'):'d',
+ ('d', 'F'):'D', ('d', 'D'):'D', ('F', 'f'):'F',
+ ('F', 'd'):'D', ('F', 'F'):'F', ('F', 'D'):'D',
+ ('D', 'f'):'D', ('D', 'd'):'d', ('D', 'F'):'D',
+ ('D', 'D'):'D'}
+_transtabl = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
+_itranstabl = {'s':'f', 'd':'d', 'c':'F', 'z':'D'}
def _convert_data(data1, data2, newtype):
if data1.dtype.char != newtype:
data1 = data1.astype(newtype)
@@ -62,6 +66,9 @@
data2 = data2.astype(newtype)
return data1, data2
+
+
+
class spmatrix(object):
""" This class provides a base class for all sparse matrices. It
cannot be instantiated. Most of the work is provided by subclasses.
@@ -476,25 +483,8 @@
# Use a double array as the source (but leave it alone)
s = s*1.0
if (rank(s) == 2):
- M, N = s.shape
- dtype = s.dtype
- func = getattr(sparsetools, _transtabl[dtype.char]+'fulltocsc')
- ierr = irow = jcol = 0
- nnz = (s != 0.0).sum()
- a = zeros((nnz,), self.dtype)
- rowa = zeros((nnz,), intc)
- ptra = zeros((N+1,), intc)
- while 1:
- a, rowa, ptra, irow, jcol, ierr = \
- func(s, a, rowa, ptra, irow, jcol, ierr)
- if (ierr == 0): break
- nnz = nnz + ALLOCSIZE
- a = resize1d(a, nnz)
- rowa = resize1d(rowa, nnz)
- self.data = a
- self.rowind = rowa
- self.indptr = ptra
- self.shape = (M, N)
+ self.shape = s.shape
+ self.indptr,self.rowind,self.data = sparsetools.densetocsr(s.shape[1],s.shape[0],s.T)
else:
raise ValueError, "dense array must have rank 1 or 2"
elif isspmatrix(arg1):
@@ -513,9 +503,7 @@
self.indptr = s.indptr
elif isinstance(s, csr_matrix):
self.shape = s.shape
- func = getattr(sparsetools, s.ftype+'transp')
- self.data, self.rowind, self.indptr = \
- func(s.shape[1], s.data, s.colind, s.indptr)
+ self.indptr,self.rowind,self.data = sparsetools.csrtocsc(s.shape[0],s.shape[1],s.indptr,s.colind,s.data)
else:
temp = s.tocsc()
self.data = temp.data
@@ -609,6 +597,10 @@
raise ValueError, \
"Last value of index list should be less than "\
"the size of data list"
+ if (self.rowind.dtype != numpy.intc):
+ raise TypeError, "rowind indices must be of type numpy.intc"
+ if (self.indptr.dtype != numpy.intc):
+ raise TypeError, "inptr indices must be of type numpy.intc"
self.nnz = nnz
self.nzmax = nzmax
self.dtype = self.data.dtype
@@ -641,18 +633,10 @@
ocs = other.tocsc()
if (ocs.shape != self.shape):
raise ValueError, "inconsistent shapes"
- dtypechar = _coerce_rules[(self.dtype.char, ocs.dtype.char)]
- nnz1, nnz2 = self.nnz, ocs.nnz
- if (nnz1 == 0): nnz1 = 1
- if (nnz2 == 0): nnz2 = 1
- data1, data2 = _convert_data(self.data[:nnz1], ocs.data[:nnz2], dtypechar)
- func = getattr(sparsetools, _transtabl[dtypechar]+'cscadd')
- c, rowc, ptrc, ierr = func(data1, self.rowind[:nnz1], self.indptr,
- data2, ocs.rowind[:nnz2], ocs.indptr)
- if ierr:
- raise RuntimeError, "ran out of space"
- M, N = self.shape
- return csc_matrix((c, rowc, ptrc), dims=(M, N))
+ indptr,rowind,data = sparsetools.cscplcsc(self.shape[0],self.shape[1], \
+ self.indptr,self.rowind,self.data,\
+ ocs.indptr,ocs.rowind,ocs.data)
+ return csc_matrix((data,rowind,indptr),self.shape)
elif isdense(other):
# Convert this matrix to a dense matrix and add them.
return self.todense() + other
@@ -667,18 +651,10 @@
ocs = other.tocsc()
if (ocs.shape != self.shape):
raise ValueError, "inconsistent shapes"
- dtypechar = _coerce_rules[(self.dtype.char, ocs.dtype.char)]
- nnz1, nnz2 = self.nnz, ocs.nnz
- if (nnz1 == 0): nnz1=1
- if (nnz2 == 0): nnz2=1
- data1, data2 = _convert_data(self.data[:nnz1], ocs.data[:nnz2], dtypechar)
- func = getattr(sparsetools, _transtabl[dtypechar]+'cscadd')
- c, rowc, ptrc, ierr = func(data1, self.rowind[:nnz1], self.indptr,
- data2, ocs.rowind[:nnz2], ocs.indptr)
- if ierr:
- raise RuntimeError, "ran out of space"
- M, N = self.shape
- return csc_matrix((c, rowc, ptrc), dims=(M, N))
+ indptr,rowind,data = sparsetools.cscplcsc(self.shape[0],self.shape[1], \
+ self.indptr,self.rowind,self.data,\
+ ocs.indptr,ocs.rowind,ocs.data)
+ return csc_matrix((data,rowind,indptr),self.shape)
elif isdense(other):
# Convert this matrix to a dense matrix and add them
return other + self.todense()
@@ -731,15 +707,10 @@
ocs = other.tocsc()
if (ocs.shape != self.shape):
raise ValueError, "inconsistent shapes"
- dtypechar = _coerce_rules[(self.dtype.char, ocs.dtype.char)]
- nnz1, nnz2 = self.nnz, ocs.nnz
- data1, data2 = _convert_data(self.data[:nnz1], ocs.data[:nnz2], dtypechar)
- func = getattr(sparsetools, _transtabl[dtypechar]+'cscmul')
- c, rowc, ptrc, ierr = func(data1, self.rowind[:nnz1], self.indptr, data2, ocs.rowind[:nnz2], ocs.indptr)
- if ierr:
- raise RuntimeError, "ran out of space"
- M, N = self.shape
- return csc_matrix((c, rowc, ptrc), dims=(M, N))
+ indptr,rowind,data = sparsetools.cscelmulcsc(self.shape[0],self.shape[1],\
+ self.indptr,self.rowind,self.data,\
+ ocs.indptr,ocs.rowind,ocs.data)
+ return csc_matrix((data,rowind,indptr),(self.shape[0],ocs.shape[1]))
def transpose(self, copy=False):
M, N = self.shape
@@ -801,8 +772,9 @@
# being created on-the-fly like with dense matrix objects.
#if len(other) != self.shape[1]:
# raise ValueError, "dimension mismatch"
- func = getattr(sparsetools, self.ftype+'cscmux')
- y = func(self.data, self.rowind, self.indptr, other, self.shape[0])
+ oth = numpy.ravel(other)
+ y = sparsetools.cscmux(self.shape[0],self.shape[1],\
+ self.indptr,self.rowind,self.data,oth)
if isinstance(other, matrix):
y = asmatrix(y)
# If 'other' was an (nx1) column vector, transpose the result
@@ -821,13 +793,12 @@
# being created on-the-fly like with dense matrix objects.
#if len(other) != self.shape[0]:
# raise ValueError, "dimension mismatch"
- func = getattr(sparsetools, self.ftype+'csrmux')
+ oth = numpy.ravel(other)
if conjugate:
cd = conj(self.data)
else:
cd = self.data
- y = func(cd, self.rowind, self.indptr, other)
-
+ y = sparsetool.csrmux(self.shape[1],self.shape[0],self.indptr,self.rowind,cd,oth)
if isinstance(other, matrix):
y = asmatrix(y)
# In the (unlikely) event that this matrix is 1x1 and 'other' was an
@@ -846,51 +817,10 @@
K2, N = other.shape
if (K1 != K2):
raise ValueError, "shape mismatch error"
- a, rowa, ptra = self.data, self.rowind, self.indptr
- if isinstance(other, csr_matrix):
- other._check()
- dtypechar = _coerce_rules[(self.dtype.char, other.dtype.char)]
- ftype = _transtabl[dtypechar]
- func = getattr(sparsetools, ftype+'cscmucsr')
- b = other.data
- rowb = other.colind
- ptrb = other.indptr
- elif isinstance(other, csc_matrix):
- other._check()
- dtypechar = _coerce_rules[(self.dtype.char, other.dtype.char)]
- ftype = _transtabl[dtypechar]
- func = getattr(sparsetools, ftype+'cscmucsc')
- b = other.data
- rowb = other.rowind
- ptrb = other.indptr
- else:
- other = other.tocsc()
- dtypechar = _coerce_rules[(self.dtype.char, other.dtype.char)]
- ftype = _transtabl[dtypechar]
- func = getattr(sparsetools, ftype+'cscmucsc')
- b = other.data
- rowb = other.rowind
- ptrb = other.indptr
- a, b = _convert_data(a, b, dtypechar)
- ptrc = zeros((N+1,), intc)
- nnzc = 2*max(ptra[-1], ptrb[-1])
- c = zeros((nnzc,), dtypechar)
- rowc = zeros((nnzc,), intc)
- ierr = irow = kcol = 0
- while True: # loop in case first call runs out of memory.
- c, rowc, ptrc, irow, kcol, ierr = func(M, a, rowa, ptra, b,
- rowb, ptrb, c, rowc,
- ptrc, irow, kcol, ierr)
- if (ierr==0): break
- # otherwise we were too small and must resize arrays
- # calculations continue where they left off...
- print "Resizing...", kcol, irow, ierr
- percent_to_go = 1- (1.0*kcol*M + irow) / (N*M)
- newnnzc = int(ceil((1+percent_to_go)*nnzc))
- c = resize1d(c, newnnzc)
- rowc = resize1d(rowc, newnnzc)
- nnzc = newnnzc
- return csc_matrix((c, rowc, ptrc), dims=(M, N))
+ other = other.tocsc()
+ indptr,rowind,data = sparsetools.cscmucsc(M,N,self.indptr,self.rowind,self.data,\
+ other.indptr,other.rowind,other.data)
+ return csc_matrix((data,rowind,indptr),(M,N))
elif isdense(other):
# This is SLOW! We need a more efficient implementation
# of sparse * dense matrix multiplication!
@@ -915,9 +845,12 @@
col = N + col
if not (0<=row<M) or not (0<=col<N):
raise IndexError, "index out of bounds"
- func = getattr(sparsetools, self.ftype+'cscgetel')
- ind, val = func(self.data, self.rowind, self.indptr, row, col)
- return val
+ #this was implemented in fortran before - is there a noticable performance advangate?
+ indxs = numpy.where(row == self.rowind[self.indptr[col]:self.indptr[col+1]])
+ if len(indxs[0]) == 0:
+ return 0
+ else:
+ return self.data[self.indptr[col]:self.indptr[col+1]][indxs[0]]
elif isintlike(key):
# Was: return self.data[key]
# If this is allowed, it should return the relevant row, as for
@@ -931,7 +864,6 @@
if isinstance(key, tuple):
row = key[0]
col = key[1]
- func = getattr(sparsetools, self.ftype+'cscsetel')
M, N = self.shape
if (row < 0):
row = M + row
@@ -946,12 +878,30 @@
if (row >= M):
M = row+1
self.shape = (M, N)
- nzmax = self.nzmax
- if (nzmax < self.nnz+1): # need more room
- alloc = max(1, self.allocsize)
- self.data = resize1d(self.data, nzmax + alloc)
- self.rowind = resize1d(self.rowind, nzmax + alloc)
- func(self.data, self.rowind, self.indptr, row, col, val)
+
+ indxs = numpy.where(row == self.rowind[self.indptr[col]:self.indptr[col+1]])
+ if len(indxs[0]) == 0:
+ #value not present
+ nzmax = self.nzmax
+ if (nzmax < self.nnz+1): # need more room
+ alloc = max(1, self.allocsize)
+ self.data = resize1d(self.data, nzmax + alloc)
+ self.rowind = resize1d(self.rowind, nzmax + alloc)
+
+ newindex = self.indptr[col]
+ self.data[newindex+1:] = self.data[newindex:-1]
+ self.rowind[newindex+1:] = self.rowind[newindex:-1]
+
+ self.data[newindex] = val
+ self.rowind[newindex] = row
+ self.indptr[col+1:] += 1
+
+ elif len(indxs[0]) == 1:
+ #value already present
+ self.data[self.indptr[col]:self.indptr[col+1]][indxs[0]] = val
+ else:
+ raise IndexError, "row index occurs more than once"
+
self._check()
else:
# We should allow slices here!
@@ -966,31 +916,18 @@
if stop <= start:
raise ValueError, "slice width must be >= 1"
startind = -1
- # Locate the first element in self.data
- # (could be made more efficient with a binary search)
- for ind in xrange(self.indptr[j], self.indptr[j+1]):
- if self.rowind[ind] >= start:
- startind = ind
- break
- if startind == -1:
- # Col has only zeros
- return csc_matrix((stop-start, 1), dtype=self.dtype, \
- nzmax=min(NZMAX, stop-start))
+
+ indices = []
- stopind = self.indptr[j+1]
- # Locate the last+1 index into self.data
- # (could be made more efficient with a binary search)
- for ind in xrange(startind, self.indptr[j+1]):
- if self.rowind[ind] >= stop:
- stopind = ind
- break
-
- data = self.data[startind : stopind]
- rowind = self.rowind[startind : stopind] - start
- indptr = numpy.array([0, stopind - startind])
- new = csc_matrix((data, rowind, indptr), dims=(stop-start, 1), \
- dtype=self.dtype)
- return new
+ for ind in xrange(self.indptr[j], self.indptr[j+1]):
+ if self.rowind[ind] >= start and self.rowind[ind] < stop:
+ indices.append(ind)
+
+ rowind = self.rowind[indices] - start
+ data = self.data[indices]
+ indptr = numpy.array([0,len(indices)])
+ return csc_matrix((data, rowind, indptr), dims=(stop-start, 1), \
+ dtype=self.dtype)
def rowcol(self, ind):
row = self.rowind[ind]
@@ -1008,19 +945,16 @@
return new
def tocoo(self):
- if self.nnz == 0:
- return coo_matrix(None, dims=self.shape, dtype=self.dtype)
- else:
- func = getattr(sparsetools, self.ftype+"csctocoo")
- data, row, col = func(self.data, self.rowind, self.indptr)
- return coo_matrix((data, (row, col)), dims=self.shape)
+ rows,cols,data = sparsetools.csctocoo(self.shape[0],self.shape[1],\
+ self.indptr,self.rowind,self.data)
+ return coo_matrix((data,(rows,cols)),self.shape)
def tocsr(self):
- return self.tocoo().tocsr()
-
+ indptr,colind,data = sparsetools.csctocsr(self.shape[0],self.shape[1],self.indptr,self.rowind,self.data)
+ return csr_matrix((data,colind,indptr),self.shape)
+
def toarray(self):
- func = getattr(sparsetools, self.ftype+'csctofull')
- return func(self.shape[0], self.data, self.rowind, self.indptr)
+ return self.tocsr().toarray()
def prune(self):
""" Remove empty space after all non-zero elements.
@@ -1097,11 +1031,6 @@
self.data = s.data
self.colind = s.colind
self.indptr = s.indptr
- elif isinstance(s, csc_matrix):
- self.shape = s.shape
- func = getattr(sparsetools, s.ftype+'transp')
- self.data, self.colind, self.indptr = \
- func(s.shape[1], s.data, s.rowind, s.indptr)
else:
try:
temp = s.tocsr()
@@ -1196,6 +1125,10 @@
raise ValueError, \
"last value of index list should be less than "\
"the size of data list"
+ if (self.colind.dtype != numpy.intc):
+ raise TypeError, "colind indices must be of type numpy.intc"
+ if (self.indptr.dtype != numpy.intc):
+ raise TypeError, "inptr indices must be of type numpy.intc"
self.nnz = nnz
self.nzmax = nzmax
self.dtype = self.data.dtype
@@ -1226,19 +1159,12 @@
raise NotImplementedError, 'adding a scalar to a CSR matrix ' \
'is not yet supported'
elif isspmatrix(other):
- ocs = other.tocsr()
- if (ocs.shape != self.shape):
+ other = other.tocsr()
+ if (other.shape != self.shape):
raise ValueError, "inconsistent shapes"
-
- dtypechar = _coerce_rules[(self.dtype.char, ocs.dtype.char)]
- data1, data2 = _convert_data(self.data, ocs.data, dtypechar)
- func = getattr(sparsetools, _transtabl[dtypechar]+'cscadd')
- c, colc, ptrc, ierr = func(data1, self.colind, self.indptr, data2,
- ocs.colind, ocs.indptr)
- if ierr:
- raise RuntimeError, "ran out of space"
- M, N = self.shape
- return csr_matrix((c, colc, ptrc), dims=(M, N))
+ indptr,colind,data = sparsetools.csrplcsr(self.shape[0],other.shape[1],\
+ self.indptr,self.colind,self.data,other.indptr,other.colind,other.data)
+ return csr_matrix((data,colind,indptr),self.shape)
elif isdense(other):
# Convert this matrix to a dense matrix and add them.
return self.todense() + other
@@ -1288,18 +1214,13 @@
new.ftype = _transtabl[new.dtype.char]
return new
elif isspmatrix(other):
- ocs = other.tocsr()
- if (ocs.shape != self.shape):
+ other = other.tocsr()
+ if (other.shape != self.shape):
raise ValueError, "inconsistent shapes"
- dtypechar = _coerce_rules[(self.dtype.char, ocs.dtype.char)]
- data1, data2 = _convert_data(self.data, ocs.data, dtypechar)
- func = getattr(sparsetools, _transtabl[dtypechar]+'cscmul')
- c, colc, ptrc, ierr = func(data1, self.colind, self.indptr, data2,
- ocs.colind, ocs.indptr)
- if ierr:
- raise RuntimeError, "ran out of space"
- M, N = self.shape
- return csr_matrix((c, colc, ptrc), dims=(M, N))
+ indptr,colind,data = sparsetools.csrelmulcsr(self.shape[0],other.shape[1],\
+ self.indptr,self.colind,self.data,\
+ other.indptr,other.colind,other.data)
+ return csr_matrix((data,colind,indptr),(self.shape[0],other.shape[1]))
else:
raise TypeError, "unsupported type for sparse matrix power"
@@ -1349,9 +1270,8 @@
# being created on-the-fly like dense matrix objects can.
#if len(other) != self.shape[1]:
# raise ValueError, "dimension mismatch"
- func = getattr(sparsetools, self.ftype+'csrmux')
- y = func(self.data, self.colind, self.indptr, other)
-
+ oth = numpy.ravel(other)
+ y = sparsetools.csrmux(self.shape[0],self.shape[1],self.indptr,self.colind,self.data,oth)
if isinstance(other, matrix):
y = asmatrix(y)
# If 'other' was an (nx1) column vector, transpose the result
@@ -1370,8 +1290,9 @@
cd = conj(self.data)
else:
cd = self.data
- y = func(cd, self.colind, self.indptr, other, self.shape[1])
+ oth = numpy.ravel(other)
+ y = sparsetools.cscmux(self.shape[0],self.shape[1],self.indptr,self.rowind,cd,oth)
if isinstance(other, matrix):
y = asmatrix(y)
# In the (unlikely) event that this matrix is 1x1 and 'other' was an
@@ -1387,69 +1308,11 @@
a, rowa, ptra = self.data, self.colind, self.indptr
if (K1 != K2):
raise ValueError, "shape mismatch error"
- if isinstance(other, csc_matrix):
- other._check()
- dtypechar = _coerce_rules[(self.dtype.char, other.dtype.char)]
- ftype = _transtabl[dtypechar]
- func = getattr(sparsetools, ftype+'csrmucsc')
- b = other.data
- colb = other.rowind
- ptrb = other.indptr
- out = 'csc'
- firstarg = ()
- elif isinstance(other, csr_matrix):
- other._check()
- dtypechar = _coerce_rules[(self.dtype.char, other.dtype.char)]
- ftype = _transtabl[dtypechar]
- func = getattr(sparsetools, ftype+'cscmucsc')
- b, colb, ptrb = a, rowa, ptra
- a, rowa, ptra = other.data, other.colind, other.indptr
- out = 'csr'
- firstarg = (N,)
- else:
- other = other.tocsc()
- dtypechar = _coerce_rules[(self.dtype.char, other.dtype.char)]
- ftype = _transtabl[dtypechar]
- func = getattr(sparsetools, ftype+'csrmucsc')
- b = other.data
- colb = other.rowind
- ptrb = other.indptr
- out = 'csc'
- firstarg = ()
- a, b = _convert_data(a, b, dtypechar)
- newshape = (M, N)
- if out == 'csr':
- ptrc = zeros((M+1,), intc)
- else:
- ptrc = zeros((N+1,), intc)
- nnzc = 2*max(ptra[-1], ptrb[-1])
- # Avoid an infinite loop when multiplying by a matrix with
- # only zeros
- if nnzc == 0:
- if out == 'csr':
- return csr_matrix(newshape, dtype=dtypechar)
- else:
- return csc_matrix(newshape, dtype=dtypechar)
- c = zeros((nnzc,), dtypechar)
- rowc = zeros((nnzc,), intc)
- ierr = irow = kcol = 0
- while 1:
- args = firstarg+(a, rowa, ptra, b, colb, ptrb, c, rowc, ptrc, irow,
- kcol, ierr)
- c, rowc, ptrc, irow, kcol, ierr = func(*args)
- if (ierr==0): break
- # otherwise we were too small and must resize
- percent_to_go = 1- (1.0*kcol) / N
- newnnzc = int(ceil((1+percent_to_go)*nnzc))
- c = resize1d(c, newnnzc)
- rowc = resize1d(rowc, newnnzc)
- nnzc = newnnzc
-
- if out == 'csr':
- # Note: 'rowc' is deliberate
- return csr_matrix((c, rowc, ptrc), dims=(M, N))
- else:
- return csc_matrix((c, rowc, ptrc), dims=(M, N))
+ other = other.tocsr()
+ indptr,colind,data = sparsetools.csrmucsr(self.shape[0],other.shape[1],\
+ self.indptr,self.colind,self.data,\
+ other.indptr,other.colind,other.data)
+ return csr_matrix((data,colind,indptr),(self.shape[0],other.shape[1]))
elif isdense(other):
# This is SLOW! We need a more efficient implementation
# of sparse * dense matrix multiplication!
@@ -1474,9 +1337,12 @@
col = N + col
if not (0<=row<M) or not (0<=col<N):
raise IndexError, "index out of bounds"
- func = getattr(sparsetools, self.ftype+'cscgetel')
- ind, val = func(self.data, self.colind, self.indptr, col, row)
- return val
+
+ indxs = numpy.where(col == self.colind[self.indptr[row]:self.indptr[row+1]])
+ if len(indxs[0]) == 0:
+ return 0
+ else:
+ return self.data[self.indptr[row]:self.indptr[row+1]][indxs[0]]
elif isintlike(key):
return self[key, :]
else:
@@ -1492,37 +1358,23 @@
if stop <= start:
raise ValueError, "slice width must be >= 1"
startind = -1
- # Locate the first element in self.data
- # (could be made more efficient with a binary search)
- for ind in xrange(self.indptr[i], self.indptr[i+1]):
- if self.colind[ind] >= start:
- startind = ind
- break
- if startind == -1:
- # Row has only zeros
- return csr_matrix((1, stop-start), dtype=self.dtype, \
- nzmax=min(NZMAX, stop-start))
+
+ indices = []
- stopind = self.indptr[i+1]
- # Locate the last+1 index into self.data
- # (could be made more efficient with a binary search)
- for ind in xrange(startind, self.indptr[i+1]):
- if self.colind[ind] >= stop:
- stopind = ind
- break
-
- data = self.data[startind : stopind]
- colind = self.colind[startind : stopind] - start
- indptr = numpy.array([0, stopind - startind])
- new = csr_matrix((data, colind, indptr), dims=(1, stop-start), \
- dtype=self.dtype)
- return new
+ for ind in xrange(self.indptr[i], self.indptr[i+1]):
+ if self.colind[ind] >= start and self.colind[ind] < stop:
+ indices.append(ind)
+
+ colind = self.colind[indices] - start
+ data = self.data[indices]
+ indptr = numpy.array([0,len(indices)])
+ return csr_matrix((data, colind, indptr), dims=(1,stop-start), \
+ dtype=self.dtype)
def __setitem__(self, key, val):
if isinstance(key, tuple):
row = key[0]
col = key[1]
- func = getattr(sparsetools, self.ftype+'cscsetel')
M, N = self.shape
if (row < 0):
row = M + row
@@ -1534,19 +1386,33 @@
self.indptr = resize1d(self.indptr, row+2)
self.indptr[M+1:] = self.indptr[M]
M = row+1
- elif (row < 0):
- row = M - row
if (col >= N):
N = col+1
- elif (col < 0):
- col = N - col
self.shape = (M, N)
- nzmax = self.nzmax
- if (nzmax < self.nnz+1): # need more room
- alloc = max(1, self.allocsize)
- self.data = resize1d(self.data, nzmax + alloc)
- self.colind = resize1d(self.colind, nzmax + alloc)
- func(self.data, self.colind, self.indptr, col, row, val)
+
+ indxs = numpy.where(col == self.colind[self.indptr[row]:self.indptr[row+1]])
+ if len(indxs[0]) == 0:
+ #value not present
+ nzmax = self.nzmax
+ if (nzmax < self.nnz+1): # need more room
+ alloc = max(1, self.allocsize)
+ self.data = resize1d(self.data, nzmax + alloc)
+ self.colind = resize1d(self.colind, nzmax + alloc)
+
+ newindex = self.indptr[row]
+ self.data[newindex+1:] = self.data[newindex:-1]
+ self.colind[newindex+1:] = self.colind[newindex:-1]
+
+ self.data[newindex] = val
+ self.colind[newindex] = col
+ self.indptr[row+1:] += 1
+
+ elif len(indxs[0]) == 1:
+ #value already present
+ self.data[self.indptr[row]:self.indptr[row+1]][indxs[0]] = val
+ else:
+ raise IndexError, "row index occurs more than once"
+
self._check()
else:
# We should allow slices here!
@@ -1567,20 +1433,19 @@
return self
def tocoo(self):
- if self.nnz == 0:
- return coo_matrix(None, dims=self.shape, dtype=self.dtype)
- else:
- func = getattr(sparsetools, self.ftype+"csctocoo")
- data, col, row = func(self.data, self.colind, self.indptr)
- return coo_matrix((data, (row, col)), dims=self.shape)
+ rows,cols,data = sparsetools.csrtocoo(self.shape[0],self.shape[1],\
+ self.indptr,self.colind,self.data)
+ return coo_matrix((data,(rows,cols)),self.shape)
def tocsc(self):
- return self.tocoo().tocsc()
+ indptr,rowind,data = sparsetools.csrtocsc(self.shape[0],self.shape[1],self.indptr,self.colind,self.data)
+ return csc_matrix((data,rowind,indptr),self.shape)
+
def toarray(self):
- func = getattr(sparsetools, self.ftype+'csctofull')
- s = func(self.shape[1], self.data, self.colind, self.indptr)
- return s.transpose()
+ data = numpy.zeros(self.shape,self.data.dtype)
+ sparsetools.csrtodense(self.shape[0],self.shape[1],self.indptr,self.colind,self.data,data)
+ return data
def prune(self):
""" Eliminate non-zero entries, leaving space for at least
@@ -1648,6 +1513,7 @@
assert M == int(M) and M > 0
assert N == int(N) and N > 0
self.shape = (int(M), int(N))
+ return
except (TypeError, ValueError, AssertionError):
raise TypeError, "dimensions must be a 2-tuple of positive"\
" integers"
@@ -2057,6 +1923,7 @@
base = dok_matrix()
ext = dok_matrix()
indx = int((columns == 1))
+ N = len(cols_or_rows)
if indx:
for key in self:
num = searchsorted(cols_or_rows, key[1])
@@ -2132,6 +1999,7 @@
ikey0 = int(key[0])
ikey1 = int(key[1])
if ikey0 != current_row:
+ N = ikey0-current_row
row_ptr[current_row+1:ikey0+1] = k
current_row = ikey0
data[k] = dict.__getitem__(self, key)
@@ -2164,6 +2032,7 @@
ikey0 = int(key[0])
ikey1 = int(key[1])
if ikey1 != current_col:
+ N = ikey1-current_col
col_ptr[current_col+1:ikey1+1] = k
current_col = ikey1
data[k] = self[key]
@@ -2272,6 +2141,10 @@
if (nnz != len(self.row)) or (nnz != len(self.col)):
raise ValueError, "row, column, and data array must all be "\
"the same length"
+ if (self.row.dtype != numpy.intc):
+ raise TypeError, "row indices must be of type numpy.intc"
+ if (self.col.dtype != numpy.intc):
+ raise TypeError, "col indices must be of type numpy.intc"
self.nnz = nnz
self.ftype = _transtabl.get(self.dtype.char,'')
@@ -2304,24 +2177,17 @@
if self.nnz == 0:
return csc_matrix(self.shape, dtype=self.dtype)
else:
- func = getattr(sparsetools, self.ftype+"cootocsc")
- data, row, col = self._normalize()
- a, rowa, ptra, ierr = func(self.shape[1], data, row, col)
- if ierr:
- raise RuntimeError, "error in conversion"
- return csc_matrix((a, rowa, ptra), dims=self.shape)
+ indptr,rowind,data = sparsetools.cootocsc(self.shape[0],self.shape[1],self.size,self.row,self.col,self.data)
+ return csc_matrix((data,rowind,indptr),self.shape)
+
def tocsr(self):
if self.nnz == 0:
return csr_matrix(self.shape, dtype=self.dtype)
else:
- func = getattr(sparsetools, self.ftype+"cootocsc")
- data, row, col = self._normalize(rowfirst=True)
- a, cola, ptra, ierr = func(self.shape[0], data, col, row)
- if ierr:
- raise RuntimeError, "error in conversion"
- return csr_matrix((a, cola, ptra), dims=self.shape)
-
+ indptr,colind,data = sparsetools.cootocsr(self.shape[0],self.shape[1],self.size,self.row,self.col,self.data)
+ return csr_matrix((data,colind,indptr),self.shape)
+
def tocoo(self, copy=False):
if copy:
return self.copy()
@@ -2811,15 +2677,10 @@
if not hasattr(offsets, '__len__' ):
offsets = (offsets,)
offsets = array(offsets, copy=False)
- mtype = diags.dtype.char
assert(len(offsets) == diags.shape[0])
- # set correct diagonal to csc conversion routine for this type
- diagfunc = eval('sparsetools.'+_transtabl[mtype]+'diatocsc')
- a, rowa, ptra, ierr = diagfunc(M, N, diags, offsets)
- if ierr:
- raise RuntimeError, "ran out of space"
- return csc_matrix((a, rowa, ptra), dims=(M, N))
-
+ indptr,rowind,data = sparsetools.spdiags(M,N,len(offsets),offsets,diags)
+ return csc_matrix((data,rowind,indptr),(M,N))
+
def spidentity(n, dtype='d'):
"""
spidentity( n ) returns the identity matrix of shape (n, n) stored
Added: trunk/Lib/sparse/sparsetools/complex_ops.h
===================================================================
--- trunk/Lib/sparse/sparsetools/complex_ops.h 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparsetools/complex_ops.h 2007-01-06 06:34:41 UTC (rev 2496)
@@ -0,0 +1,159 @@
+#ifndef COMPLEX_OPS_H
+#define COMPLEX_OPS_H
+
+/*
+ * Functions to handle arithmetic operations on NumPy complex values
+ */
+
+
+
+
+/*
+ * Addition
+ */
+inline npy_cfloat operator+(const npy_cfloat& A, const npy_cfloat& B){
+ npy_cfloat result;
+ result.real = A.real + B.real;
+ result.imag = A.imag + B.imag;
+ return result;
+}
+inline npy_cdouble operator+(const npy_cdouble& A, const npy_cdouble& B){
+ npy_cdouble result;
+ result.real = A.real + B.real;
+ result.imag = A.imag + B.imag;
+ return result;
+}
+inline npy_clongdouble operator+(const npy_clongdouble& A, const npy_clongdouble& B){
+ npy_clongdouble result;
+ result.real = A.real + B.real;
+ result.imag = A.imag + B.imag;
+ return result;
+}
+
+inline npy_cfloat& operator+=(npy_cfloat& A, const npy_cfloat& B){
+ A.real += B.real;
+ A.imag += B.imag;
+ return A;
+}
+inline npy_cdouble& operator+=(npy_cdouble& A, const npy_cdouble& B){
+ A.real += B.real;
+ A.imag += B.imag;
+ return A;
+}
+inline npy_clongdouble& operator+=(npy_clongdouble& A, const npy_clongdouble& B){
+ A.real += B.real;
+ A.imag += B.imag;
+ return A;
+}
+
+/*
+ * Multiplication
+ */
+inline npy_cfloat operator*(const npy_cfloat& A, const npy_cfloat& B){
+ npy_cfloat result;
+ result.real = A.real * B.real - A.imag * B.imag;
+ result.imag = A.real * B.imag + A.imag * B.real;
+ return result;
+}
+inline npy_cdouble operator*(const npy_cdouble& A, const npy_cdouble& B){
+ npy_cdouble result;
+ result.real = A.real * B.real - A.imag * B.imag;
+ result.imag = A.real * B.imag + A.imag * B.real;
+ return result;
+}
+inline npy_clongdouble operator*(const npy_clongdouble& A, const npy_clongdouble& B){
+ npy_clongdouble result;
+ result.real = A.real * B.real - A.imag * B.imag;
+ result.imag = A.real * B.imag + A.imag * B.real;
+ return result;
+}
+
+inline npy_cfloat& operator*=(npy_cfloat& A, const npy_cfloat& B){
+ npy_float temp = A.real * B.real - A.imag * B.imag;
+ A.imag = A.real * B.imag + A.imag * B.real;
+ A.real = temp;
+ return A;
+}
+inline npy_cdouble& operator*=(npy_cdouble& A, const npy_cdouble& B){
+ npy_double temp = A.real * B.real - A.imag * B.imag;
+ A.imag = A.real * B.imag + A.imag * B.real;
+ A.real = temp;
+ return A;
+}
+inline npy_clongdouble& operator*=(npy_clongdouble& A, const npy_clongdouble& B){
+ npy_longdouble temp = A.real * B.real - A.imag * B.imag;
+ A.imag = A.real * B.imag + A.imag * B.real;
+ A.real = temp;
+ return A;
+}
+
+/*
+ * Equality (complex==complex)
+ */
+inline bool operator==(const npy_cfloat& A, const npy_cfloat& B){
+ return A.real == B.real && A.imag == B.imag;
+}
+inline bool operator==(const npy_cdouble& A, const npy_cdouble& B){
+ return A.real == B.real && A.imag == B.imag;
+}
+inline bool operator==(const npy_clongdouble& A, const npy_clongdouble& B){
+ return A.real == B.real && A.imag == B.imag;
+}
+
+inline bool operator!=(const npy_cfloat& A, const npy_cfloat& B){
+ return A.real != B.real || A.imag != B.imag;
+}
+inline bool operator!=(const npy_cdouble& A, const npy_cdouble& B){
+ return A.real != B.real || A.imag != B.imag;
+}
+inline bool operator!=(const npy_clongdouble& A, const npy_clongdouble& B){
+ return A.real != B.real || A.imag != B.imag;
+}
+
+/*
+ * Equality (complex==scalar)
+ */
+inline bool operator==(const npy_cfloat& A, const npy_float& B){
+ return A.real == B && A.imag == 0;
+}
+inline bool operator==(const npy_cdouble& A, const npy_double& B){
+ return A.real == B && A.imag == 0;
+}
+inline bool operator==(const npy_clongdouble& A, const npy_longdouble& B){
+ return A.real == B && A.imag == 0;
+}
+
+inline bool operator!=(const npy_cfloat& A, const npy_float& B){
+ return A.real != B || A.imag != 0;
+}
+inline bool operator!=(const npy_cdouble& A, const npy_double& B){
+ return A.real != B || A.imag != 0;
+}
+inline bool operator!=(const npy_clongdouble& A, const npy_longdouble& B){
+ return A.real != B || A.imag != 0;
+}
+
+/*
+ * Equality (scalar==complex)
+ */
+inline bool operator==(const npy_float& A, const npy_cfloat& B){
+ return A == B.real && 0 == B.imag;
+}
+inline bool operator==(const npy_double& A, const npy_cdouble& B){
+ return A == B.real && 0 == B.imag;
+}
+inline bool operator==(const npy_longdouble& A, const npy_clongdouble& B){
+ return A == B.real && 0 == B.imag;
+}
+
+inline bool operator!=(const npy_float& A, const npy_cfloat& B){
+ return A != B.real || 0 != B.imag;
+}
+inline bool operator!=(const npy_double& A, const npy_cdouble& B){
+ return A != B.real || 0 != B.imag;
+}
+inline bool operator!=(const npy_longdouble& A, const npy_clongdouble& B){
+ return A != B.real || 0 != B.imag;
+}
+
+#endif
Added: trunk/Lib/sparse/sparsetools/numpy.i
===================================================================
--- trunk/Lib/sparse/sparsetools/numpy.i 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparsetools/numpy.i 2007-01-06 06:34:41 UTC (rev 2496)
@@ -0,0 +1,563 @@
+/* -*- C -*- (not really, but good for syntax highlighting) */
+%{
+#ifndef SWIG_FILE_WITH_INIT
+# define NO_IMPORT_ARRAY
+#endif
+#include "stdio.h"
+#include <numpy/arrayobject.h>
+
+/* The following code originally appeared in enthought/kiva/agg/src/numeric.i,
+ * author unknown. It was translated from C++ to C by John Hunter. Bill
+ * Spotz has modified it slightly to fix some minor bugs, add some comments
+ * and some functionality.
+ */
+
+/* Macros to extract array attributes.
+ */
+#define is_array(a) ((a) && PyArray_Check((PyArrayObject *)a))
+#define array_type(a) (int)(PyArray_TYPE(a))
+#define array_dimensions(a) (((PyArrayObject *)a)->nd)
+#define array_size(a,i) (((PyArrayObject *)a)->dimensions[i])
+#define array_is_contiguous(a) (PyArray_ISCONTIGUOUS(a))
+
+/* Given a PyObject, return a string describing its type.
+ */
+char* pytype_string(PyObject* py_obj) {
+ if (py_obj == NULL ) return "C NULL value";
+ if (PyCallable_Check(py_obj)) return "callable" ;
+ if (PyString_Check( py_obj)) return "string" ;
+ if (PyInt_Check( py_obj)) return "int" ;
+ if (PyFloat_Check( py_obj)) return "float" ;
+ if (PyDict_Check( py_obj)) return "dict" ;
+ if (PyList_Check( py_obj)) return "list" ;
+ if (PyTuple_Check( py_obj)) return "tuple" ;
+ if (PyFile_Check( py_obj)) return "file" ;
+ if (PyModule_Check( py_obj)) return "module" ;
+ if (PyInstance_Check(py_obj)) return "instance" ;
+
+ return "unkown type";
+}
+
+/* Given a Numeric typecode, return a string describing the type.
+ */
+char* typecode_string(int typecode) {
+ char* type_names[20] = {"char","unsigned byte","byte","short",
+ "unsigned short","int","unsigned int","long",
+ "float","double","complex float","complex double",
+ "object","ntype","unkown"};
+ return type_names[typecode];
+}
+
+/* Make sure input has correct numeric type. Allow character and byte
+ * to match. Also allow int and long to match.
+ */
+int type_match(int actual_type, int desired_type) {
+ return PyArray_EquivTypenums(actual_type, desired_type);
+}
+
+/* Given a PyObject pointer, cast it to a PyArrayObject pointer if
+ * legal. If not, set the python error string appropriately and
+ * return NULL./
+ */
+PyArrayObject* obj_to_array_no_conversion(PyObject* input, int typecode) {
+ PyArrayObject* ary = NULL;
+ if (is_array(input) && (typecode == PyArray_NOTYPE ||
+ PyArray_EquivTypenums(array_type(input),
+ typecode))) {
+ ary = (PyArrayObject*) input;
+ }
+ else if is_array(input) {
+ char* desired_type = typecode_string(typecode);
+ char* actual_type = typecode_string(array_type(input));
+ PyErr_Format(PyExc_TypeError,
+ "Array of type '%s' required. Array of type '%s' given",
+ desired_type, actual_type);
+ ary = NULL;
+ }
+ else {
+ char * desired_type = typecode_string(typecode);
+ char * actual_type = pytype_string(input);
+ PyErr_Format(PyExc_TypeError,
+ "Array of type '%s' required. A %s was given",
+ desired_type, actual_type);
+ ary = NULL;
+ }
+ return ary;
+}
+
+/* Convert the given PyObject to a Numeric array with the given
+ * typecode. On Success, return a valid PyArrayObject* with the
+ * correct type. On failure, the python error string will be set and
+ * the routine returns NULL.
+ */
+PyArrayObject* obj_to_array_allow_conversion(PyObject* input, int typecode,
+ int* is_new_object)
+{
+ PyArrayObject* ary = NULL;
+ PyObject* py_obj;
+ if (is_array(input) && (typecode == PyArray_NOTYPE || type_match(array_type(input),typecode))) {
+ ary = (PyArrayObject*) input;
+ *is_new_object = 0;
+ }
+ else {
+ py_obj = PyArray_FromObject(input, typecode, 0, 0);
+ /* If NULL, PyArray_FromObject will have set python error value.*/
+ ary = (PyArrayObject*) py_obj;
+ *is_new_object = 1;
+ }
+ return ary;
+}
+
+/* Given a PyArrayObject, check to see if it is contiguous. If so,
+ * return the input pointer and flag it as not a new object. If it is
+ * not contiguous, create a new PyArrayObject using the original data,
+ * flag it as a new object and return the pointer.
+ */
+PyArrayObject* make_contiguous(PyArrayObject* ary, int* is_new_object,
+ int min_dims, int max_dims)
+{
+ PyArrayObject* result;
+ if (array_is_contiguous(ary)) {
+ result = ary;
+ *is_new_object = 0;
+ }
+ else {
+ result = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*)ary,
+ array_type(ary),
+ min_dims,
+ max_dims);
+ *is_new_object = 1;
+ }
+ return result;
+}
+
+/* Convert a given PyObject to a contiguous PyArrayObject of the
+ * specified type. If the input object is not a contiguous
+ * PyArrayObject, a new one will be created and the new object flag
+ * will be set.
+ */
+PyArrayObject* obj_to_array_contiguous_allow_conversion(PyObject* input,
+ int typecode,
+ int* is_new_object) {
+ int is_new1 = 0;
+ int is_new2 = 0;
+ PyArrayObject* ary2;
+ PyArrayObject* ary1 = obj_to_array_allow_conversion(input, typecode,
+ &is_new1);
+ if (ary1) {
+ ary2 = make_contiguous(ary1, &is_new2, 0, 0);
+ if ( is_new1 && is_new2) {
+ Py_DECREF(ary1);
+ }
+ ary1 = ary2;
+ }
+ *is_new_object = is_new1 || is_new2;
+ return ary1;
+}
+
+/* Test whether a python object is contiguous. If array is
+ * contiguous, return 1. Otherwise, set the python error string and
+ * return 0.
+ */
+int require_contiguous(PyArrayObject* ary) {
+ int contiguous = 1;
+ if (!array_is_contiguous(ary)) {
+ PyErr_SetString(PyExc_TypeError, "Array must be contiguous. A discontiguous array was given");
+ contiguous = 0;
+ }
+ return contiguous;
+}
+
+/* Require the given PyArrayObject to have a specified number of
+ * dimensions. If the array has the specified number of dimensions,
+ * return 1. Otherwise, set the python error string and return 0.
+ */
+int require_dimensions(PyArrayObject* ary, int exact_dimensions) {
+ int success = 1;
+ if (array_dimensions(ary) != exact_dimensions) {
+ PyErr_Format(PyExc_TypeError,
+ "Array must be have %d dimensions. Given array has %d dimensions",
+ exact_dimensions, array_dimensions(ary));
+ success = 0;
+ }
+ return success;
+}
+
+/* Require the given PyArrayObject to have one of a list of specified
+ * number of dimensions. If the array has one of the specified number
+ * of dimensions, return 1. Otherwise, set the python error string
+ * and return 0.
+ */
+int require_dimensions_n(PyArrayObject* ary, int* exact_dimensions, int n) {
+ int success = 0;
+ int i;
+ char dims_str[255] = "";
+ char s[255];
+ for (i = 0; i < n && !success; i++) {
+ if (array_dimensions(ary) == exact_dimensions[i]) {
+ success = 1;
+ }
+ }
+ if (!success) {
+ for (i = 0; i < n-1; i++) {
+ sprintf(s, "%d, ", exact_dimensions[i]);
+ strcat(dims_str,s);
+ }
+ sprintf(s, " or %d", exact_dimensions[n-1]);
+ strcat(dims_str,s);
+ PyErr_Format(PyExc_TypeError,
+ "Array must be have %s dimensions. Given array has %d dimensions",
+ dims_str, array_dimensions(ary));
+ }
+ return success;
+}
+
+/* Require the given PyArrayObject to have a specified shape. If the
+ * array has the specified shape, return 1. Otherwise, set the python
+ * error string and return 0.
+ */
+int require_size(PyArrayObject* ary, int* size, int n) {
+ int i;
+ int success = 1;
+ int len;
+ char desired_dims[255] = "[";
+ char s[255];
+ char actual_dims[255] = "[";
+ for(i=0; i < n;i++) {
+ if (size[i] != -1 && size[i] != array_size(ary,i)) {
+ success = 0;
+ }
+ }
+ if (!success) {
+ for (i = 0; i < n; i++) {
+ if (size[i] == -1) {
+ sprintf(s, "*,");
+ }
+ else
+ {
+ sprintf(s, "%d,", size[i]);
+ }
+ strcat(desired_dims,s);
+ }
+ len = strlen(desired_dims);
+ desired_dims[len-1] = ']';
+ for (i = 0; i < n; i++) {
+ sprintf(s, "%d,", array_size(ary,i));
+ strcat(actual_dims,s);
+ }
+ len = strlen(actual_dims);
+ actual_dims[len-1] = ']';
+ PyErr_Format(PyExc_TypeError,
+ "Array must be have shape of %s. Given array has shape of %s",
+ desired_dims, actual_dims);
+ }
+ return success;
+}
+/* End John Hunter translation (with modifications by Bill Spotz) */
+
+%}
+
+/* TYPEMAP_IN macros
+ *
+ * This family of typemaps allows pure input C arguments of the form
+ *
+ * (type* IN_ARRAY1, int DIM1)
+ * (type* IN_ARRAY2, int DIM1, int DIM2)
+ *
+ * where "type" is any type supported by the Numeric module, to be
+ * called in python with an argument list of a single array (or any
+ * python object that can be passed to the Numeric.array constructor
+ * to produce an arrayof te specified shape). This can be applied to
+ * a existing functions using the %apply directive:
+ *
+ * %apply (double* IN_ARRAY1, int DIM1) {double* series, int length}
+ * %apply (double* IN_ARRAY2, int DIM1, int DIM2) {double* mx, int rows, int cols}
+ * double sum(double* series, int length);
+ * double max(double* mx, int rows, int cols);
+ *
+ * or with
+ *
+ * double sum(double* IN_ARRAY1, int DIM1);
+ * double max(double* IN_ARRAY2, int DIM1, int DIM2);
+ */
+
+/* One dimensional input arrays */
+%define TYPEMAP_IN1(type,typecode)
+ %typemap(in) type* IN_ARRAY1 (PyArrayObject* array=NULL, int is_new_object) {
+ int size[1] = {-1};
+ array = obj_to_array_contiguous_allow_conversion($input, typecode, &is_new_object);
+ if (!array || !require_dimensions(array,1) || !require_size(array,size,1)) SWIG_fail;
+ $1 = (type*) array->data;
+}
+%typemap(freearg) type* IN_ARRAY1 {
+ if (is_new_object$argnum && array$argnum) Py_DECREF(array$argnum);
+}
+%enddef
+
+/* Define concrete examples of the TYPEMAP_IN1 macros */
+TYPEMAP_IN1(char, PyArray_CHAR )
+TYPEMAP_IN1(unsigned char, PyArray_UBYTE )
+TYPEMAP_IN1(signed char, PyArray_SBYTE )
+TYPEMAP_IN1(short, PyArray_SHORT )
+TYPEMAP_IN1(int, PyArray_INT )
+TYPEMAP_IN1(long, PyArray_LONG )
+TYPEMAP_IN1(float, PyArray_FLOAT )
+TYPEMAP_IN1(double, PyArray_DOUBLE)
+TYPEMAP_IN1(long double, PyArray_LONGDOUBLE)
+TYPEMAP_IN1(npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_IN1(npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_IN1(npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_IN1(const char, PyArray_CHAR )
+TYPEMAP_IN1(const unsigned char, PyArray_UBYTE )
+TYPEMAP_IN1(const signed char, PyArray_SBYTE )
+TYPEMAP_IN1(const short, PyArray_SHORT )
+TYPEMAP_IN1(const int, PyArray_INT )
+TYPEMAP_IN1(const long, PyArray_LONG )
+TYPEMAP_IN1(const float, PyArray_FLOAT )
+TYPEMAP_IN1(const double, PyArray_DOUBLE)
+TYPEMAP_IN1(const long double, PyArray_LONGDOUBLE)
+TYPEMAP_IN1(const npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_IN1(const npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_IN1(const npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_IN1(PyObject, PyArray_OBJECT)
+
+#undef TYPEMAP_IN1
+
+
+
+
+ /* Two dimensional input arrays */
+%define TYPEMAP_IN2(type,typecode)
+ %typemap(in) (type* IN_ARRAY2)
+ (PyArrayObject* array=NULL, int is_new_object) {
+ int size[2] = {-1,-1};
+ array = obj_to_array_contiguous_allow_conversion($input, typecode, &is_new_object);
+ if (!array || !require_dimensions(array,2) || !require_size(array,size,1)) SWIG_fail;
+ $1 = (type*) array->data;
+}
+%typemap(freearg) (type* IN_ARRAY2) {
+ if (is_new_object$argnum && array$argnum) Py_DECREF(array$argnum);
+}
+%enddef
+
+/* Define concrete examples of the TYPEMAP_IN2 macros */
+TYPEMAP_IN2(char, PyArray_CHAR )
+TYPEMAP_IN2(unsigned char, PyArray_UBYTE )
+TYPEMAP_IN2(signed char, PyArray_SBYTE )
+TYPEMAP_IN2(short, PyArray_SHORT )
+TYPEMAP_IN2(int, PyArray_INT )
+TYPEMAP_IN2(long, PyArray_LONG )
+TYPEMAP_IN2(float, PyArray_FLOAT )
+TYPEMAP_IN2(double, PyArray_DOUBLE)
+TYPEMAP_IN2(long double, PyArray_LONGDOUBLE)
+TYPEMAP_IN2(npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_IN2(npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_IN2(npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_IN2(const char, PyArray_CHAR )
+TYPEMAP_IN2(const unsigned char, PyArray_UBYTE )
+TYPEMAP_IN2(const signed char, PyArray_SBYTE )
+TYPEMAP_IN2(const short, PyArray_SHORT )
+TYPEMAP_IN2(const int, PyArray_INT )
+TYPEMAP_IN2(const long, PyArray_LONG )
+TYPEMAP_IN2(const float, PyArray_FLOAT )
+TYPEMAP_IN2(const double, PyArray_DOUBLE)
+TYPEMAP_IN2(const long double, PyArray_LONGDOUBLE)
+TYPEMAP_IN2(const npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_IN2(const npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_IN2(const npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_IN2(PyObject, PyArray_OBJECT)
+
+#undef TYPEMAP_IN2
+
+
+/* TYPEMAP_INPLACE macros
+ *
+ * This family of typemaps allows input/output C arguments of the form
+ *
+ * (type* INPLACE_ARRAY1, int DIM1)
+ * (type* INPLACE_ARRAY2, int DIM1, int DIM2)
+ *
+ * where "type" is any type supported by the Numeric module, to be
+ * called in python with an argument list of a single contiguous
+ * Numeric array. This can be applied to an existing function using
+ * the %apply directive:
+ *
+ * %apply (double* INPLACE_ARRAY1, int DIM1) {double* series, int length}
+ * %apply (double* INPLACE_ARRAY2, int DIM1, int DIM2) {double* mx, int rows, int cols}
+ * void negate(double* series, int length);
+ * void normalize(double* mx, int rows, int cols);
+ *
+ *
+ * or with
+ *
+ * void sum(double* INPLACE_ARRAY1, int DIM1);
+ * void sum(double* INPLACE_ARRAY2, int DIM1, int DIM2);
+ */
+
+ /* One dimensional input/output arrays */
+%define TYPEMAP_INPLACE1(type,typecode)
+%typemap(in) (type* INPLACE_ARRAY) (PyArrayObject* temp=NULL) {
+ int i;
+ temp = obj_to_array_no_conversion($input,typecode);
+ if (!temp || !require_contiguous(temp)) SWIG_fail;
+ $1 = (type*) temp->data;
+}
+%enddef
+
+/* Define concrete examples of the TYPEMAP_INPLACE1 macro */
+TYPEMAP_INPLACE1(char, PyArray_CHAR )
+TYPEMAP_INPLACE1(unsigned char, PyArray_UBYTE )
+TYPEMAP_INPLACE1(signed char, PyArray_SBYTE )
+TYPEMAP_INPLACE1(short, PyArray_SHORT )
+TYPEMAP_INPLACE1(int, PyArray_INT )
+TYPEMAP_INPLACE1(long, PyArray_LONG )
+TYPEMAP_INPLACE1(float, PyArray_FLOAT )
+TYPEMAP_INPLACE1(double, PyArray_DOUBLE)
+TYPEMAP_INPLACE1(long double, PyArray_LONGDOUBLE)
+TYPEMAP_INPLACE1(npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_INPLACE1(npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_INPLACE1(npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_INPLACE1(const char, PyArray_CHAR )
+TYPEMAP_INPLACE1(const unsigned char, PyArray_UBYTE )
+TYPEMAP_INPLACE1(const signed char, PyArray_SBYTE )
+TYPEMAP_INPLACE1(const short, PyArray_SHORT )
+TYPEMAP_INPLACE1(const int, PyArray_INT )
+TYPEMAP_INPLACE1(const long, PyArray_LONG )
+TYPEMAP_INPLACE1(const float, PyArray_FLOAT )
+TYPEMAP_INPLACE1(const double, PyArray_DOUBLE)
+TYPEMAP_INPLACE1(const long double, PyArray_LONGDOUBLE)
+TYPEMAP_INPLACE1(const npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_INPLACE1(const npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_INPLACE1(const npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_INPLACE1(PyObject, PyArray_OBJECT)
+
+#undef TYPEMAP_INPLACE1
+
+
+
+
+
+ /* Two dimensional input/output arrays */
+%define TYPEMAP_INPLACE2(type,typecode)
+ %typemap(in) (type* INPLACE_ARRAY2) (PyArrayObject* temp=NULL) {
+ temp = obj_to_array_no_conversion($input,typecode);
+ if (!temp || !require_contiguous(temp)) SWIG_fail;
+ $1 = (type*) temp->data;
+}
+%enddef
+
+/* Define concrete examples of the TYPEMAP_INPLACE2 macro */
+TYPEMAP_INPLACE2(char, PyArray_CHAR )
+TYPEMAP_INPLACE2(unsigned char, PyArray_UBYTE )
+TYPEMAP_INPLACE2(signed char, PyArray_SBYTE )
+TYPEMAP_INPLACE2(short, PyArray_SHORT )
+TYPEMAP_INPLACE2(int, PyArray_INT )
+TYPEMAP_INPLACE2(long, PyArray_LONG )
+TYPEMAP_INPLACE2(float, PyArray_FLOAT )
+TYPEMAP_INPLACE2(double, PyArray_DOUBLE)
+TYPEMAP_INPLACE2(long double, PyArray_LONGDOUBLE)
+TYPEMAP_INPLACE2(npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_INPLACE2(npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_INPLACE2(npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_INPLACE2(const char, PyArray_CHAR )
+TYPEMAP_INPLACE2(const unsigned char, PyArray_UBYTE )
+TYPEMAP_INPLACE2(const signed char, PyArray_SBYTE )
+TYPEMAP_INPLACE2(const short, PyArray_SHORT )
+TYPEMAP_INPLACE2(const int, PyArray_INT )
+TYPEMAP_INPLACE2(const long, PyArray_LONG )
+TYPEMAP_INPLACE2(const float, PyArray_FLOAT )
+TYPEMAP_INPLACE2(const double, PyArray_DOUBLE)
+TYPEMAP_INPLACE2(const long double, PyArray_LONGDOUBLE)
+TYPEMAP_INPLACE2(const npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_INPLACE2(const npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_INPLACE2(const npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_INPLACE2(PyObject, PyArray_OBJECT)
+
+#undef TYPEMAP_INPLACE2
+
+
+
+
+/* TYPEMAP_ARGOUT macros
+ *
+ * This family of typemaps allows output C arguments of the form
+ *
+ * (type* ARGOUT_ARRAY[ANY])
+ * (type* ARGOUT_ARRAY[ANY][ANY])
+ *
+ * where "type" is any type supported by the Numeric module, to be
+ * called in python with an argument list of a single contiguous
+ * Numeric array. This can be applied to an existing function using
+ * the %apply directive:
+ *
+ * %apply (double* ARGOUT_ARRAY[ANY] {double series, int length}
+ * %apply (double* ARGOUT_ARRAY[ANY][ANY]) {double* mx, int rows, int cols}
+ * void negate(double* series, int length);
+ * void normalize(double* mx, int rows, int cols);
+ *
+ *
+ * or with
+ *
+ * void sum(double* ARGOUT_ARRAY[ANY]);
+ * void sum(double* ARGOUT_ARRAY[ANY][ANY]);
+ */
+
+ /* One dimensional input/output arrays */
+%define TYPEMAP_ARGOUT1(type,typecode)
+%typemap(in,numinputs=0) type ARGOUT_ARRAY[ANY] {
+ $1 = (type*) malloc($1_dim0*sizeof(type));
+ if (!$1) {
+ PyErr_SetString(PyExc_RuntimeError, "Failed to allocate memory");
+ SWIG_fail;
+ }
+}
+%typemap(argout) ARGOUT_ARRAY[ANY] {
+ int dimensions[1] = {$1_dim0};
+ PyObject* outArray = PyArray_FromDimsAndData(1, dimensions, typecode, (char*)$1);
+}
+%enddef
+
+/* Define concrete examples of the TYPEMAP_ARGOUT1 macro */
+TYPEMAP_ARGOUT1(char, PyArray_CHAR )
+TYPEMAP_ARGOUT1(unsigned char, PyArray_UBYTE )
+TYPEMAP_ARGOUT1(signed char, PyArray_SBYTE )
+TYPEMAP_ARGOUT1(short, PyArray_SHORT )
+TYPEMAP_ARGOUT1(int, PyArray_INT )
+TYPEMAP_ARGOUT1(long, PyArray_LONG )
+TYPEMAP_ARGOUT1(float, PyArray_FLOAT )
+TYPEMAP_ARGOUT1(double, PyArray_DOUBLE)
+TYPEMAP_ARGOUT1(long double, PyArray_LONGDOUBLE)
+TYPEMAP_ARGOUT1(npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_ARGOUT1(npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_ARGOUT1(npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_ARGOUT1(PyObject, PyArray_OBJECT)
+
+#undef TYPEMAP_ARGOUT1
+
+ /* Two dimensional input/output arrays */
+%define TYPEMAP_ARGOUT2(type,typecode)
+ %typemap(in) (type* ARGOUT_ARRAY2, int DIM1, int DIM2) (PyArrayObject* temp=NULL) {
+ temp = obj_to_array_no_conversion($input,typecode);
+ if (!temp || !require_contiguous(temp)) SWIG_fail;
+ $1 = (type*) temp->data;
+ $2 = temp->dimensions[0];
+ $3 = temp->dimensions[1];
+}
+%enddef
+
+/* Define concrete examples of the TYPEMAP_ARGOUT2 macro */
+TYPEMAP_ARGOUT2(char, PyArray_CHAR )
+TYPEMAP_ARGOUT2(unsigned char, PyArray_UBYTE )
+TYPEMAP_ARGOUT2(signed char, PyArray_SBYTE )
+TYPEMAP_ARGOUT2(short, PyArray_SHORT )
+TYPEMAP_ARGOUT2(int, PyArray_INT )
+TYPEMAP_ARGOUT2(long, PyArray_LONG )
+TYPEMAP_ARGOUT2(float, PyArray_FLOAT )
+TYPEMAP_ARGOUT2(double, PyArray_DOUBLE)
+TYPEMAP_ARGOUT2(long double, PyArray_LONGDOUBLE)
+TYPEMAP_ARGOUT2(npy_cfloat, PyArray_CFLOAT )
+TYPEMAP_ARGOUT2(npy_cdouble, PyArray_CDOUBLE)
+TYPEMAP_ARGOUT2(npy_clongdouble, PyArray_CLONGDOUBLE)
+TYPEMAP_ARGOUT2(PyObject, PyArray_OBJECT)
+
+#undef TYPEMAP_ARGOUT2
Added: trunk/Lib/sparse/sparsetools/sparsetools.h
===================================================================
--- trunk/Lib/sparse/sparsetools/sparsetools.h 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparsetools/sparsetools.h 2007-01-06 06:34:41 UTC (rev 2496)
@@ -0,0 +1,863 @@
+#ifndef SPARSETOOLS_H
+#define SPARSETOOLS_H
+
+
+/*
+ * sparsetools.h
+ * A collection of CSR/CSC/COO matrix conversion and arithmetic functions.
+ *
+ * Authors:
+ * Nathan Bell
+ *
+ * Revisions:
+ * 01/06/2007 - initial inclusion into SciPy
+ *
+ */
+
+
+
+
+#include <vector>
+
+
+
+
+
+/*
+ * Return zero of the appropriate type
+ */
+template <class T>
+T ZERO(){
+ T temp = {0};
+ return temp;
+}
+
+
+
+/*
+ * Compute B = A for CSR matrix A, CSC matrix B
+ *
+ * Also, with the appropriate arguments can also be used to:
+ * - compute B = A^t for CSR matrix A, CSR matrix B
+ * - compute B = A^t for CSC matrix A, CSC matrix B
+ * - convert CSC->CSR
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ *
+ * Output Arguments:
+ * vec<int> Bp - row pointer
+ * vec<int> Bj - column indices
+ * vec<T> Bx - nonzeros
+ *
+ * Note:
+ * Output arrays Bp,Bj,Bx will be allocated within in the method
+ *
+ * Note:
+ * Input: column indices *are not* assumed to be in sorted order
+ * Output: row indices *will be* in sorted order
+ *
+ * Complexity: Linear. Specifically O(nnz(A) + max(n_row,n_col))
+ *
+ */
+template <class T>
+void csrtocsc(const int n_row,
+ const int n_col,
+ const int Ap[],
+ const int Aj[],
+ const T Ax[],
+ std::vector<int>* Bp,
+ std::vector<int>* Bi,
+ std::vector<T>* Bx)
+{
+ int NNZ = Ap[n_row];
+
+ *Bp = std::vector<int>(n_col+1);
+ *Bi = std::vector<int>(NNZ);
+ *Bx = std::vector<T>(NNZ);
+
+ std::vector<int> nnz_per_col(n_col,0); //temp array
+
+ //compute number of non-zero entries per column of A
+ for (int i = 0; i < NNZ; i++){
+ nnz_per_col[Aj[i]]++;
+ }
+
+ //cumsum the nnz_per_col to get Bp[]
+ for(int i = 0, cumsum = 0; i < n_col; i++){
+ (*Bp)[i] = cumsum;
+ cumsum += nnz_per_col[i];
+ nnz_per_col[i] = 0; //reset count
+ }
+ (*Bp)[n_col] = NNZ;
+
+ for(int i = 0; i < n_row; i++){
+ int row_start = Ap[i];
+ int row_end = Ap[i+1];
+ for(int j = row_start; j < row_end; j++){
+ int col = Aj[j];
+ int k = (*Bp)[col] + nnz_per_col[col];
+
+ (*Bi)[k] = i;
+ (*Bx)[k] = Ax[j];
+
+ nnz_per_col[col]++;
+ }
+ }
+}
+
+
+
+
+/*
+ * Compute B = A for CSR matrix A, COO matrix B
+ *
+ * Also, with the appropriate arguments can also be used to:
+ * - convert CSC->COO
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ *
+ * Output Arguments:
+ * vec<int> Bi - row indices
+ * vec<int> Bj - column indices
+ * vec<T> Bx - nonzeros
+ *
+ * Note:
+ * Output arrays Bi,Bj,Bx will be allocated within in the method
+ *
+ * Note:
+ * Complexity: Linear.
+ *
+ */
+template<class T>
+void csrtocoo(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Aj[],
+ const T Ax[],
+ std::vector<int>* Bi,
+ std::vector<int>* Bj,
+ std::vector<T>* Bx)
+{
+ for(int i = 0; i < n_row; i++){
+ int row_start = Ap[i];
+ int row_end = Ap[i+1];
+ for(int jj = row_start; jj < row_end; jj++){
+ Bi->push_back(i);
+ Bj->push_back(Aj[jj]);
+ Bx->push_back(Ax[jj]);
+ }
+ }
+}
+
+
+
+/*
+ * Compute C = A*B for CSR matrices A,B
+ *
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in B (hence C is n_row by n_col)
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ * int Bp[?] - row pointer
+ * int Bj[nnz(B)] - column indices
+ * T Bx[nnz(B)] - nonzeros
+ * Output Arguments:
+ * vec<int> Cp - row pointer
+ * vec<int> Cj - column indices
+ * vec<T> Cx - nonzeros
+ *
+ * Note:
+ * Output arrays Cp,Cj, and Cx will be allocated within in the method
+ *
+ * Note:
+ * Input: A and B column indices *are not* assumed to be in sorted order
+ * Output: C column indices *are not* assumed to be in sorted order
+ * Cx will not contain any zero entries
+ *
+ * Complexity: O(n_row*K^2 + max(n_row,n_col))
+ * where K is the maximum nnz in a row of A
+ * and column of B.
+ *
+ *
+ * This implementation closely follows the SMMP algorithm:
+ *
+ * "Sparse Matrix Multiplication Package (SMMP)"
+ * Randolph E. Bank and Craig C. Douglas
+ *
+ * http://citeseer.ist.psu.edu/445062.html
+ * http://www.mgnet.org/~douglas/ccd-codes.html
+ *
+ */
+template<class T>
+void csrmucsr(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Aj[],
+ const T Ax[],
+ const int Bp[],
+ const int Bj[],
+ const T Bx[],
+ std::vector<int>* Cp,
+ std::vector<int>* Cj,
+ std::vector<T>* Cx)
+{
+ *Cp = std::vector<int>(n_row+1,0);
+ Cj->clear();
+ Cx->clear();
+
+ const T zero = ZERO<T>();
+
+ std::vector<int> index(n_col,-1);
+ std::vector<T> sums(n_col,zero);
+
+ for(int i = 0; i < n_row; i++){
+ int istart = -1;
+ int length = 0;
+
+ for(int jj = Ap[i]; jj < Ap[i+1]; jj++){
+ int j = Aj[jj];
+ for(int kk = Bp[j]; kk < Bp[j+1]; kk++){
+ int k = Bj[kk];
+
+ sums[k] += Ax[jj]*Bx[kk];
+
+ if(index[k] == -1){
+ index[k] = istart;
+ istart = k;
+ length++;
+ }
+ }
+ }
+
+ for(int jj = 0; jj < length; jj++){
+ if(sums[istart] != zero){
+ Cj->push_back(istart);
+ Cx->push_back(sums[istart]);
+ }
+
+ int temp = istart;
+ istart = index[istart];
+
+ index[temp] = -1; //clear arrays
+ sums[temp] = zero;
+ }
+
+ (*Cp)[i+1] = Cx->size();
+ }
+}
+
+
+
+
+/*
+ * Compute C = A+B for CSR matrices A,B
+ *
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A (and B)
+ * int n_col - number of columns in A (and B)
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ * int Bp[?] - row pointer
+ * int Bj[nnz(B)] - column indices
+ * T Bx[nnz(B)] - nonzeros
+ * Output Arguments:
+ * vec<int> Cp - row pointer
+ * vec<int> Cj - column indices
+ * vec<T> Cx - nonzeros
+ *
+ * Note:
+ * Output arrays Cp,Cj, and Cx will be allocated within in the method
+ *
+ * Note:
+ * Input: A and B column indices *are not* assumed to be in sorted order
+ * Output: C column indices *are not* assumed to be in sorted order
+ * Cx will not contain any zero entries
+ *
+ */
+template <class T>
+void csrplcsr(const int n_row,
+ const int n_col,
+ const int Ap[],
+ const int Aj[],
+ const T Ax[],
+ const int Bp[],
+ const int Bj[],
+ const T Bx[],
+ std::vector<int>* Cp,
+ std::vector<int>* Cj,
+ std::vector<T> * Cx)
+{
+
+ *Cp = std::vector<int>(n_row+1,0);
+ Cj->clear();
+ Cx->clear();
+
+ const T zero = ZERO<T>();
+
+ std::vector<int> index(n_col,-1);
+ std::vector<T> sums(n_col,zero);
+
+ for(int i = 0; i < n_row; i++){
+ int istart = -1;
+ int length = 0;
+
+ //add a row of A to sums
+ for(int jj = Ap[i]; jj < Ap[i+1]; jj++){
+ int j = Aj[jj];
+ sums[j] += Ax[jj];
+
+ if(index[j] == -1){
+ index[j] = istart;
+ istart = j;
+ length++;
+ }
+ }
+
+ //add a row of B to sums
+ for(int jj = Bp[i]; jj < Bp[i+1]; jj++){
+ int j = Bj[jj];
+ sums[j] += Bx[jj];
+
+ if(index[j] == -1){
+ index[j] = istart;
+ istart = j;
+ length++;
+ }
+ }
+
+
+ for(int jj = 0; jj < length; jj++){
+ if(sums[istart] != zero){
+ Cj->push_back(istart);
+ Cx->push_back(sums[istart]);
+ }
+
+ int temp = istart;
+ istart = index[istart];
+
+ index[temp] = -1;
+ sums[temp] = zero;
+ }
+
+ (*Cp)[i+1] = Cx->size();
+ }
+}
+
+/*
+ * Compute C = A (elmul) B for CSR matrices A,B
+ *
+ * (elmul) - elementwise multiplication
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A (and B)
+ * int n_col - number of columns in A (and B)
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ * int Bp[?] - row pointer
+ * int Bj[nnz(B)] - column indices
+ * T Bx[nnz(B)] - nonzeros
+ * Output Arguments:
+ * vec<int> Cp - row pointer
+ * vec<int> Cj - column indices
+ * vec<T> Cx - nonzeros
+ *
+ * Note:
+ * Output arrays Cp,Cj, and Cx will be allocated within in the method
+ *
+ * Note:
+ * Input: A and B column indices *are not* assumed to be in sorted order
+ * Output: C column indices *are not* assumed to be in sorted order
+ * Cx will not contain any zero entries
+ *
+ */
+template <class T>
+void csrelmulcsr(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Aj[],
+ const T Ax[],
+ const int Bp[],
+ const int Bj[],
+ const T Bx[],
+ std::vector<int>* Cp,
+ std::vector<int>* Cj,
+ std::vector<T>* Cx)
+{
+ *Cp = std::vector<int>(n_row+1,0);
+ Cj->clear();
+ Cx->clear();
+
+ const T zero = ZERO<T>();
+
+ std::vector<int> index(n_col,-1);
+ std::vector<T> A_row(n_col,zero);
+ std::vector<T> B_row(n_col,zero);
+
+ for(int i = 0; i < n_row; i++){
+ int istart = -1;
+ int length = 0;
+
+ //add a row of A to A_row
+ for(int jj = Ap[i]; jj < Ap[i+1]; jj++){
+ int j = Aj[jj];
+
+ A_row[j] += Ax[jj];
+
+ if(index[j] == -1){
+ index[j] = istart;
+ istart = j;
+ length++;
+ }
+ }
+
+ //add a row of B to B_row
+ for(int jj = Bp[i]; jj < Bp[i+1]; jj++){
+ int j = Bj[jj];
+
+ B_row[j] += Bx[jj];
+
+ if(index[j] == -1){
+ index[j] = istart;
+ istart = j;
+ length++;
+ }
+ }
+
+
+ for(int jj = 0; jj < length; jj++){
+ T prod = A_row[istart] * B_row[istart];
+
+ if(prod != zero){
+ Cj->push_back(istart);
+ Cx->push_back(prod);
+ }
+
+ int temp = istart;
+ istart = index[istart];
+
+ index[temp] = -1;
+ A_row[temp] = zero;
+ B_row[temp] = zero;
+ }
+
+ (*Cp)[i+1] = Cx->size();
+ }
+}
+
+
+/*
+ * Compute B = A for COO matrix A, CSR matrix B
+ *
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int Ai[nnz(A)] - row indices
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ * Output Arguments:
+ * vec<int> Bp - row pointer
+ * vec<int> Bj - column indices
+ * vec<T> Bx - nonzeros
+ *
+ * Note:
+ * Output arrays Bp,Bj,Bx will be allocated within in the method
+ *
+ * Note:
+ * Input: row and column indices *are not* assumed to be ordered
+ * duplicate (i,j) entries will be summed together
+ *
+ * Output: CSR column indices *will be* in sorted order
+ *
+ * Complexity: Linear. Specifically O(nnz(A) + max(n_row,n_col))
+ *
+ */
+template<class T>
+void cootocsr(const int n_row,
+ const int n_col,
+ const int NNZ,
+ const int Ai[],
+ const int Aj[],
+ const T Ax[],
+ std::vector<int>* Bp,
+ std::vector<int>* Bj,
+ std::vector<T>* Bx)
+{
+ std::vector<int> tempBp(n_row+1,0);
+ std::vector<int> tempBj(NNZ);
+ std::vector<T> tempBx(NNZ);
+
+ std::vector<int> nnz_per_row(n_row,0); //temp array
+
+ //compute nnz per row, then compute Bp
+ for(int i = 0; i < NNZ; i++){
+ nnz_per_row[Ai[i]]++;
+ }
+ for(int i = 0, cumsum = 0; i < n_row; i++){
+ tempBp[i] = cumsum;
+ cumsum += nnz_per_row[i];
+ nnz_per_row[i] = 0; //reset count
+ }
+ tempBp[n_row] = NNZ;
+
+
+ //write Aj,Ax into tempBj,tempBx
+ for(int i = 0; i < NNZ; i++){
+ int row = Ai[i];
+ int n = tempBp[row] + nnz_per_row[row];
+
+ tempBj[n] = Aj[i];
+ tempBx[n] = Ax[i];
+
+ nnz_per_row[row]++;
+ }
+ //now tempBp,tempBj,tempBx form a CSR representation (with duplicates)
+
+
+ //use (tempB + 0) to sum duplicates
+ std::vector<int> Xp(n_row+1,0); //row pointer for an empty matrix
+
+ csrplcsr<T>(n_row,n_col,
+ &tempBp[0],&tempBj[0],&tempBx[0],
+ &Xp[0],NULL,NULL,
+ Bp,Bj,Bx);
+}
+
+
+
+
+
+/*
+ * Compute Y = A*X for CSR matrix A and dense vectors X,Y
+ *
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[n_col] - nonzeros
+ * T Xx[n_col] - nonzeros
+ *
+ * Output Arguments:
+ * vec<T> Yx - nonzeros (real part)
+ *
+ * Note:
+ * Output array Xx will be allocated within in the method
+ *
+ * Complexity: Linear. Specifically O(nnz(A) + max(n_row,n_col))
+ *
+ */
+template <class T>
+void csrmux(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Aj[],
+ const T Ax[],
+ const T Xx[],
+ std::vector<T>* Yx)
+{
+ const T zero = ZERO<T>();
+
+ *Yx = std::vector<T>(n_row,zero);
+
+ for(int i = 0; i < n_row; i++){
+ int row_start = Ap[i];
+ int row_end = Ap[i+1];
+
+ T& Yx_i = (*Yx)[i];
+ for(int jj = row_start; jj < row_end; jj++){
+ Yx_i += Ax[jj] * Xx[Aj[jj]];
+ }
+ }
+}
+
+
+
+/*
+ * Compute Y = A*X for CSC matrix A and dense vectors X,Y
+ *
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int Ap[n_row+1] - column pointer
+ * int Ai[nnz(A)] - row indices
+ * T Ax[n_col] - nonzeros (real part)
+ * T Xx[n_col] - nonzeros (real part)
+ * bool do_complex - switch scalar/complex modes
+ *
+ * Output Arguments:
+ * vec<T> Yx - nonzeros (real part)
+ *
+ * Note:
+ * Output arrays Xx will be allocated within in the method
+ *
+ * Complexity: Linear. Specifically O(nnz(A) + max(n_row,n_col))
+ *
+ */
+template <class T>
+void cscmux(const int n_row,
+ const int n_col,
+ const int Ap[],
+ const int Ai[],
+ const T Ax[],
+ const T Xx[],
+ std::vector<T>* Yx)
+{
+ const T zero = ZERO<T>();
+
+ *Yx = std::vector<T>(n_row,zero);
+
+ for(int j = 0; j < n_col; j++){
+ int col_start = Ap[j];
+ int col_end = Ap[j+1];
+
+ for(int ii = col_start; ii < col_end; ii++){
+ int row = Ai[ii];
+ (*Yx)[row] += Ax[ii] * Xx[j];
+ }
+ }
+}
+
+
+
+
+/*
+ * Construct CSC matrix A from diagonals
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int n_diags - number of diagonals
+ * int diags_indx[n_diags] - where to place each diagonal
+ * T diags[n_diags][min(n_row,n_col)] - diagonals
+ *
+ * Output Arguments:
+ * vec<int> Ap - row pointer
+ * vec<int> Aj - column indices
+ * vec<T> Ax - nonzeros
+ *
+ * Note:
+ * Output arrays Ap,Aj,Ax will be allocated within in the method
+ *
+ * Note:
+ * Output: row indices are not in sorted order
+ *
+ * Complexity: Linear
+ *
+ */
+template<class T>
+void spdiags(const int n_row,
+ const int n_col,
+ const int n_diag,
+ const int offsets[],
+ const T diags[],
+ std::vector<int> * Ap,
+ std::vector<int> * Ai,
+ std::vector<T> * Ax)
+{
+ const int diags_length = std::min(n_row,n_col);
+ Ap->push_back(0);
+
+ for(int i = 0; i < n_col; i++){
+ for(int j = 0; j < n_diag; j++){
+ if(offsets[j] <= 0){ //sub-diagonal
+ int row = i - offsets[j];
+ if (row >= n_row){ continue; }
+
+ Ai->push_back(row);
+ Ax->push_back(diags[j*diags_length + i]);
+
+ } else { //super-diagonal
+ int row = i - offsets[j];
+ if (row < 0 || row >= n_row){ continue; }
+
+ Ai->push_back(row);
+ Ax->push_back(diags[j*diags_length + row]);
+ }
+ }
+ Ap->push_back(Ai->size());
+ }
+}
+
+
+
+/*
+ * Compute M = A for CSR matrix A, dense matrix M
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ * T Mx[n_row*n_col] - dense matrix
+ *
+ * Note:
+ * Output arrays Mx are assumed to be allocated and
+ * initialized to 0 by the caller.
+ *
+ */
+template<class T>
+void csrtodense(const int n_row,
+ const int n_col,
+ const int Ap[],
+ const int Aj[],
+ const T Ax[],
+ T Mx[])
+{
+ int row_base = 0;
+ for(int i = 0; i < n_row; i++){
+ int row_start = Ap[i];
+ int row_end = Ap[i+1];
+ for(int jj = row_start; jj < row_end; jj++){
+ int j = Aj[jj];
+
+ Mx[row_base + j] = Ax[jj];
+ }
+ row_base += n_col;
+ }
+}
+
+
+
+/*
+ * Compute A = M for CSR matrix A, dense matrix M
+ *
+ * Input Arguments:
+ * int n_row - number of rows in A
+ * int n_col - number of columns in A
+ * T Mx[n_row*n_col] - dense matrix
+ * int Ap[n_row+1] - row pointer
+ * int Aj[nnz(A)] - column indices
+ * T Ax[nnz(A)] - nonzeros
+ *
+ * Note:
+ * Output arrays Ap,Aj,Ax will be allocated within the method
+ *
+ */
+template<class T>
+void densetocsr(const int n_row,
+ const int n_col,
+ const T Mx[],
+ std::vector<int>* Ap,
+ std::vector<int>* Aj,
+ std::vector<T>* Ax)
+{
+ const T zero = ZERO<T>();
+ const T* x_ptr = Mx;
+
+ Ap->push_back(0);
+ for(int i = 0; i < n_row; i++){
+ for(int j = 0; j < n_col; j++){
+ if(*x_ptr != zero){
+ Aj->push_back(j);
+ Ax->push_back(*x_ptr);
+ }
+ x_ptr++;
+ }
+ Ap->push_back(Aj->size());
+ }
+}
+
+
+
+/*
+ * Derived methods
+ */
+template <class T>
+void csctocsr(const int n_row,
+ const int n_col,
+ const int Ap[],
+ const int Ai[],
+ const T Ax[],
+ std::vector<int>* Bp,
+ std::vector<int>* Bj,
+ std::vector<T>* Bx)
+{ csrtocsc<T>(n_col,n_row,Ap,Ai,Ax,Bp,Bj,Bx); }
+
+template<class T>
+void csctocoo(const int n_row,
+ const int n_col,
+ const int Ap[],
+ const int Ai[],
+ const T Ax[],
+ std::vector<int>* Bi,
+ std::vector<int>* Bj,
+ std::vector<T>* Bx)
+{ csrtocoo<T>(n_col,n_row,Ap,Ai,Ax,Bj,Bi,Bx); }
+
+template<class T>
+void cscmucsc(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Ai[],
+ const T Ax[],
+ const int Bp[],
+ const int Bi[],
+ const T Bx[],
+ std::vector<int>* Cp,
+ std::vector<int>* Ci,
+ std::vector<T> * Cx)
+{ csrmucsr<T>(n_col,n_row,Bp,Bi,Bx,Ap,Ai,Ax,Cp,Ci,Cx); }
+
+template <class T>
+void cscplcsc(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Ai[],
+ const T Ax[],
+ const int Bp[],
+ const int Bi[],
+ const T Bx[],
+ std::vector<int>* Cp,
+ std::vector<int>* Ci,
+ std::vector<T>* Cx)
+{ csrplcsr<T>(n_col,n_row,Ap,Ai,Ax,Bp,Bi,Bx,Cp,Ci,Cx); }
+
+template <class T>
+void cscelmulcsc(const int n_row,
+ const int n_col,
+ const int Ap [],
+ const int Ai[],
+ const T Ax[],
+ const int Bp[],
+ const int Bi[],
+ const T Bx[],
+ std::vector<int>* Cp,
+ std::vector<int>* Ci,
+ std::vector<T>* Cx)
+{ csrelmulcsr<T>(n_col,n_row,Ap,Ai,Ax,Bp,Bi,Bx,Cp,Ci,Cx); }
+
+template<class T>
+void cootocsc(const int n_row,
+ const int n_col,
+ const int NNZ,
+ const int Ai[],
+ const int Aj[],
+ const T Ax[],
+ std::vector<int>* Bp,
+ std::vector<int>* Bi,
+ std::vector<T>* Bx)
+{ cootocsr<T>(n_col,n_row,NNZ,Aj,Ai,Ax,Bp,Bi,Bx); }
+
+
+
+#endif
Added: trunk/Lib/sparse/sparsetools/sparsetools.i
===================================================================
--- trunk/Lib/sparse/sparsetools/sparsetools.i 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparsetools/sparsetools.i 2007-01-06 06:34:41 UTC (rev 2496)
@@ -0,0 +1,440 @@
+/* -*- C -*- (not really, but good for syntax highlighting) */
+%module sparsetools
+
+ /* why does SWIG complain about int arrays? a typecheck is provided */
+#pragma SWIG nowarn=467
+
+%{
+#define SWIG_FILE_WITH_INIT
+#include "numpy/arrayobject.h"
+#include "complex_ops.h"
+#include "sparsetools.h"
+%}
+
+%feature("autodoc", "1");
+
+%include "numpy.i"
+
+%init %{
+ import_array();
+%}
+
+
+%{
+/*!
+ Appends @a what to @a where. On input, @a where need not to be a tuple, but on
+ return it always is.
+
+ @par Revision history:
+ - 17.02.2005, c
+*/
+PyObject *helper_appendToTuple( PyObject *where, PyObject *what ) {
+ PyObject *o2, *o3;
+
+ if ((!where) || (where == Py_None)) {
+ where = what;
+ } else {
+ if (!PyTuple_Check( where )) {
+ o2 = where;
+ where = PyTuple_New( 1 );
+ PyTuple_SetItem( where, 0, o2 );
+ }
+ o3 = PyTuple_New( 1 );
+ PyTuple_SetItem( o3, 0, what );
+ o2 = where;
+ where = PySequence_Concat( o2, o3 );
+ Py_DECREF( o2 );
+ Py_DECREF( o3 );
+ }
+ return where;
+}
+
+
+
+
+%}
+
+
+
+/*
+ * Use STL vectors for ARGOUTs
+ */
+%define VEC_ARRAY_ARGOUT( ctype, atype )
+%typemap( in, numinputs=0 ) std::vector<ctype>* array_argout( std::vector<ctype>* tmp ) {
+ tmp = new std::vector<ctype>();
+ $1 = tmp;
+};
+%typemap( argout ) std::vector<ctype>* array_argout {
+ int length = ($1)->size();
+ PyObject *obj = PyArray_FromDims(1, &length, PyArray_##atype);
+ memcpy(PyArray_DATA(obj),&((*($1))[0]),sizeof(ctype)*length);
+ delete $1;
+ $result = helper_appendToTuple( $result, (PyObject *)obj );
+};
+%enddef
+
+
+
+
+%include "typemaps.i"
+
+%typemap(typecheck) int *, const int *, int [], const int []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_INT)) ? 1 : 0;
+};
+
+%typemap(typecheck) float *, const float *, float [], const float []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_FLOAT)) ? 1 : 0;
+};
+
+%typemap(typecheck) double *, const double *, double [], const double []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_DOUBLE)) ? 1 : 0;
+};
+
+%typemap(typecheck) long double *, const long double *, long double [], const long double []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_LONGDOUBLE)) ? 1 : 0;
+};
+
+%typemap(typecheck) npy_cfloat *, const npy_cfloat *, npy_cfloat [], const npy_cfloat []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_CFLOAT)) ? 1 : 0;
+};
+
+%typemap(typecheck) npy_cdouble *, const npy_cdouble *, npy_cdouble [], const npy_cdouble []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_CDOUBLE)) ? 1 : 0;
+};
+
+%typemap(typecheck) npy_clongdouble *, const npy_clongdouble *, npy_clongdouble [], const npy_clongdouble []
+{
+ $1 = (is_array($input) && PyArray_CanCastSafely(PyArray_TYPE($input),PyArray_CLONGDOUBLE)) ? 1 : 0;
+};
+
+
+
+
+
+
+/*
+ * IN types
+ */
+%apply int * IN_ARRAY1 {
+ const int Ap [ ],
+ const int Ai [ ],
+ const int Aj [ ],
+ const int Bp [ ],
+ const int Bi [ ],
+ const int Bj [ ],
+ const int offsets [ ]
+};
+
+%apply float * IN_ARRAY1 {
+ const float Ax [ ],
+ const float Bx [ ],
+ const float Xx [ ],
+ const float Yx [ ]
+};
+
+%apply double * IN_ARRAY1 {
+ const double Ax [ ],
+ const double Bx [ ],
+ const double Xx [ ],
+ const double Yx [ ]
+};
+
+%apply long double * IN_ARRAY1 {
+ const long double Ax [ ],
+ const long double Bx [ ],
+ const long double Xx [ ],
+ const long double Yx [ ]
+};
+
+%apply npy_cfloat * IN_ARRAY1 {
+ const npy_cfloat Ax [ ],
+ const npy_cfloat Bx [ ],
+ const npy_cfloat Xx [ ],
+ const npy_cfloat Yx [ ]
+};
+
+%apply npy_cdouble * IN_ARRAY1 {
+ const npy_cdouble Ax [ ],
+ const npy_cdouble Bx [ ],
+ const npy_cdouble Xx [ ],
+ const npy_cdouble Yx [ ]
+};
+
+%apply npy_clongdouble * IN_ARRAY1 {
+ const npy_clongdouble Ax [ ],
+ const npy_clongdouble Bx [ ],
+ const npy_clongdouble Xx [ ],
+ const npy_clongdouble Yx [ ]
+};
+
+
+%apply float * IN_ARRAY2 { const float Mx[] }
+%apply double * IN_ARRAY2 { const double Mx[] }
+%apply long double * IN_ARRAY2 { const long double Mx[] }
+%apply npy_cfloat * IN_ARRAY2 { const npy_cfloat Mx[] }
+%apply npy_cdouble * IN_ARRAY2 { const npy_cdouble Mx[] }
+%apply npy_clongdouble * IN_ARRAY2 { const npy_longdouble Mx[] }
+
+
+%apply float * IN_ARRAY2 { const float diags[] }
+%apply double * IN_ARRAY2 { const double diags[] }
+%apply long double * IN_ARRAY2 { const long double diags[] }
+%apply npy_cfloat * IN_ARRAY2 { const npy_cfloat diags[] }
+%apply npy_cdouble * IN_ARRAY2 { const npy_cdouble diags[] }
+%apply npy_clongdouble * IN_ARRAY2 { const npy_longdouble diags[] }
+
+
+
+/*
+ * OUT types
+ */
+VEC_ARRAY_ARGOUT( int, INT )
+%apply std::vector<int>* array_argout {
+ std::vector<int>* Ap,
+ std::vector<int>* Ai,
+ std::vector<int>* Aj,
+ std::vector<int>* Bp,
+ std::vector<int>* Bi,
+ std::vector<int>* Bj,
+ std::vector<int>* Cp,
+ std::vector<int>* Ci,
+ std::vector<int>* Cj
+};
+
+VEC_ARRAY_ARGOUT( float, FLOAT )
+%apply std::vector<float>* array_argout {
+ std::vector<float>* Ax,
+ std::vector<float>* Bx,
+ std::vector<float>* Cx,
+ std::vector<float>* Xx,
+ std::vector<float>* Yx
+};
+
+VEC_ARRAY_ARGOUT( double, DOUBLE )
+%apply std::vector<double>* array_argout {
+ std::vector<double>* Ax,
+ std::vector<double>* Bx,
+ std::vector<double>* Cx,
+ std::vector<double>* Xx,
+ std::vector<double>* Yx
+};
+
+
+VEC_ARRAY_ARGOUT( long double, LONGDOUBLE )
+%apply std::vector<long double>* array_argout {
+ std::vector<long double>* Ax,
+ std::vector<long double>* Bx,
+ std::vector<long double>* Cx,
+ std::vector<long double>* Xx,
+ std::vector<long double>* Yx
+};
+
+VEC_ARRAY_ARGOUT( npy_cfloat, CFLOAT )
+%apply std::vector<npy_cfloat>* array_argout {
+ std::vector<npy_cfloat>* Ax,
+ std::vector<npy_cfloat>* Bx,
+ std::vector<npy_cfloat>* Cx,
+ std::vector<npy_cfloat>* Xx,
+ std::vector<npy_cfloat>* Yx
+};
+
+
+VEC_ARRAY_ARGOUT( npy_cdouble, CDOUBLE )
+%apply std::vector<npy_cdouble>* array_argout {
+ std::vector<npy_cdouble>* Ax,
+ std::vector<npy_cdouble>* Bx,
+ std::vector<npy_cdouble>* Cx,
+ std::vector<npy_cdouble>* Xx,
+ std::vector<npy_cdouble>* Yx
+};
+
+
+VEC_ARRAY_ARGOUT( npy_clongdouble, CLONGDOUBLE )
+%apply std::vector<npy_clongdouble>* array_argout {
+ std::vector<npy_clongdouble>* Ax,
+ std::vector<npy_clongdouble>* Bx,
+ std::vector<npy_clongdouble>* Cx,
+ std::vector<npy_clongdouble>* Xx,
+ std::vector<npy_clongdouble>* Yx
+};
+
+
+
+/*
+ * INOUT types
+ */
+%apply float * INPLACE_ARRAY2 { float Mx [] }
+%apply double * INPLACE_ARRAY2 { double Mx [] }
+%apply long double * INPLACE_ARRAY2 { long double Mx[] }
+%apply npy_cfloat * INPLACE_ARRAY2 { npy_cfloat Mx[] }
+%apply npy_cdouble * INPLACE_ARRAY2 { npy_cdouble Mx[] }
+%apply npy_clongdouble * INPLACE_ARRAY2 { npy_longdouble Mx[] }
+
+
+
+
+
+
+%include "sparsetools.h"
+
+
+
+ /*
+ * Order may be important here, list float before double
+ */
+
+/*
+ * CSR->CSC or CSC->CSR or CSR = CSR^T or CSC = CSC^T
+ */
+%template(csrtocsc) csrtocsc<float>;
+%template(csrtocsc) csrtocsc<double>;
+%template(csrtocsc) csrtocsc<long double>;
+%template(csrtocsc) csrtocsc<npy_cfloat>;
+%template(csrtocsc) csrtocsc<npy_cdouble>;
+%template(csrtocsc) csrtocsc<npy_clongdouble>;
+
+%template(csctocsr) csctocsr<float>;
+%template(csctocsr) csctocsr<double>;
+%template(csctocsr) csctocsr<long double>;
+%template(csctocsr) csctocsr<npy_cfloat>;
+%template(csctocsr) csctocsr<npy_cdouble>;
+%template(csctocsr) csctocsr<npy_clongdouble>;
+
+
+/*
+ * CSR<->COO and CSC<->COO
+ */
+%template(csrtocoo) csrtocoo<float>;
+%template(csrtocoo) csrtocoo<double>;
+%template(csrtocoo) csrtocoo<long double>;
+%template(csrtocoo) csrtocoo<npy_cfloat>;
+%template(csrtocoo) csrtocoo<npy_cdouble>;
+%template(csrtocoo) csrtocoo<npy_clongdouble>;
+
+%template(cootocsr) cootocsr<float>;
+%template(cootocsr) cootocsr<double>;
+%template(cootocsr) cootocsr<long double>;
+%template(cootocsr) cootocsr<npy_cfloat>;
+%template(cootocsr) cootocsr<npy_cdouble>;
+%template(cootocsr) cootocsr<npy_clongdouble>;
+
+%template(csctocoo) csctocoo<float>;
+%template(csctocoo) csctocoo<double>;
+%template(csctocoo) csctocoo<long double>;
+%template(csctocoo) csctocoo<npy_cfloat>;
+%template(csctocoo) csctocoo<npy_cdouble>;
+%template(csctocoo) csctocoo<npy_clongdouble>;
+
+%template(cootocsc) cootocsc<float>;
+%template(cootocsc) cootocsc<double>;
+%template(cootocsc) cootocsc<long double>;
+%template(cootocsc) cootocsc<npy_cfloat>;
+%template(cootocsc) cootocsc<npy_cdouble>;
+%template(cootocsc) cootocsc<npy_clongdouble>;
+
+
+/*
+ * CSR+CSR and CSC+CSC
+ */
+%template(csrplcsr) csrplcsr<float>;
+%template(csrplcsr) csrplcsr<double>;
+%template(csrplcsr) csrplcsr<long double>;
+%template(csrplcsr) csrplcsr<npy_cfloat>;
+%template(csrplcsr) csrplcsr<npy_cdouble>;
+%template(csrplcsr) csrplcsr<npy_clongdouble>;
+
+%template(cscplcsc) cscplcsc<float>;
+%template(cscplcsc) cscplcsc<double>;
+%template(cscplcsc) cscplcsc<long double>;
+%template(cscplcsc) cscplcsc<npy_cfloat>;
+%template(cscplcsc) cscplcsc<npy_cdouble>;
+%template(cscplcsc) cscplcsc<npy_clongdouble>;
+
+
+
+/*
+ * CSR*CSR and CSC*CSC
+ */
+%template(csrmucsr) csrmucsr<float>;
+%template(csrmucsr) csrmucsr<double>;
+%template(csrmucsr) csrmucsr<long double>;
+%template(csrmucsr) csrmucsr<npy_cfloat>;
+%template(csrmucsr) csrmucsr<npy_cdouble>;
+%template(csrmucsr) csrmucsr<npy_clongdouble>;
+
+%template(cscmucsc) cscmucsc<float>;
+%template(cscmucsc) cscmucsc<double>;
+%template(cscmucsc) cscmucsc<long double>;
+%template(cscmucsc) cscmucsc<npy_cfloat>;
+%template(cscmucsc) cscmucsc<npy_cdouble>;
+%template(cscmucsc) cscmucsc<npy_clongdouble>;
+
+/*
+ * CSR*x and CSC*x
+ */
+%template(csrmux) csrmux<float>;
+%template(csrmux) csrmux<double>;
+%template(csrmux) csrmux<long double>;
+%template(csrmux) csrmux<npy_cfloat>;
+%template(csrmux) csrmux<npy_cdouble>;
+%template(csrmux) csrmux<npy_clongdouble>;
+
+%template(cscmux) cscmux<float>;
+%template(cscmux) cscmux<double>;
+%template(cscmux) cscmux<long double>;
+%template(cscmux) cscmux<npy_cfloat>;
+%template(cscmux) cscmux<npy_cdouble>;
+%template(cscmux) cscmux<npy_clongdouble>;
+
+
+/*
+ * CSR(elmul)CSR and CSC(elmul)CSC
+ */
+%template(csrelmulcsr) csrelmulcsr<float>;
+%template(csrelmulcsr) csrelmulcsr<double>;
+%template(csrelmulcsr) csrelmulcsr<long double>;
+%template(csrelmulcsr) csrelmulcsr<npy_cfloat>;
+%template(csrelmulcsr) csrelmulcsr<npy_cdouble>;
+%template(csrelmulcsr) csrelmulcsr<npy_clongdouble>;
+
+%template(cscelmulcsc) cscelmulcsc<float>;
+%template(cscelmulcsc) cscelmulcsc<double>;
+%template(cscelmulcsc) cscelmulcsc<long double>;
+%template(cscelmulcsc) cscelmulcsc<npy_cfloat>;
+%template(cscelmulcsc) cscelmulcsc<npy_cdouble>;
+%template(cscelmulcsc) cscelmulcsc<npy_clongdouble>;
+
+
+/*
+ * spdiags->CSC
+ */
+%template(spdiags) spdiags<float>;
+%template(spdiags) spdiags<double>;
+%template(spdiags) spdiags<long double>;
+%template(spdiags) spdiags<npy_cfloat>;
+%template(spdiags) spdiags<npy_cdouble>;
+%template(spdiags) spdiags<npy_clongdouble>;
+
+
+/*
+ * CSR<->Dense
+ */
+%template(csrtodense) csrtodense<float>;
+%template(csrtodense) csrtodense<double>;
+%template(csrtodense) csrtodense<long double>;
+%template(csrtodense) csrtodense<npy_cfloat>;
+%template(csrtodense) csrtodense<npy_cdouble>;
+%template(csrtodense) csrtodense<npy_clongdouble>;
+
+%template(densetocsr) densetocsr<float>;
+%template(densetocsr) densetocsr<double>;
+%template(densetocsr) densetocsr<long double>;
+%template(densetocsr) densetocsr<npy_cfloat>;
+%template(densetocsr) densetocsr<npy_cdouble>;
+%template(densetocsr) densetocsr<npy_clongdouble>;
Added: trunk/Lib/sparse/sparsetools/sparsetools.py
===================================================================
--- trunk/Lib/sparse/sparsetools/sparsetools.py 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparsetools/sparsetools.py 2007-01-06 06:34:41 UTC (rev 2496)
@@ -0,0 +1,398 @@
+# This file was created automatically by SWIG 1.3.28.
+# Don't modify this file, modify the SWIG interface instead.
+# This file is compatible with both classic and new-style classes.
+
+import _sparsetools
+import new
+new_instancemethod = new.instancemethod
+def _swig_setattr_nondynamic(self,class_type,name,value,static=1):
+ if (name == "thisown"): return self.this.own(value)
+ if (name == "this"):
+ if type(value).__name__ == 'PySwigObject':
+ self.__dict__[name] = value
+ return
+ method = class_type.__swig_setmethods__.get(name,None)
+ if method: return method(self,value)
+ if (not static) or hasattr(self,name):
+ self.__dict__[name] = value
+ else:
+ raise AttributeError("You cannot add attributes to %s" % self)
+
+def _swig_setattr(self,class_type,name,value):
+ return _swig_setattr_nondynamic(self,class_type,name,value,0)
+
+def _swig_getattr(self,class_type,name):
+ if (name == "thisown"): return self.this.own()
+ method = class_type.__swig_getmethods__.get(name,None)
+ if method: return method(self)
+ raise AttributeError,name
+
+import types
+try:
+ _object = types.ObjectType
+ _newclass = 1
+except AttributeError:
+ class _object : pass
+ _newclass = 0
+del types
+
+
+
+
+def csrtocsc(*args):
+ """
+ csrtocsc(int n_row, int n_col, int Ap, int Aj, float Ax, std::vector<(int)> Bp,
+ std::vector<(int)> Bi, std::vector<(float)> Bx)
+ csrtocsc(int n_row, int n_col, int Ap, int Aj, double Ax, std::vector<(int)> Bp,
+ std::vector<(int)> Bi, std::vector<(double)> Bx)
+ csrtocsc(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(long double)> Bx)
+ csrtocsc(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(npy_cfloat)> Bx)
+ csrtocsc(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(npy_cdouble)> Bx)
+ csrtocsc(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(npy_clongdouble)> Bx)
+ """
+ return _sparsetools.csrtocsc(*args)
+
+def csctocsr(*args):
+ """
+ csctocsr(int n_row, int n_col, int Ap, int Ai, float Ax, std::vector<(int)> Bp,
+ std::vector<(int)> Bj, std::vector<(float)> Bx)
+ csctocsr(int n_row, int n_col, int Ap, int Ai, double Ax, std::vector<(int)> Bp,
+ std::vector<(int)> Bj, std::vector<(double)> Bx)
+ csctocsr(int n_row, int n_col, int Ap, int Ai, long double Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(long double)> Bx)
+ csctocsr(int n_row, int n_col, int Ap, int Ai, npy_cfloat Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(npy_cfloat)> Bx)
+ csctocsr(int n_row, int n_col, int Ap, int Ai, npy_cdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(npy_cdouble)> Bx)
+ csctocsr(int n_row, int n_col, int Ap, int Ai, npy_clongdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(npy_clongdouble)> Bx)
+ """
+ return _sparsetools.csctocsr(*args)
+
+def csrtocoo(*args):
+ """
+ csrtocoo(int n_row, int n_col, int Ap, int Aj, float Ax, std::vector<(int)> Bi,
+ std::vector<(int)> Bj, std::vector<(float)> Bx)
+ csrtocoo(int n_row, int n_col, int Ap, int Aj, double Ax, std::vector<(int)> Bi,
+ std::vector<(int)> Bj, std::vector<(double)> Bx)
+ csrtocoo(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(long double)> Bx)
+ csrtocoo(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(npy_cfloat)> Bx)
+ csrtocoo(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(npy_cdouble)> Bx)
+ csrtocoo(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(npy_clongdouble)> Bx)
+ """
+ return _sparsetools.csrtocoo(*args)
+
+def cootocsr(*args):
+ """
+ cootocsr(int n_row, int n_col, int NNZ, int Ai, int Aj, float Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(float)> Bx)
+ cootocsr(int n_row, int n_col, int NNZ, int Ai, int Aj, double Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(double)> Bx)
+ cootocsr(int n_row, int n_col, int NNZ, int Ai, int Aj, long double Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(long double)> Bx)
+ cootocsr(int n_row, int n_col, int NNZ, int Ai, int Aj, npy_cfloat Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(npy_cfloat)> Bx)
+ cootocsr(int n_row, int n_col, int NNZ, int Ai, int Aj, npy_cdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(npy_cdouble)> Bx)
+ cootocsr(int n_row, int n_col, int NNZ, int Ai, int Aj, npy_clongdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bj,
+ std::vector<(npy_clongdouble)> Bx)
+ """
+ return _sparsetools.cootocsr(*args)
+
+def csctocoo(*args):
+ """
+ csctocoo(int n_row, int n_col, int Ap, int Ai, float Ax, std::vector<(int)> Bi,
+ std::vector<(int)> Bj, std::vector<(float)> Bx)
+ csctocoo(int n_row, int n_col, int Ap, int Ai, double Ax, std::vector<(int)> Bi,
+ std::vector<(int)> Bj, std::vector<(double)> Bx)
+ csctocoo(int n_row, int n_col, int Ap, int Ai, long double Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(long double)> Bx)
+ csctocoo(int n_row, int n_col, int Ap, int Ai, npy_cfloat Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(npy_cfloat)> Bx)
+ csctocoo(int n_row, int n_col, int Ap, int Ai, npy_cdouble Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(npy_cdouble)> Bx)
+ csctocoo(int n_row, int n_col, int Ap, int Ai, npy_clongdouble Ax,
+ std::vector<(int)> Bi, std::vector<(int)> Bj,
+ std::vector<(npy_clongdouble)> Bx)
+ """
+ return _sparsetools.csctocoo(*args)
+
+def cootocsc(*args):
+ """
+ cootocsc(int n_row, int n_col, int NNZ, int Ai, int Aj, float Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(float)> Bx)
+ cootocsc(int n_row, int n_col, int NNZ, int Ai, int Aj, double Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(double)> Bx)
+ cootocsc(int n_row, int n_col, int NNZ, int Ai, int Aj, long double Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(long double)> Bx)
+ cootocsc(int n_row, int n_col, int NNZ, int Ai, int Aj, npy_cfloat Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(npy_cfloat)> Bx)
+ cootocsc(int n_row, int n_col, int NNZ, int Ai, int Aj, npy_cdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(npy_cdouble)> Bx)
+ cootocsc(int n_row, int n_col, int NNZ, int Ai, int Aj, npy_clongdouble Ax,
+ std::vector<(int)> Bp, std::vector<(int)> Bi,
+ std::vector<(npy_clongdouble)> Bx)
+ """
+ return _sparsetools.cootocsc(*args)
+
+def csrplcsr(*args):
+ """
+ csrplcsr(int n_row, int n_col, int Ap, int Aj, float Ax, int Bp,
+ int Bj, float Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(float)> Cx)
+ csrplcsr(int n_row, int n_col, int Ap, int Aj, double Ax, int Bp,
+ int Bj, double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(double)> Cx)
+ csrplcsr(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ int Bp, int Bj, long double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(long double)> Cx)
+ csrplcsr(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ int Bp, int Bj, npy_cfloat Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_cfloat)> Cx)
+ csrplcsr(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ int Bp, int Bj, npy_cdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_cdouble)> Cx)
+ csrplcsr(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ int Bp, int Bj, npy_clongdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_clongdouble)> Cx)
+ """
+ return _sparsetools.csrplcsr(*args)
+
+def cscplcsc(*args):
+ """
+ cscplcsc(int n_row, int n_col, int Ap, int Ai, float Ax, int Bp,
+ int Bi, float Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(float)> Cx)
+ cscplcsc(int n_row, int n_col, int Ap, int Ai, double Ax, int Bp,
+ int Bi, double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(double)> Cx)
+ cscplcsc(int n_row, int n_col, int Ap, int Ai, long double Ax,
+ int Bp, int Bi, long double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(long double)> Cx)
+ cscplcsc(int n_row, int n_col, int Ap, int Ai, npy_cfloat Ax,
+ int Bp, int Bi, npy_cfloat Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_cfloat)> Cx)
+ cscplcsc(int n_row, int n_col, int Ap, int Ai, npy_cdouble Ax,
+ int Bp, int Bi, npy_cdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_cdouble)> Cx)
+ cscplcsc(int n_row, int n_col, int Ap, int Ai, npy_clongdouble Ax,
+ int Bp, int Bi, npy_clongdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_clongdouble)> Cx)
+ """
+ return _sparsetools.cscplcsc(*args)
+
+def csrmucsr(*args):
+ """
+ csrmucsr(int n_row, int n_col, int Ap, int Aj, float Ax, int Bp,
+ int Bj, float Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(float)> Cx)
+ csrmucsr(int n_row, int n_col, int Ap, int Aj, double Ax, int Bp,
+ int Bj, double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(double)> Cx)
+ csrmucsr(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ int Bp, int Bj, long double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(long double)> Cx)
+ csrmucsr(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ int Bp, int Bj, npy_cfloat Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_cfloat)> Cx)
+ csrmucsr(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ int Bp, int Bj, npy_cdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_cdouble)> Cx)
+ csrmucsr(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ int Bp, int Bj, npy_clongdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_clongdouble)> Cx)
+ """
+ return _sparsetools.csrmucsr(*args)
+
+def cscmucsc(*args):
+ """
+ cscmucsc(int n_row, int n_col, int Ap, int Ai, float Ax, int Bp,
+ int Bi, float Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(float)> Cx)
+ cscmucsc(int n_row, int n_col, int Ap, int Ai, double Ax, int Bp,
+ int Bi, double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(double)> Cx)
+ cscmucsc(int n_row, int n_col, int Ap, int Ai, long double Ax,
+ int Bp, int Bi, long double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(long double)> Cx)
+ cscmucsc(int n_row, int n_col, int Ap, int Ai, npy_cfloat Ax,
+ int Bp, int Bi, npy_cfloat Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_cfloat)> Cx)
+ cscmucsc(int n_row, int n_col, int Ap, int Ai, npy_cdouble Ax,
+ int Bp, int Bi, npy_cdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_cdouble)> Cx)
+ cscmucsc(int n_row, int n_col, int Ap, int Ai, npy_clongdouble Ax,
+ int Bp, int Bi, npy_clongdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_clongdouble)> Cx)
+ """
+ return _sparsetools.cscmucsc(*args)
+
+def csrmux(*args):
+ """
+ csrmux(int n_row, int n_col, int Ap, int Aj, float Ax, float Xx,
+ std::vector<(float)> Yx)
+ csrmux(int n_row, int n_col, int Ap, int Aj, double Ax, double Xx,
+ std::vector<(double)> Yx)
+ csrmux(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ long double Xx, std::vector<(long double)> Yx)
+ csrmux(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ npy_cfloat Xx, std::vector<(npy_cfloat)> Yx)
+ csrmux(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ npy_cdouble Xx, std::vector<(npy_cdouble)> Yx)
+ csrmux(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ npy_clongdouble Xx, std::vector<(npy_clongdouble)> Yx)
+ """
+ return _sparsetools.csrmux(*args)
+
+def cscmux(*args):
+ """
+ cscmux(int n_row, int n_col, int Ap, int Ai, float Ax, float Xx,
+ std::vector<(float)> Yx)
+ cscmux(int n_row, int n_col, int Ap, int Ai, double Ax, double Xx,
+ std::vector<(double)> Yx)
+ cscmux(int n_row, int n_col, int Ap, int Ai, long double Ax,
+ long double Xx, std::vector<(long double)> Yx)
+ cscmux(int n_row, int n_col, int Ap, int Ai, npy_cfloat Ax,
+ npy_cfloat Xx, std::vector<(npy_cfloat)> Yx)
+ cscmux(int n_row, int n_col, int Ap, int Ai, npy_cdouble Ax,
+ npy_cdouble Xx, std::vector<(npy_cdouble)> Yx)
+ cscmux(int n_row, int n_col, int Ap, int Ai, npy_clongdouble Ax,
+ npy_clongdouble Xx, std::vector<(npy_clongdouble)> Yx)
+ """
+ return _sparsetools.cscmux(*args)
+
+def csrelmulcsr(*args):
+ """
+ csrelmulcsr(int n_row, int n_col, int Ap, int Aj, float Ax, int Bp,
+ int Bj, float Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(float)> Cx)
+ csrelmulcsr(int n_row, int n_col, int Ap, int Aj, double Ax, int Bp,
+ int Bj, double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(double)> Cx)
+ csrelmulcsr(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ int Bp, int Bj, long double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(long double)> Cx)
+ csrelmulcsr(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ int Bp, int Bj, npy_cfloat Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_cfloat)> Cx)
+ csrelmulcsr(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ int Bp, int Bj, npy_cdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_cdouble)> Cx)
+ csrelmulcsr(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ int Bp, int Bj, npy_clongdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Cj, std::vector<(npy_clongdouble)> Cx)
+ """
+ return _sparsetools.csrelmulcsr(*args)
+
+def cscelmulcsc(*args):
+ """
+ cscelmulcsc(int n_row, int n_col, int Ap, int Ai, float Ax, int Bp,
+ int Bi, float Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(float)> Cx)
+ cscelmulcsc(int n_row, int n_col, int Ap, int Ai, double Ax, int Bp,
+ int Bi, double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(double)> Cx)
+ cscelmulcsc(int n_row, int n_col, int Ap, int Ai, long double Ax,
+ int Bp, int Bi, long double Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(long double)> Cx)
+ cscelmulcsc(int n_row, int n_col, int Ap, int Ai, npy_cfloat Ax,
+ int Bp, int Bi, npy_cfloat Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_cfloat)> Cx)
+ cscelmulcsc(int n_row, int n_col, int Ap, int Ai, npy_cdouble Ax,
+ int Bp, int Bi, npy_cdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_cdouble)> Cx)
+ cscelmulcsc(int n_row, int n_col, int Ap, int Ai, npy_clongdouble Ax,
+ int Bp, int Bi, npy_clongdouble Bx, std::vector<(int)> Cp,
+ std::vector<(int)> Ci, std::vector<(npy_clongdouble)> Cx)
+ """
+ return _sparsetools.cscelmulcsc(*args)
+
+def spdiags(*args):
+ """
+ spdiags(int n_row, int n_col, int n_diag, int offsets, float diags,
+ std::vector<(int)> Ap, std::vector<(int)> Ai,
+ std::vector<(float)> Ax)
+ spdiags(int n_row, int n_col, int n_diag, int offsets, double diags,
+ std::vector<(int)> Ap, std::vector<(int)> Ai,
+ std::vector<(double)> Ax)
+ spdiags(int n_row, int n_col, int n_diag, int offsets, long double diags,
+ std::vector<(int)> Ap, std::vector<(int)> Ai,
+ std::vector<(long double)> Ax)
+ spdiags(int n_row, int n_col, int n_diag, int offsets, npy_cfloat diags,
+ std::vector<(int)> Ap, std::vector<(int)> Ai,
+ std::vector<(npy_cfloat)> Ax)
+ spdiags(int n_row, int n_col, int n_diag, int offsets, npy_cdouble diags,
+ std::vector<(int)> Ap, std::vector<(int)> Ai,
+ std::vector<(npy_cdouble)> Ax)
+ spdiags(int n_row, int n_col, int n_diag, int offsets, npy_clongdouble diags,
+ std::vector<(int)> Ap, std::vector<(int)> Ai,
+ std::vector<(npy_clongdouble)> Ax)
+ """
+ return _sparsetools.spdiags(*args)
+
+def csrtodense(*args):
+ """
+ csrtodense(int n_row, int n_col, int Ap, int Aj, float Ax, float Mx)
+ csrtodense(int n_row, int n_col, int Ap, int Aj, double Ax, double Mx)
+ csrtodense(int n_row, int n_col, int Ap, int Aj, long double Ax,
+ long double Mx)
+ csrtodense(int n_row, int n_col, int Ap, int Aj, npy_cfloat Ax,
+ npy_cfloat Mx)
+ csrtodense(int n_row, int n_col, int Ap, int Aj, npy_cdouble Ax,
+ npy_cdouble Mx)
+ csrtodense(int n_row, int n_col, int Ap, int Aj, npy_clongdouble Ax,
+ npy_clongdouble Mx)
+ """
+ return _sparsetools.csrtodense(*args)
+
+def densetocsr(*args):
+ """
+ densetocsr(int n_row, int n_col, float Mx, std::vector<(int)> Ap,
+ std::vector<(int)> Aj, std::vector<(float)> Ax)
+ densetocsr(int n_row, int n_col, double Mx, std::vector<(int)> Ap,
+ std::vector<(int)> Aj, std::vector<(double)> Ax)
+ densetocsr(int n_row, int n_col, long double Mx, std::vector<(int)> Ap,
+ std::vector<(int)> Aj, std::vector<(long double)> Ax)
+ densetocsr(int n_row, int n_col, npy_cfloat Mx, std::vector<(int)> Ap,
+ std::vector<(int)> Aj, std::vector<(npy_cfloat)> Ax)
+ densetocsr(int n_row, int n_col, npy_cdouble Mx, std::vector<(int)> Ap,
+ std::vector<(int)> Aj, std::vector<(npy_cdouble)> Ax)
+ densetocsr(int n_row, int n_col, npy_clongdouble Mx, std::vector<(int)> Ap,
+ std::vector<(int)> Aj, std::vector<(npy_clongdouble)> Ax)
+ """
+ return _sparsetools.densetocsr(*args)
+
Deleted: trunk/Lib/sparse/sparsetools/sparsetools.pyf.src
===================================================================
--- trunk/Lib/sparse/sparsetools/sparsetools.pyf.src 2007-01-06 06:01:15 UTC (rev 2495)
+++ trunk/Lib/sparse/sparsetools/sparsetools.pyf.src 2007-01-06 06:34:41 UTC (rev 2496)
@@ -1,240 +0,0 @@
-! -*- f90 -*-
-python module sparsetools ! in
- interface ! in :sparsetools
- subroutine <_c>cscadd(n,a,rowa,ptra,nnzamax,b,rowb,ptrb,nnzbmax,c,rowc,ptrc,nnzcmax,ierr) ! in :sparsetools:dspblas.f
- integer intent(hide),depend(ptra) :: n=(len(ptra)-1)
- <_t> dimension(nnzamax) :: a
- integer dimension(nnzamax),depend(nnzamax) :: rowa
- integer dimension(n + 1) :: ptra
- integer intent(hide),depend(a) :: nnzamax=len(a)
- <_t> dimension(nnzbmax) :: b
- integer dimension(nnzbmax),depend(nnzbmax) :: rowb
- integer dimension(n + 1) :: ptrb
- integer intent(hide),depend(b) :: nnzbmax=len(b)
- <_t> intent(out), dimension(nnzcmax), depend(nnzcmax) :: c
- integer intent(out),dimension(nnzcmax),depend(nnzcmax) :: rowc
- integer intent(out),dimension(n + 1),depend(n) :: ptrc
- integer intent(hide),depend(nnzamax,nnzbmax) :: nnzcmax=nnzamax+nnzbmax
- integer intent(out) :: ierr
- end subroutine <_c>cscadd
-
- subroutine <_c>cscmul(n,a,rowa,ptra,nnzamax,b,rowb,ptrb,nnzbmax,c,rowc,ptrc,nnzcmax,ierr) ! in :sparsetools:dspblas.f
- integer intent(hide),depend(ptra) :: n=(len(ptra)-1)
- <_t> dimension(nnzamax) :: a
- integer dimension(nnzamax),depend(nnzamax) :: rowa
- integer dimension(n + 1) :: ptra
- integer intent(hide),depend(a) :: nnzamax=len(a)
- <_t> dimension(nnzbmax) :: b
- integer dimension(nnzbmax),depend(nnzbmax) :: rowb
- integer dimension(n + 1) :: ptrb
- integer intent(hide),depend(b) :: nnzbmax=len(b)
- <_t> intent(out), dimension(nnzcmax), depend(nnzcmax) :: c
- integer intent(out),dimension(nnzcmax),depend(nnzcmax) :: rowc
- integer intent(out),dimension(n + 1),depend(n) :: ptrc
- integer intent(hide),depend(nnzamax,nnzbmax) :: nnzcmax=nnzamax+nnzbmax
- integer intent(out) :: ierr
- end subroutine <_c>cscmul
-
- subroutine <_c>cscmux(a,rowa,ptra,nnzamax,ncol,x,mrow,y) ! in :sparsetools:dspblas.f
- <_t> dimension(nnzamax) :: a
- integer dimension(nnzamax),depend(nnzamax) :: rowa
- integer dimension(ncol + 1) :: ptra
- integer intent(hide),depend(a) :: nnzamax=len(a)
- integer intent(hide),depend(ptra) :: ncol=(len(ptra)-1)
- <_t> dimension(ncol),depend(ncol) :: x
- integer intent(in) :: mrow
- <_t> intent(out), dimension(mrow), depend(mrow) :: y
- end subroutine <_c>cscmux
-
- subroutine <_c>csrmux(a,cola,ptra,nnzamax,ncol,x,mrow,y) ! in :sparsetools:dspblas.f
- <_t> dimension(nnzamax) :: a
- integer dimension(nnzamax),depend(nnzamax) :: cola
- integer dimension(mrow + 1) :: ptra
- integer intent(hide),depend(a) :: nnzamax=len(a)
- integer intent(hide),depend(x) :: ncol=len(x)
- <_t> dimension(ncol) :: x
- integer intent(hide),depend(ptra) :: mrow=(len(ptra)-1)
- <_t> intent(out),dimension(mrow),depend(mrow) :: y
- end subroutine <_c>csrmux
-
- subroutine <_c>cscmucsr(m,k,n,a,rowa,ptra,nnzamax,b,colb,ptrb,nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr) ! in :sparsetools:dspblas.f
- integer :: m
- integer intent(hide),depend(ptra) :: k=(len(ptra)-1)
- integer intent(hide),depend(ptrc) :: n=(len(ptrc)-1)
- <_t> dimension(nnzamax) :: a
- integer dimension(nnzamax),depend(nnzamax) :: rowa
- integer dimension(k + 1) :: ptra
- integer intent(hide),depend(a) :: nnzamax=len(a)
- <_t> dimension(nnzbmax) :: b
- integer dimension(nnzbmax),depend(nnzbmax) :: colb
- integer dimension(k + 1),depend(k) :: ptrb
- integer intent(hide),depend(b) :: nnzbmax=len(b)
- <_t> intent(in, out), dimension(nnzcmax) :: c
- integer intent(in, out), dimension(nnzcmax),depend(nnzcmax) :: rowc
- integer intent(in, out), dimension(n + 1) :: ptrc
- integer intent(hide), depend(c) :: nnzcmax=len(c)
- integer intent(in, out) :: irow
- integer intent(in, out) :: kcol
- integer intent(in, out) :: ierr
- end subroutine <_c>cscmucsr
-
- subroutine <_c>csrmucsc(m,n,a,rowa,ptra,nnzamax,b,colb,ptrb,nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr) ! in :sparsetools:dspblas.f
- integer intent(hide), depend(ptra) :: m=(len(ptra)-1)
- integer intent(hide), depend(ptrb) :: n=(len(ptrb)-1)
- <_t> dimension(nnzamax) :: a
- integer dimension(nnzamax), depend(nnzamax) :: rowa
- integer dimension(m + 1) :: ptra
- integer intent(hide), depend(a) :: nnzamax=len(a)
- <_t> dimension(nnzbmax) :: b
- integer dimension(nnzbmax),depend(nnzbmax) :: colb
- integer dimension(n + 1) :: ptrb
- integer intent(hide), depend(b) :: nnzbmax=len(b)
- <_t> intent(in,out),dimension(nnzcmax) :: c
- integer intent(in, out), dimension(nnzcmax),depend(nnzcmax) :: rowc
- integer intent(in, out), dimension(n + 1), depend(n) :: ptrc
- integer intent(hide), depend(c) :: nnzcmax=len(c)
- integer intent(in, out) :: irow
- integer intent(in, out) :: kcol
- integer intent(in, out) :: ierr
- end subroutine <_c>csrmucsc
-
- subroutine <_c>cscmucsc(m,k,n,a,rowa,ptra,nnzamax,b,colb,ptrb,nnzbmax,c,rowc,ptrc,nnzcmax,irow,kcol,ierr) ! in :sparsetools:dspblas.f
- integer :: m
- integer intent(hid