updates

base/PyNucleus_base/CSR_LinearOperator_{SCALAR}.pxi

@@ -517,6 +517,31 @@ cdef class {SCALAR_label}CSR_VectorLinearOperator({SCALAR_label}VectorLinearOper
                     y[i, k] += self.data[jj, k]*x[j]
         return 0
 
+    cdef INDEX_t matvecTrans({SCALAR_label}CSR_VectorLinearOperator self,
+                             {SCALAR}_t[::1] x,
+                             {SCALAR}_t[:, ::1] y) except -1:
+        cdef:
+            INDEX_t i, j, jj, k
+        y[:, :] = 0.
+        for i in range(self.indptr.shape[0]-1):
+            for jj in range(self.indptr[i], self.indptr[i+1]):
+                j = self.indices[jj]
+                for k in range(self.data.shape[1]):
+                    y[j, k] += self.data[jj, k]*x[i]
+        return 0
+
+    cdef INDEX_t matvecTrans_no_overwrite({SCALAR_label}CSR_VectorLinearOperator self,
+                                          {SCALAR}_t[::1] x,
+                                          {SCALAR}_t[:, ::1] y) except -1:
+        cdef:
+            INDEX_t i, j, jj, k
+        for i in range(self.indptr.shape[0]-1):
+            for jj in range(self.indptr[i], self.indptr[i+1]):
+                j = self.indices[jj]
+                for k in range(self.data.shape[1]):
+                    y[j, k] += self.data[jj, k]*x[i]
+        return 0
+
     def isSparse(self):
         return True
 

base/PyNucleus_base/LinearOperator_decl_{SCALAR}.pxi

@@ -110,6 +110,12 @@ cdef class {SCALAR_label}VectorLinearOperator:
     cdef INDEX_t matvec_no_overwrite(self,
                                      {SCALAR}_t[::1] x,
                                      {SCALAR}_t[:, ::1] y) except -1
+    cdef INDEX_t matvecTrans(self,
+                             {SCALAR}_t[::1] x,
+                             {SCALAR}_t[:, ::1] y) except -1
+    cdef INDEX_t matvecTrans_no_overwrite(self,
+                                          {SCALAR}_t[::1] x,
+                                          {SCALAR}_t[:, ::1] y) except -1
     cdef void addToEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val)
     cdef void getEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val)
     cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val)

base/PyNucleus_base/LinearOperator_{SCALAR}.pxi

@@ -46,7 +46,10 @@ cdef class {SCALAR_label}LinearOperator:
             else:
                 self.matvec(x, y)
         else:
-            self.matvecTrans(x, y)
+            if no_overwrite:
+                self.matvecTrans_no_overwrite(x, y)
+            else:
+                self.matvecTrans(x, y)
 
     def dot(self, {SCALAR}_t[::1] x):
         cdef:
@@ -545,14 +548,31 @@ cdef class {SCALAR_label}VectorLinearOperator:
                                      {SCALAR}_t[:, ::1] y) except -1:
         return -1
 
+    cdef INDEX_t matvecTrans(self,
+                             {SCALAR}_t[::1] x,
+                             {SCALAR}_t[:, ::1] y) except -1:
+        return -1
+
+    cdef INDEX_t matvecTrans_no_overwrite(self,
+                                          {SCALAR}_t[::1] x,
+                                          {SCALAR}_t[:, ::1] y) except -1:
+        return -1
+
     def __call__(self,
                  {SCALAR}_t[::1] x,
                  {SCALAR}_t[:, ::1] y,
-                 BOOL_t no_overwrite=False):
-        if no_overwrite:
-            self.matvec_no_overwrite(x, y)
+                 BOOL_t no_overwrite=False,
+                 BOOL_t trans=False):
+        if not trans:
+            if no_overwrite:
+                self.matvec_no_overwrite(x, y)
+            else:
+                self.matvec(x, y)
         else:
