Ticket #12177: trac_12177-coeff_matrices_template.patch

File trac_12177-coeff_matrices_template.patch, 29.6 KB (added by burcin, 10 years ago)
  • module_list.py

    # HG changeset patch
    # User Burcin Erocal <burcin@erocal.org>
    # Date 1324468679 -3600
    
    diff --git a/module_list.py b/module_list.py
    a b  
    951951              libraries = ['gmp', 'linbox', BLAS, BLAS2],
    952952              extra_compile_args = ['-DDISABLE_COMMENTATOR']),
    953953
     954    Extension('sage.matrix.matrix_modq_dense_float',
     955              sources = ['sage/matrix/matrix_modq_dense_float.pyx'],
     956              language="c++",
     957              libraries = ['linbox', BLAS, BLAS2],
     958              extra_compile_args = ['-DDISABLE_COMMENTATOR']),
     959
     960    Extension('sage.matrix.matrix_modq_dense_double',
     961              sources = ['sage/matrix/matrix_modq_dense_double.pyx'],
     962              language="c++",
     963              libraries = ['linbox', BLAS, BLAS2],
     964              extra_compile_args = ['-DDISABLE_COMMENTATOR']),
     965
    954966    Extension('sage.matrix.matrix_modn_sparse',
    955967              sources = ['sage/matrix/matrix_modn_sparse.pyx'],
    956968              libraries = ['gmp']),
  • sage/libs/linbox/fflas.pxd

    diff --git a/sage/libs/linbox/fflas.pxd b/sage/libs/linbox/fflas.pxd
    a b  
    5555                    ModDoubleFieldElement beta, ModDoubleFieldElement* C,
    5656                    size_t C_stride)
    5757
     58    void ModDouble_faxpy "LinBox::FFLAS::faxpy"  (ModDoubleField F,
     59            size_t N, ModDoubleFieldElement alpha,
     60            ModDoubleFieldElement* X, size_t incX,
     61            ModDoubleFieldElement* Y, size_t incY)
     62
     63    #ModDoubleElement * ModDouble_fsquare "LinBox::FFLAS::fsquare" \
     64    #        (ModDoubleField F, fflas_trans_enum transA, size_t n,
     65    #        ModDoubleElement alpha, ModDoubleElement* A, size_t lda,
     66    #        ModDoubleElement beta, ModDoubleElement* C, size_t ldc)
     67
     68    #void ModDouble_ftrmm "Linbox::FFLAS::ftrmm" (ModDoubleField F,
     69    #        int, #const FFLAS_SIDE    Side,
     70    #        int, #const FFLAS_UPLO    Uplo,
     71    #        fflas_trans_enum transa,
     72    #        int, #const FFLAS_DIAG    Diag,
     73    #        size_t M, size_t N, ModDoubleFieldElement alpha,
     74    #        ModDoubleFieldElement* A, size_t lda,
     75    #        ModDoubleFieldElement* B, size_t ldb)
    5876
    5977    # float
    6078    void ModFloat_fgemv "LinBox::FFLAS::fgemv" \
     
    7492                    ModFloatFieldElement beta, ModFloatFieldElement* C,
    7593                    size_t C_stride)
    7694
     95    void ModFloat_faxpy "LinBox::FFLAS::faxpy"  (ModFloatField F,
     96            size_t N, ModFloatFieldElement alpha,
     97            ModFloatFieldElement* X, size_t incX,
     98            ModFloatFieldElement* Y, size_t incY)
     99
     100    #ModFloatElement * ModFloat_fsquare "LinBox::FFLAS::fsquare" \
     101    #        (ModFloatField F, fflas_trans_enum transA, size_t n,
     102    #        ModFloatElement alpha, ModFloatElement* A, size_t lda,
     103    #        ModFloatElement beta, ModFloatElement* C, size_t ldc)
     104
     105    #void ModFloat_ftrmm "Linbox::FFLAS::ftrmm" (ModFloatField F,
     106    #        int, #const FFLAS_SIDE    Side,
     107    #        int, #const FFLAS_UPLO    Uplo,
     108    #        fflas_trans_enum transa,
     109    #        int, #const FFLAS_DIAG    Diag,
     110    #        size_t M, size_t N, ModFloatFieldElement alpha,
     111    #        ModFloatFieldElement* A, size_t lda,
     112    #        ModFloatFieldElement* B, size_t ldb)
     113
    77114cdef extern from "linbox/ffpack/ffpack.h":
    78115    # double
    79116    bint ModDouble_is_singular "LinBox::FFPACK::IsSingular" (ModDoubleField F,
  • new file sage/matrix/matrix_dense_slice_template.pxi

    diff --git a/sage/matrix/matrix_dense_slice_template.pxi b/sage/matrix/matrix_dense_slice_template.pxi
    new file mode 100644
    - +  
     1###############################################################################
     2#   Sage: Open Source Mathematical Software
     3#       Copyright (C) 2010 Burcin Erocal <burcin@erocal.org>
     4#  Distributed under the terms of the GNU General Public License (GPL),
     5#  version 2 or any later version.  The full text of the GPL is available at:
     6#                  http://www.gnu.org/licenses/
     7###############################################################################
     8"""
     9This is a template for dense matrices which can be represented as a polynomial
     10with matrix coefficients whose entries are C types such as int, float. Some examples of these structures are matrices over GF(p^k), or GF(p)[x].
     11
     12
     13# FIXME: this needs to be updated to mention the matrix modn dense templates
     14#
     15#The templating system works in a similar fashion to the template for polynomial
     16#classes provided in ``sage/rings/polynomial/polynomial_template.pxi``. In
     17#order to use this template, write a wrapper file similar to
     18#``sage/libs/flint/fmpz_poly_linkage.pxi``, defining the basic functions
     19#used to manipulate the underlying C structs. In your `pyx` file defining the
     20#element class for your object, first include the linkage file, then this
     21#template file. In the corresponding pxd file, include the header template
     22#``sage/matrix/matrix_dense_template_header.pxi``.
     23
     24#See ``sage/matrix/matrix_ZZx_dense.pyx`` and
     25#``sage/matrix/matrix_ZZx_dense.pxd`` for an example.
     26"""
     27
     28from sage.structure.element cimport Element
     29 
     30cdef class Matrix_dense_slice_template(matrix_dense.Matrix_dense):
     31    """
     32    TESTS::
     33
     34    """
     35    def __cinit__(self, parent, entries, copy=True, coerce=True):
     36        """
     37        Allocate memory for the matrix.
     38
     39        TESTS::
     40
     41        """
     42        self._nrows = parent.nrows()
     43        self._ncols = parent.ncols()
     44        self._char = parent.base_ring().characteristic()
     45        self._degree = parent.base_ring().degree()
     46
     47        cdef Py_ssize_t i, j
     48
     49        # allocate space for the entries
     50        _sig_on
     51        self._coeffs_block = <celement*>sage_malloc(sizeof(celement)*self._degree*self._nrows*self._ncols)
     52        if self._coeffs_block == NULL:
     53            raise MemoryError, "out of memory allocating a matrix"
     54        self._coeffs = <celement**>sage_malloc(self._degree*sizeof(celement*))
     55        if self._coeffs == NULL:
     56            raise MemoryError, "out of memory allocating a matrix"
     57        for i in range(self._degree):
     58            self._coeffs[i] = self._coeffs_block + i*self._nrows*self._ncols
     59        _sig_off
     60
     61        for i in range(self._degree*self._nrows*self._ncols):
     62                self._coeffs_block[i] = <celement>0
     63
     64    def __init__(self, parent, entries, copy=True, coerce=True):
     65        """
     66        TESTS::
     67
     68        """
     69        cdef Py_ssize_t i
     70        matrix_dense.Matrix_dense.__init__(self, parent)
     71        cdef celement* tvector
     72        tvector = <celement*>sage_malloc(sizeof(celement)*self._degree)
     73
     74        if entries == 0:
     75            pass
     76
     77        elif isinstance(entries, (list, tuple)):
     78            if len(entries) != self._nrows*self._ncols:
     79                raise ValueError, "entries has wrong length"
     80            i = 0
     81            for e in entries:
     82                if coerce: # the argument 'coerce' is meant to be convert
     83                    selement = self._base_ring(e)
     84                else:
     85                    selement = e
     86
     87                sage_element_to_celement_vector(tvector, selement,
     88                        self._char, self._degree)
     89                for j in range(self._degree):
     90                    self._coeffs[j][i] = tvector[j]
     91                i += 1
     92        else:
     93            if self._ncols != self._nrows:
     94                raise ValueError, "matrix is not square"
     95            if coerce:
     96                selement = self._base_ring.coerce(entries)
     97            else:
     98                selement = entries
     99            sage_element_to_celement_vector(tvector, selement,
     100                    self._char, self._degree)
     101            for i in range(self._degree):
     102                for j in range(self._nrows):
     103                    self._coeffs[i][j*self._nrows + j] = tvector[i]
     104
     105        sage_free(tvector)
     106
     107
     108    def __dealloc__(self):
     109        """
     110        TESTS::
     111
     112            sage: from sage.matrix.matrix_modq_dense_float import Matrix_modq_dense_float
     113            sage: K.<a> = GF(5^2)
     114            sage: MS = MatrixSpace(K, 2000, 2000)
     115            sage: M = Matrix_modq_dense_float(MS, 1)
     116            sage: M.randomize()
     117            sage: del M
     118        """
     119        cdef Py_ssize_t i
     120        if self._coeffs_block != NULL:
     121            sage_free(self._coeffs_block)
     122            sage_free(self._coeffs)
     123
     124    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, x) except +:
     125        # no bounds or any other checks; assumes x is in self._base_ring
     126        cdef Py_ssize_t k
     127        cdef celement* tvector
     128        tvector = <celement*>sage_malloc(sizeof(celement)*self._degree)
     129        sage_element_to_celement_vector(tvector, x, self._char, self._degree)
     130
     131        for k in range(self._degree):
     132            self._coeffs[k][i*self._nrows + j] = tvector[k]
     133        sage_free(tvector)
     134
     135    cdef object get_unsafe(Matrix_dense_slice_template self,
     136            Py_ssize_t i, Py_ssize_t j) except +:
     137        """
     138        """
     139        # no checks
     140        cdef Py_ssize_t k
     141        res = [int(self._coeffs[k][self._nrows*i + j])
     142                for k in range(self._degree)]
     143        return self.base_ring()(res)
     144
     145    def _print_coeff_mat(self, deg):
     146        for i in range(self._nrows):
     147            s= ' '.join([str(self._coeffs[deg][i*self._nrows + j])
     148                    for j in range(self._ncols)])
     149            print s
     150
     151
     152cdef class Matrix_dense_slice_template_GFq(Matrix_dense_slice_template):
     153    def __init__(self, parent, entries, copy=True, coerce=True):
     154        Matrix_dense_slice_template.__init__(self, parent, entries, copy, coerce)
     155        # these should really be a part of the parent
     156        # unfortunately, matrix spaces are implemented in python and we cannot
     157        # subclass them and attach C objects to these
     158        self._linbox_field = new ModField(<long>self._char)
     159        poly = self.base_ring().polynomial()
     160        self._minpoly_vector  = <ModFieldElement*> \
     161                sage_malloc(sizeof(ModFieldElement)*self._degree)
     162        for i in range(self._degree):
     163            print "i: {i}, coeff: {coeff}".format(i=i, coeff=poly[i])
     164            t = -poly[i]
     165            print "t:",t
     166            self._linbox_field[0].init(self._minpoly_vector[i], <int>t)
     167
     168        # stores constants for zero and one in that order
     169        self._linbox_consts = <ModFieldElement*>\
     170                sage_malloc(sizeof(ModFieldElement)*2)
     171        self._linbox_field[0].init(self._linbox_consts[0], <int>0)
     172        self._linbox_field[0].init(self._linbox_consts[1], <int>1)
     173
     174    def __dealloc__(self):
     175        sage_free(self._minpoly_vector)
     176        sage_free(self._linbox_consts)
     177        del self._linbox_field
     178
     179    def _print_minpoly_vec(self):
     180        s= ' '.join([str(<celement>self._minpoly_vector[i]) \
     181                for i in range(self._degree)])
     182        print s
     183
     184    def _multiply_toom(self, Matrix_dense_slice_template other):
     185        cdef Py_ssize_t i, j, k
     186        cdef Matrix_dense_slice_template res = self.__class__(self.parent(), 0)
     187
     188        cdef Py_ssize_t l = 2*self._degree - 1
     189
     190        # allocate temporary memory
     191        cdef celement* res_coeff_matrices_block = <celement*>sage_malloc(
     192                sizeof(celement)*l*self._nrows*self._ncols)
     193        if res_coeff_matrices_block == NULL:
     194            raise MemoryError("cannot allocate temporary memory")
     195        cdef celement** res_coeff_matrices = <celement**>sage_malloc(
     196                sizeof(celement*)*l)
     197        if res_coeff_matrices == NULL:
     198            sage_free(res_coeff_matrices_block)
     199            raise MemoryError("cannot allocate temporary memory")
     200        for i in range(l):
     201            res_coeff_matrices[i] = res_coeff_matrices_block + \
     202                    i*self._nrows*self._ncols
     203
     204        # clear temporary matrices
     205        #for i in range(l):
     206        #    for j in range(self._nrows*self._ncols):
     207        #        res_coeff_matrices[i][j] = <celement>0
     208       
     209        cdef ModFieldElement* eval_helper = <ModFieldElement*>sage_malloc(
     210                sizeof(ModFieldElement)*l*l)
     211        if eval_helper == NULL:
     212            sage_free(res_coeff_matrices_block)
     213            # FIXME free res_coeff_matrices
     214            raise MemoryError("cannot allocate temporary memory")
     215
     216        cdef ModFieldElement* eval_helper_inverse = \
     217                <ModFieldElement*>sage_malloc(sizeof(ModFieldElement)*l*l)
     218        if eval_helper_inverse == NULL:
     219            sage_free(res_coeff_matrices_block)
     220            sage_free(res_coeff_matrices)
     221            sage_free(eval_helper)
     222            raise MemoryError("cannot allocate temporary memory")
     223
     224        # FWD = Matrix(K, l, l, [K(i**j) for i in range(l) for j in range(l)])
     225        for i in range(l):
     226            for j in range(l):
     227                self._linbox_field[0].init(eval_helper[i*l + j],
     228                        <int>i**j)
     229
     230        #print "eval_helper:"
     231        #for q in range(l):
     232        #    s= ' '.join([str(<celement>(eval_helper[q*l + m]))
     233        #            for m in range(l)])
     234        #    print s
     235
     236        cdef celement* self_evaluated = <celement*>sage_malloc(
     237                sizeof(celement)*self._nrows*self._ncols)
     238        if self_evaluated == NULL:
     239            sage_free(res_coeff_matrices_block)
     240            # FIXME: free all eval_helper and res_coeff_matrices
     241            raise MemoryError("cannot allocate temporary memory")
     242
     243        cdef celement* other_evaluated = <celement*>sage_malloc(
     244                sizeof(celement)*self._nrows*self._ncols)
     245        if other_evaluated == NULL:
     246            sage_free(res_coeff_matrices_block)
     247            # FIXME: free self_evaluated, eval_helper and res_coeff_matrices
     248            raise MemoryError("cannot allocate temporary memory")
     249
     250        # AY = [sum(FWD[i,j]*A[j] for j in range(len(A))) for i in range(l)]
     251        # BY = [sum(FWD[i,j]*B[j] for j in range(len(B))) for i in range(l)]
     252        # Y = [AY[i]*BY[i] for i in range(l)]
     253        # these 3 operations broken to a loop over i
     254        for i in range(l):
     255            #for j in range(self._nrows*self._ncols):
     256            #    self_evaluated[j] = <celement>0
     257            #for j in range(self._nrows*self._ncols):
     258            #    other_evaluated[j] = <celement>0
     259            Mod_fgemv(self._linbox_field[0], fflas_trans,
     260                        self._degree, self._nrows*self._ncols,
     261                        self._linbox_consts[1],
     262                        <ModFieldElement*>self._coeffs_block,
     263                        self._nrows*self._ncols,
     264                        eval_helper + i*l, 1,
     265                        self._linbox_consts[0],
     266                        <ModFieldElement*>self_evaluated, 1)
     267
     268            Mod_fgemv(self._linbox_field[0], fflas_trans,
     269                        self._degree, self._nrows*self._ncols,
     270                        self._linbox_consts[1],
     271                        <ModFieldElement*>other._coeffs_block,
     272                        self._nrows*self._ncols,
     273                        eval_helper + i*l, 1,
     274                        self._linbox_consts[0],
     275                        <ModFieldElement*>other_evaluated, 1)
     276
     277            Mod_fgemm(self._linbox_field[0],
     278                    fflas_no_trans, fflas_no_trans,
     279                    self._nrows, self._ncols, self._nrows,
     280                    self._linbox_consts[1],
     281                    <ModFieldElement*>self_evaluated, self._nrows,
     282                    <ModFieldElement*>other_evaluated, self._nrows,
     283                    self._linbox_consts[0],
     284                    <ModFieldElement*>res_coeff_matrices[i],
     285                    self._nrows)
     286
     287        # BCK = ~FWD
     288        cdef int tmp=0
     289        Mod_invert_in_place(self._linbox_field[0], l,
     290                <ModFieldElement*>eval_helper, l, tmp)
     291
     292        #print "eval_helper inverse:"
     293        #for q in range(l):
     294        #    s= ' '.join([str(<celement>(eval_helper[q*l + m]))
     295        #            for m in range(l)])
     296        #    print s
     297
     298        # construct a temporary matrix in the memory allocated for the result
     299        # to represent the rollover for higher degree products
     300        for i in range(self._degree):
     301            res._coeffs_block[i*l + i] = <celement>1
     302        for j in range(self._degree-1):
     303            res._coeffs_block[j*l + self._degree] = <celement>self._minpoly_vector[j]
     304        cdef ModFieldElement coeff
     305        cdef ModFieldElement tmp_coeff, tmp_coeff2
     306        for i in range(1, self._degree-1):
     307            for j in range(1, self._degree+1):
     308                res._coeffs_block[j*l + i + self._degree] = res._coeffs_block[(j-1)*l + i - 1 + self._degree]
     309            if not self._linbox_field[0].isZero(<ModFieldElement>
     310                    res._coeffs_block[self._degree*l + i + self._degree]):
     311                #print "nonzero"
     312                coeff = <ModFieldElement>res._coeffs_block[
     313                        self._degree*l + i + self._degree]
     314                #print "coeff:",<celement>coeff
     315                for j in range(self._degree-1):
     316                    self._linbox_field[0].mul(tmp_coeff, coeff,
     317                            self._minpoly_vector[j])
     318                    #print "assigning:", <celement>tmp_coeff
     319                    tmp_coeff2 = <ModFieldElement>(res._coeffs_block[
     320                                j*l + i + self._degree])
     321                    self._linbox_field[0].addin(tmp_coeff2,
     322                            tmp_coeff)
     323
     324                    res._coeffs_block[j*l + i + self._degree] = \
     325                            <celement>tmp_coeff2
     326
     327        #print "eval_helper helper:"
     328        #for q in range(self._degree):
     329        #    s= ' '.join([str(<celement>(res._coeffs_block[q*l + m]))
     330        #            for m in range(l)])
     331        #    print s
     332
     333        Mod_fgemm(self._linbox_field[0],
     334                fflas_no_trans, fflas_no_trans,
     335                l, l, l, self._linbox_consts[1],
     336                <ModFieldElement*>res._coeffs_block, l,
     337                <ModFieldElement*>eval_helper, l,
     338                self._linbox_consts[0],
     339                <ModFieldElement*>eval_helper_inverse, l)
     340
     341        #print "eval_helper_inverse with rollback:"
     342        #for q in range(self._degree):
     343        #    s= ' '.join([str(<celement>(eval_helper_inverse[q*l + m]))
     344        #            for m in range(l)])
     345        #    print s
     346
     347
     348        # Y = [sum(BCK[i,j]*Y[j] for j in range(l)) for i in range(l)]
     349        Mod_fgemm(self._linbox_field[0],
     350                fflas_no_trans, fflas_no_trans,
     351                self._degree, self._nrows*self._ncols, l,
     352                self._linbox_consts[1],
     353                <ModFieldElement*>eval_helper_inverse, l,
     354                <ModFieldElement*>res_coeff_matrices_block,
     355                self._nrows*self._ncols,
     356                self._linbox_consts[0],
     357                <ModFieldElement*>res._coeffs_block,
     358                self._nrows*self._ncols)
     359
     360        #print "res deg 0:"
     361        #res._print_coeff_mat(0)
     362        #print "res deg 1:"
     363        #res._print_coeff_mat(1)
     364
     365        sage_free(res_coeff_matrices_block)
     366        sage_free(res_coeff_matrices)
     367        sage_free(self_evaluated)
     368        sage_free(other_evaluated)
     369        sage_free(eval_helper)
     370        sage_free(eval_helper_inverse)
     371
     372        return res
     373
     374    def _multiply_classical(self, Matrix_dense_slice_template other):
     375        """
     376            sage: from sage.matrix.matrix_modq_dense_float import Matrix_modq_dense_float
     377            sage: K.<a> = GF(5^2)
     378            sage: MS = MatrixSpace(K, 3, 3)
     379            sage: M = Matrix_modq_dense_float(MS, 1)
     380            sage: M*M
     381
     382
     383            sage: p = 17
     384            sage: K.<a> = GF(p^2)
     385            sage: MS = MatrixSpace(K, 2000, 2000)
     386            age: M = Matrix_modq_dense_float(MS, 1)
     387
     388        """
     389        cdef Py_ssize_t i, j, k
     390        cdef Matrix_dense_slice_template res = self.__class__(self.parent(), 0)
     391
     392        cdef ModFieldElement coeff
     393
     394        #print "before allocating temp"
     395
     396        # temporary matrix to store results of coefficient matrix
     397        # multiplications
     398        cdef celement* tmp_matrix= <celement*>sage_malloc(\
     399                sizeof(celement)*self._nrows*self._ncols)
     400        # pointer we'll pass to FFLAS calls to specify where to put the result
     401        cdef celement* res_matrix
     402        # coefficient of result
     403        # this is usually set to 1 to get an add_mul operation
     404        # when computing the coefficient matrix that will be rolled over
     405        # we won't care about the existing data, so this is set to zero
     406        cdef Py_ssize_t res_coeff_ind
     407
     408        # FIXME: sig_on()
     409        for i in range(self._degree):
     410            for j in range(self._degree):
     411                print "i: {i}, j: {j}".format(i=i, j=j)
     412                if i+j >= self._degree:
     413                    res_matrix = tmp_matrix
     414                    res_coeff_ind = 0
     415                else:
     416                    res_matrix = res._coeffs[i+j]
     417                    res_coeff_ind = 1
     418               
     419                Mod_fgemm(self._linbox_field[0],
     420                        fflas_no_trans, fflas_no_trans,
     421                        self._nrows, self._ncols, self._nrows,
     422                        self._linbox_consts[1],
     423                        <ModFieldElement*>self._coeffs[i], self._nrows,
     424                        <ModFieldElement*>other._coeffs[j], self._nrows,
     425                        self._linbox_consts[res_coeff_ind],
     426                        <ModFieldElement*>res_matrix,
     427                        self._nrows)
     428
     429                print "after fgemm"
     430
     431                # FIXME: rollback code is wrong!
     432                if i+j >= self._degree:
     433                    #print "tmp_matrix:"
     434                    #for l in range(self._nrows):
     435                    #    s= ' '.join([str(tmp_matrix[l*self._nrows + m])
     436                    #            for m in range(self._ncols)])
     437                    #    print s
     438
     439                    # for each nonzero coefficient in the minimal polynomial
     440                    # of the base field add the result matrix * coefficient
     441                    # to the corresponding coefficient matrix
     442                    for k in range(self._degree):
     443                        print "k:",k
     444                        if self._linbox_field[0].isZero(\
     445                                self._minpoly_vector[k]):
     446                            continue
     447                        print "not zero"
     448                        print
     449
     450                        Mod_faxpy(self._linbox_field[0],
     451                                self._nrows*self._ncols,
     452                                self._minpoly_vector[k],
     453                                <ModFieldElement*>tmp_matrix, 1,
     454                                <ModFieldElement*>\
     455                                        res._coeffs[(i+j+k) % self._degree], 1)
     456                        print "after faxpy"
     457
     458                        #print "res coeffs 0"
     459                        #res._print_coeff_mat(0)
     460                        #print "res coeffs 1"
     461                        #res._print_coeff_mat(1)
     462                        #print "tmp_matrix:"
     463                        #for l in range(self._nrows):
     464                        #    s= ' '.join([str(tmp_matrix[l*self._nrows + m])
     465                        #            for m in range(self._ncols)])
     466                        #    print s
     467
     468                #print "res coeffs 0"
     469                #res._print_coeff_mat(0)
     470                #print "res coeffs 1"
     471                #res._print_coeff_mat(1)
     472
     473        # FIXME: sig off
     474        sage_free(tmp_matrix)
     475
     476        return res
     477
     478
  • new file sage/matrix/matrix_dense_slice_template_header.pxi

    diff --git a/sage/matrix/matrix_dense_slice_template_header.pxi b/sage/matrix/matrix_dense_slice_template_header.pxi
    new file mode 100644
    - +  
     1"""
     2Dense Matrix Template for C/C++ Library Interfaces
     3"""
     4
     5cimport matrix_dense
     6
     7cdef class Matrix_dense_slice_template(matrix_dense.Matrix_dense):
     8    cdef celement* _coeffs_block
     9    cdef celement** _coeffs
     10    cdef unsigned long _char
     11    cdef unsigned long _degree
     12
     13cdef class Matrix_dense_slice_template_GFq(Matrix_dense_slice_template):
     14    cdef ModField* _linbox_field
     15    cdef ModFieldElement* _minpoly_vector
     16    cdef ModFieldElement* _linbox_consts
  • new file sage/matrix/matrix_modq_dense_double.pxd

    diff --git a/sage/matrix/matrix_modq_dense_double.pxd b/sage/matrix/matrix_modq_dense_double.pxd
    new file mode 100644
    - +  
     1include "../ext/cdefs.pxi"
     2
     3ctypedef double celement
     4
     5from sage.libs.linbox.modular cimport ModDoubleField as ModField, ModDoubleFieldElement as ModFieldElement
     6
     7include "matrix_dense_slice_template_header.pxi"
     8
     9cdef class Matrix_modq_dense_double(Matrix_dense_slice_template_GFq):
     10    pass
  • new file sage/matrix/matrix_modq_dense_double.pyx

    diff --git a/sage/matrix/matrix_modq_dense_double.pyx b/sage/matrix/matrix_modq_dense_double.pyx
    new file mode 100644
    - +  
     1###############################################################################
     2#   Sage: Open Source Mathematical Software
     3#       Copyright (C) 2011 Burcin Erocal <burcin@erocal.org>
     4#  Distributed under the terms of the GNU General Public License (GPL),
     5#  version 2 or any later version.  The full text of the GPL is available at:
     6#                  http://www.gnu.org/licenses/
     7###############################################################################
     8
     9include "../ext/stdsage.pxi"
     10include "../ext/interrupt.pxi"
     11
     12from sage.structure.element cimport Element, FieldElement
     13from sage.rings.finite_rings.element_givaro cimport FiniteField_givaroElement
     14from sage.rings.finite_rings.element_ext_pari import FiniteField_ext_pariElement
     15
     16
     17from sage.libs.linbox.fflas cimport ModDouble_fgemm as Mod_fgemm,\
     18        ModDouble_faxpy as Mod_faxpy,\
     19        ModDouble_fgemv as Mod_fgemv,\
     20        ModDouble_invert_in_place as Mod_invert_in_place
     21from sage.libs.linbox.fflas cimport fflas_trans_enum, fflas_no_trans, fflas_trans
     22
     23cdef void sage_element_to_celement_vector(celement* dest, Element el,
     24        unsigned long characteristic, unsigned int degree) except +:
     25    """
     26    Converts a given Sage element to a vector of ctypes.
     27    """
     28    cdef Py_ssize_t ind
     29    if isinstance(el, FiniteField_givaroElement):
     30        quo = el.log_to_int()
     31
     32        for ind in range(degree):
     33            quo, coeff = quo.quo_rem(characteristic)
     34            dest[ind] = <celement>coeff
     35    elif isinstance(el, FiniteField_ext_pariElement):
     36        v = el.vector()
     37        for ind in range(degree):
     38            dest[ind] = <celement>v[ind]
     39    else:
     40        raise NotImplementedError("cannot handle elements of class {cls}".format(cls=type(el)))
     41
     42include "matrix_dense_slice_template.pxi"
     43
     44
     45cdef class Matrix_modq_dense_double(Matrix_dense_slice_template_GFq):
     46    """
     47    TESTS::
     48
     49        sage: from sage.matrix.matrix_modq_dense_float import Matrix_modq_dense_float
     50        sage: K.<a> = GF(5^2)
     51        sage: MS = MatrixSpace(K, 3, 3)
     52        sage: M = Matrix_modq_dense_float(MS, range(9))
     53        sage: TestSuite(M).run() is None
     54        True
     55    """
     56    pass
  • new file sage/matrix/matrix_modq_dense_float.pxd

    diff --git a/sage/matrix/matrix_modq_dense_float.pxd b/sage/matrix/matrix_modq_dense_float.pxd
    new file mode 100644
    - +  
     1include "../ext/cdefs.pxi"
     2
     3ctypedef float celement
     4
     5from sage.libs.linbox.modular cimport ModFloatField as ModField, ModFloatFieldElement as ModFieldElement
     6
     7include "matrix_dense_slice_template_header.pxi"
     8
     9cdef class Matrix_modq_dense_float(Matrix_dense_slice_template_GFq):
     10    pass
  • new file sage/matrix/matrix_modq_dense_float.pyx

    diff --git a/sage/matrix/matrix_modq_dense_float.pyx b/sage/matrix/matrix_modq_dense_float.pyx
    new file mode 100644
    - +  
     1###############################################################################
     2#   Sage: Open Source Mathematical Software
     3#       Copyright (C) 2011 Burcin Erocal <burcin@erocal.org>
     4#  Distributed under the terms of the GNU General Public License (GPL),
     5#  version 2 or any later version.  The full text of the GPL is available at:
     6#                  http://www.gnu.org/licenses/
     7###############################################################################
     8
     9include "../ext/stdsage.pxi"
     10include "../ext/interrupt.pxi"
     11
     12from sage.structure.element cimport Element, FieldElement
     13from sage.rings.finite_rings.element_givaro cimport FiniteField_givaroElement
     14from sage.rings.finite_rings.element_ext_pari import FiniteField_ext_pariElement
     15
     16
     17from sage.libs.linbox.fflas cimport ModFloat_fgemm as Mod_fgemm,\
     18        ModFloat_faxpy as Mod_faxpy,\
     19        ModFloat_fgemv as Mod_fgemv,\
     20        ModFloat_invert_in_place as Mod_invert_in_place
     21from sage.libs.linbox.fflas cimport fflas_trans_enum, fflas_no_trans, fflas_trans
     22
     23cdef void sage_element_to_celement_vector(celement* dest, Element el,
     24        unsigned long characteristic, unsigned int degree) except +:
     25    """
     26    Converts a given Sage element to a vector of ctypes.
     27    """
     28    cdef Py_ssize_t ind
     29    if isinstance(el, FiniteField_givaroElement):
     30        quo = el.log_to_int()
     31
     32        for ind in range(degree):
     33            quo, coeff = quo.quo_rem(characteristic)
     34            dest[ind] = <celement>coeff
     35    elif isinstance(el, FiniteField_ext_pariElement):
     36        v = el.vector()
     37        for ind in range(degree):
     38            dest[ind] = <celement>v[ind]
     39    else:
     40        raise NotImplementedError("cannot handle elements of class {cls}".format(cls=type(el)))
     41
     42include "matrix_dense_slice_template.pxi"
     43
     44
     45cdef class Matrix_modq_dense_float(Matrix_dense_slice_template_GFq):
     46    """
     47    TESTS::
     48
     49        sage: from sage.matrix.matrix_modq_dense_float import Matrix_modq_dense_float
     50        sage: K.<a> = GF(5^2)
     51        sage: MS = MatrixSpace(K, 3, 3)
     52        sage: M = Matrix_modq_dense_float(MS, range(9))
     53        sage: TestSuite(M).run() is None
     54        True
     55    """
     56    pass