-            self.matvec(x, y)
+            if no_overwrite:
+                self.matvecTrans_no_overwrite(x, y)
+            else:
+                self.matvecTrans(x, y)
 
     def __mul__(self, x):
         cdef:

base/PyNucleus_base/SSS_LinearOperator_{SCALAR}.pxi

@@ -392,6 +392,16 @@ cdef class {SCALAR_label}SSS_VectorLinearOperator({SCALAR_label}VectorLinearOper
                 y[i, l] += self.temp[l]
         return 0
 
+    cdef INDEX_t matvecTrans({SCALAR_label}SSS_VectorLinearOperator self,
+                             {SCALAR}_t[::1] x,
+                             {SCALAR}_t[:, ::1] y) except -1:
+        return self.matvec(x, y)
+
+    cdef INDEX_t matvecTrans_no_overwrite({SCALAR_label}SSS_VectorLinearOperator self,
+                                          {SCALAR}_t[::1] x,
+                                          {SCALAR}_t[:, ::1] y) except -1:
+        return self.matvec_no_overwrite(x, y)
+
     cdef void setEntry({SCALAR_label}SSS_VectorLinearOperator self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val):
         cdef:
             INDEX_t i, low, mid, high, l

base/PyNucleus_base/intTuple.pyx

@@ -7,7 +7,7 @@
 
 import numpy as np
 cimport numpy as np
-from libc.stdlib cimport malloc
+from cpython.mem cimport PyMem_Malloc
 from libc.string cimport memcpy
 from . myTypes import INDEX
 
@@ -19,7 +19,7 @@ cdef enum:
 cdef class intTuple:
     cdef void set(self, INDEX_t * t, int size):
         self.size = size
-        self.entries = <INDEX_t *>malloc(size*INDEX_SIZE)
+        self.entries = <INDEX_t *>PyMem_Malloc(size*INDEX_SIZE)
         memcpy(&self.entries[0], &t[0], size*INDEX_SIZE)
 
     cdef void assign(self, INDEX_t * t):
@@ -55,7 +55,7 @@ cdef class intTuple:
         cdef:
             intTuple t = intTuple()
         t.size = 2
-        t.entries = <INDEX_t *>malloc(2*INDEX_SIZE)
+        t.entries = <INDEX_t *>PyMem_Malloc(2*INDEX_SIZE)
         t.entries[0] = a
         t.entries[1] = b
         return t
@@ -69,7 +69,7 @@ cdef class intTuple:
         cdef:
             intTuple t = intTuple()
         t.size = 3
-        t.entries = <INDEX_t *>malloc(3*INDEX_SIZE)
+        t.entries = <INDEX_t *>PyMem_Malloc(3*INDEX_SIZE)
         t.entries[0] = a
         t.entries[1] = b
         t.entries[2] = c

base/PyNucleus_base/opt_false_blas.pxi

@@ -54,7 +54,7 @@ cdef void scaleScalar(SCALAR_t[::1] x, SCALAR_t alpha):
         for i in range(x.shape[0]):
             x[i] = x[i]*alpha
 
-cdef SCALAR_t mydot(SCALAR_t[::1] v0, SCALAR_t[::1] v1) nogil:
+cdef SCALAR_t mydot(SCALAR_t[::1] v0, SCALAR_t[::1] v1):
     cdef:
         int i
         SCALAR_t s = 0.0
@@ -110,6 +110,25 @@ cdef void gemv(SCALAR_t[:, ::1] A, SCALAR_t[::1] x, SCALAR_t[::1] y, SCALAR_t be
                 y[i] = s
 
 
+cdef void gemvF(SCALAR_t[::1, :] A, SCALAR_t[::1] x, SCALAR_t[::1] y, SCALAR_t beta=0.):
+    cdef:
+        INDEX_t i, j
+    if SCALAR_t is COMPLEX_t:
+        if beta != 0.:
+            for i in range(A.shape[0]):
+                y[i] = beta*y[i]
+        for j in range(A.shape[1]):
+            for i in range(A.shape[0]):
+                y[i] = y[i] + A[i, j]*x[j]
+    else:
+        if beta != 0.:
+            for i in range(A.shape[0]):
+                y[i] *= beta
+        for j in range(A.shape[1]):
+            for i in range(A.shape[0]):
+                y[i] += A[i, j]*x[j]
+
+
 cdef void gemvT(SCALAR_t[:, ::1] A, SCALAR_t[::1] x, SCALAR_t[::1] y, SCALAR_t beta=0.):
     cdef:
         INDEX_t i, j
@@ -119,17 +138,17 @@ cdef void gemvT(SCALAR_t[:, ::1] A, SCALAR_t[::1] x, SCALAR_t[::1] y, SCALAR_t b
                 y[i] = y[i]*beta
         else:
             y[:] = 0.
-        for j in range(A.shape[1]):
-            for i in range(A.shape[0]):
+        for i in range(A.shape[1]):
+            for j in range(A.shape[0]):
                 y[i] = y[i]+A[j, i]*x[j]
     else:
         if beta != 0.:
             for i in range(A.shape[0]):
                 y[i] *= beta
         else:
             y[:] = 0.
-        for j in range(A.shape[1]):
-            for i in range(A.shape[0]):
+        for i in range(A.shape[1]):
+            for j in range(A.shape[0]):
                 y[i] += A[j, i]*x[j]
 
 

base/PyNucleus_base/opt_true_blas.pxi

@@ -215,10 +215,10 @@ cdef void gemvT(SCALAR_t[:, ::1] A, SCALAR_t[::1] x, SCALAR_t[::1] y, SCALAR_t b
         double complex* A_ptr_c
         double complex* x_ptr_c
         double complex* y_ptr_c
-        int m = A.shape[0]
-        int n = A.shape[1]
+        int m = A.shape[1]
+        int n = A.shape[0]
         SCALAR_t alpha = 1.
-        int lda = A.shape[0]
+        int lda = A.shape[1]
         int incx = 1
         int incy = 1
     if SCALAR_t is COMPLEX_t:

base/PyNucleus_base/sparsityPattern.pyx

@@ -5,7 +5,7 @@
 # If you want to use this code, please refer to the README.rst and LICENSE files. #
 ###################################################################################
 
-from libc.stdlib cimport malloc, realloc, free
+from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
 import numpy as np
 cimport numpy as np
 from . myTypes import INDEX
@@ -21,10 +21,10 @@ cdef class sparsityPattern:
         self.nnz = 0
         self.counts = np.zeros((num_dofs), dtype=INDEX)
         self.lengths = initial_length*np.ones((num_dofs), dtype=np.uint16)
-        self.indexL = <INDEX_t **>malloc(num_dofs*sizeof(INDEX_t *))
+        self.indexL = <INDEX_t **>PyMem_Malloc(num_dofs*sizeof(INDEX_t *))
         # reserve initial memory for array of variable column size
         for i in range(num_dofs):
-            self.indexL[i] = <INDEX_t *>malloc(self.initial_length*sizeof(INDEX_t))
+            self.indexL[i] = <INDEX_t *>PyMem_Malloc(self.initial_length*sizeof(INDEX_t))
 
     cdef inline BOOL_t findIndex(self, INDEX_t I, INDEX_t J):
         cdef:
@@ -68,7 +68,7 @@ cdef class sparsityPattern:
             # J was not present
             # Do we need more space?
             if self.counts[I] == self.lengths[I]:
-                self.indexL[I] = <INDEX_t *>realloc(self.indexL[I], (self.lengths[I]+self.initial_length)*sizeof(INDEX_t))
+                self.indexL[I] = <INDEX_t *>PyMem_Realloc(self.indexL[I], (self.lengths[I]+self.initial_length)*sizeof(INDEX_t))
                 self.lengths[I] += self.initial_length
             # where should we insert?
             m = self.index
@@ -95,8 +95,8 @@ cdef class sparsityPattern:
             for j in range(self.counts[i]):
                 indices[k] = self.indexL[i][j]
                 k += 1
-            free(self.indexL[i])
-        free(self.indexL)
+            PyMem_Free(self.indexL[i])
+        PyMem_Free(self.indexL)
 
         # fill indptr array
         indptr_mem = uninitialized((self.num_dofs+1), dtype=INDEX)

base/PyNucleus_base/tupleDict.pxd

@@ -15,6 +15,7 @@ ctypedef np.uint64_t MEM_t
 
 cdef class indexSet:
     cdef BOOL_t inSet(self, INDEX_t i)
+    cdef INDEX_t position(self, INDEX_t i)
     cpdef void fromSet(self, set s)
     cpdef set toSet(self)
     cpdef INDEX_t[::1] toArray(self)

base/PyNucleus_base/tupleDict.pyx

@@ -6,7 +6,7 @@
 ###################################################################################
 
 import numpy as np
-from libc.stdlib cimport malloc, realloc, free
+from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
 from libc.stdlib cimport qsort
 from . myTypes import INDEX
 
@@ -25,6 +25,9 @@ cdef class indexSet:
     def inSet_py(self, INDEX_t i):
         return self.inSet(i)
 
+    cdef INDEX_t position(self, INDEX_t i):
+        raise NotImplementedError()
+
     cpdef void fromSet(self, set s):
         raise NotImplementedError()
 
@@ -176,6 +179,27 @@ cdef class arrayIndexSet(indexSet):
                     high = mid
             return True
 
+    cdef INDEX_t position(self, INDEX_t i):
+        cdef:
+            INDEX_t low = 0
+            INDEX_t high = self.indexArray.shape[0]
+            INDEX_t mid
+        if high-low < 20:
+            for mid in range(low, high):
+                if self.indexArray[mid] == i:
+                    return mid
+            return -1
+        else:
+            while self.indexArray[low] != i:
+                if high-low <= 1:
+                    return -1
+                mid = (low+high) >> 1
+                if self.indexArray[mid] <= i:
+                    low = mid
+                else:
+                    high = mid
+            return low
+
     cpdef void fromSet(self, set s):
         cdef:
             INDEX_t i, k
@@ -385,6 +409,14 @@ cdef class unsortedArrayIndexSet(arrayIndexSet):
                 return True
         return False
 
+    cdef INDEX_t position(self, INDEX_t i):
+        cdef:
+            INDEX_t j
+        for j in range(self.indexArray.shape[0]):
+            if self.indexArray[j] == i:
+                return j
+        return -1
+
     cpdef void fromSet(self, set s):
         cdef:
             INDEX_t i, k
@@ -446,7 +478,7 @@ cdef class arrayIndexSetIterator(indexSetIterator):
 cdef class bitArray(indexSet):
     def __init__(self, size_t hintMaxLength=1, INDEX_t maxElement=0):
         self.length = max(hintMaxLength, maxElement/(sizeof(MEM_t)*8)+1)
-        self.a = <MEM_t *>malloc(self.length*sizeof(MEM_t))
+        self.a = <MEM_t *>PyMem_Malloc(self.length*sizeof(MEM_t))
         for j in range(self.length):
             self.a[j] = 0
 
@@ -459,7 +491,7 @@ cdef class bitArray(indexSet):
         if k >= self.length:
             l = self.length
             self.length = k+1
-            self.a = <MEM_t *>realloc(self.a, self.length * sizeof(MEM_t))
+            self.a = <MEM_t *>PyMem_Realloc(self.a, self.length * sizeof(MEM_t))
             for j in range(l, self.length):
                 self.a[j] = 0
         self.a[k] |= one << n
@@ -510,7 +542,7 @@ cdef class bitArray(indexSet):
             self.a[j] = 0
 
     def __dealloc__(self):
-        free(self.a)
+        PyMem_Free(self.a)
 
     cdef indexSetIterator getIter(self):
         return bitArrayIterator(self)