Changeset 1230 in lmdz_wrf for trunk/tools


Ignore:
Timestamp:
Oct 25, 2016, 7:55:03 PM (8 years ago)
Author:
lfita
Message:

Adding `operations': Function to perform different operations to a pair of matrices of values

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/generic_tools.py

    r1229 r1230  
    7777# num_split: Function to split a string at each numeric value keeping the number as the last character of each cut
    7878# oper_submatrix: Function to perform an operation of a given matrix along a sub-set of values along its dimensions
     79# operations: Function to perform different operations to a pair of matrices of values
    7980# period_information_360d: Function to provide the information of a given period idate, edate (dates in
    8081#  [YYYY][MM][DD][HH][MI][SS] format) in a 360 years calendar
    … …  
    94479448            if searchInlist(firstchars,noc): firstchars.remove(noc)
    94489449
    9449     newtext = text + ''
     9450    newtext = text
    94509451    for char in firstchars:
     9452#        newtext = newtext.replace(char.decode('iso-8859-15'),chars[char].decode('iso-8859-15'))
    94519453        newtext = newtext.replace(char,chars[char])
    94529454        chars.pop(char)
    … …  
    94889490    newln = newln.replace('ï', '\\"i')
    94899491    newln = newln.replace('ö', '\\"o')
     9492    newln = newln.replace('\xf6', '\\"o')
    94909493    newln = newln.replace('ÃŒ', '\\"u')
    94919494
    … …  
    98219824    timestep.sort()
    98229825    return timestep
     9826
     9827def operations(opN,matA,matB=None):
     9828    """ Function to perform different operations to a pair of matrices of values
     9829      opN: name of the operation
     9830        'add': adding matB to matA (matA + matB)
     9831        'centerderiv',[N],[ord],[dim]: un-scaled center [N]-derivative of order [ord] along dimension [dim] of matA
     9832        'div': dividing by matB (matA / matB)
     9833        'forwrdderiv',[N],[ord],[dim]: un-scaled forward [N]-derivative of order [ord] along dimension [dim] of matA
     9834        'inv': inverting matA (1/matA)
     9835        'mul': multiplying ny matB (matA * matB)
     9836        'pot': powering with matB (matA ** matB)
     9837        'sub': substracting matB to matA (matA - matB)
     9838      matA: initial values
     9839      matB: values to operate with
     9840    >>> operations('centerderiv,1,8,0',np.arange(81.).reshape(9,9))
     9841    [[ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9842     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9843     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9844     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9845     [ 9.  9.  9.  9.  9.  9.  9.  9.  9.]
     9846     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9847     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9848     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9849     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]]
     9850    >>> operations('forwrdderiv,1,6,0',np.arange(81.).reshape(9,9))
     9851    [[ 9.  9.  9.  9.  9.  9.  9.  9.  9.]
     9852     [ 9.  9.  9.  9.  9.  9.  9.  9.  9.]
     9853     [ 9.  9.  9.  9.  9.  9.  9.  9.  9.]
     9854     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9855     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9856     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9857     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9858     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9859     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]]
     9860    >>> operations('forwrdderiv,2,2,0',np.arange(81.).reshape(9,9))
     9861    [[ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9862     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9863     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9864     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9865     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9866     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9867     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9868     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]
     9869     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]]
     9870    """
     9871    fname = 'operations'
     9872    newmat = None
     9873
     9874    matAshape = list(matA.shape)
     9875    print fname + '; Lluis matAshape:', matAshape
     9876
     9877    # From: https://en.wikipedia.org/wiki/Finite_difference_coefficient
     9878    # Coefficients for center derivatives at different orders (as key)
     9879    centerderiv = np.zeros((6,8,9), dtype=np.float)
     9880    # Different available orders for each derivative
     9881    centerorder = {}
     9882    centerorder[0] = [1,3,5,7]
     9883    centerorder[1] = [1,3,5,7]
     9884    centerorder[2] = [1,3,5]
     9885    centerorder[3] = [1,3,5]
     9886    centerorder[4] = [1]
     9887    centerorder[5] = [1]
     9888
     9889    # First derivative
     9890    centerderiv[0,1,:]= [        0.,        0.,        0.,    -1./2.,        0.,     1./2.,        0.,        0.,        0.]
     9891    centerderiv[0,3,:]= [        0.,        0.,    1./12.,    -2./3.,        0.,     2./3.,   -1./12.,        0.,        0.]
     9892    centerderiv[0,5,:]= [        0.,   -1./60.,    3./20.,    -3./4.,        0.,     3./4.,   -3./20.,     1./60,        0.]
     9893    centerderiv[0,7,:]= [   1./280.,  -4./105.,     1./5.,    -4./5.,        0.,     4./5.,    -1./5.,   4./105.,  -1./280.]
     9894    # Second derivative
     9895    centerderiv[1,1,:]= [        0.,        0.,       0.,         1.,       -2.,        1.,        0.,        0.,        0.]
     9896    centerderiv[1,3,:]= [        0.,        0.,   -1./12.,     4./3.,    -5./2.,     4./3.,   -1./12.,        0.,        0.]
     9897    centerderiv[1,5,:]= [        0.,    1./90.,   -3./20.,     3./2.,  -49./18.,     3./2.,   -.3/20.,    1./90.,        0.]
     9898    centerderiv[1,7,:]= [  -1./560.,   8./315.,    -1./5.,     8./5., -205./72.,     8./5.,    -1./5.,   8./315.,  -1./560.]
     9899    # 3rd derivative
     9900    centerderiv[2,1,:]= [        0.,        0.,    -1./2.,        1.,        0.,       -1.,     1./2.,        0.,        0.]
     9901    centerderiv[2,3,:]= [        0.,     1./8.,       -1.,    13./8.,        0.,   -13./8.,        1.,    -1./8.,        0.]
     9902    centerderiv[2,5,:]= [  -7./240.,    3./10.,-169./120.,   61./30.,        0.,  -61./30., 169./120.,   -3./10.,   7./240.]
     9903    # 4th derivative
     9904    centerderiv[3,1,:]= [        0.,        0.,        1.,       -4.,        6.,       -4.,        1.,        0.,        0.]
     9905    centerderiv[3,3,:]= [        0.,    -1./6.,        2.,   -13./2.,    28./3.,   -13./2.,        2.,    -1./6.,        0.]
     9906    centerderiv[3,5,:]= [   7./240.,    -2./5.,  169./60., -122./15.,    91./8., -122./15.,  169./60.,    -2./5.,   7./240.]
     9907    # 5th derivative
     9908    centerderiv[4,1,:]= [        0.,    -1./2.,        2.,    -5./2.,        0.,     5./2.,       -2.,     1./2.,        0.]
     9909    # 6th derivative
     9910    centerderiv[5,1,:]= [        0.,        1.,       -6.,       15.,      -20.,       15.,       -6.,        1.,        0.]
     9911
     9912    # Coefficients for forward derivatives at different orders (as key)
     9913    forwrdderiv = np.zeros((6,6,9), dtype=np.float)
     9914    # Different available orders for each derivative
     9915    forwrdorder = {}
     9916    forwrdorder[0] = range(6)
     9917    forwrdorder[1] = range(6)
     9918    forwrdorder[2] = range(6)
     9919    forwrdorder[3] = range(4)
     9920
     9921    # First derivative
     9922    forwrdderiv[0,0,:]= [       -1.,        1.,        0.,        0.,        0.,        0.,        0.,        0.,        0.]
     9923    forwrdderiv[0,1,:]= [    -3./2.,        2.,    -1./2.,        0.,        0.,        0.,        0.,        0.,        0.]
     9924    forwrdderiv[0,2,:]= [   -11./6.,        3.,    -3./2.,     1./3.,        0.,        0.,        0.,        0.,        0.]
     9925    forwrdderiv[0,3,:]= [  -25./12.,        4.,       -3.,     4./3.,    -1./4.,        0.,        0.,        0.,        0.]
     9926    forwrdderiv[0,4,:]= [ -137./60.,        5.,       -5.,    10./3.,    -5./4.,     1./5.,        0.,        0.,        0.]
     9927    forwrdderiv[0,5,:]= [  -49./20.,        6.,   -15./2.,    20./3.,   -15./4.,     6./5.,    -1./6.,        0.,        0.]
     9928    # Second derivative
     9929    forwrdderiv[1,0,:]= [        1.,       -2.,        1.,        0.,        0.,        0.,        0.,        0.,        0.]
     9930    forwrdderiv[1,1,:]= [        2.,       -5.,        4.,       -1.,        0.,        0.,        0.,        0.,        0.]
     9931    forwrdderiv[1,2,:]= [   35./12.,   -26./3.,    19./2.,   -14./3.,   11./12.,        0.,        0.,        0.,        0.]
     9932    forwrdderiv[1,3,:]= [    15./4.,   -77./6.,   107./6.,      -13.,   61./12.,    -5./6.,        0.,        0.,        0.]
     9933    forwrdderiv[1,4,:]= [  203./45.,   -87./5.,   117./4.,  -254./9.,    33./2.,   -27./5., 137./180.,        0.,        0.]
     9934    forwrdderiv[1,5,:]= [  469./90., -223./10.,  879./20., -949./18.,       41., -201./10.,1019./180.,   -7./10.,        0.]
     9935    # 3rd derivative
     9936    forwrdderiv[2,0,:]= [       -1.,        3.,       -3.,        1.,        0.,        0.,        0.,        0.,        0.]
     9937    forwrdderiv[2,1,:]= [    -5./2.,        9.,      -12.,        7.,    -3./2.,        0.,        0.,        0.,        0.]
     9938    forwrdderiv[2,2,:]= [   -17./4.,    71./4.,   -59./2.,    49./2.,   -41./4.,     7./4.,        0.,        0.,        0.]
     9939    forwrdderiv[2,3,:]= [   -49./8.,       29.,  -461./8.,       62.,  -307./8.,       13.,   -15./8.,        0.,        0.]
     9940    forwrdderiv[2,4,:]= [-967./120.,  638./15.,-3929./40.,   389./3.,-2545./24.,   268./5.,-1849./120,   29./15.,        0.]
     9941    forwrdderiv[2,5,:]= [ -801./80.,   349./6.,-18353./120,2391./10., -1457./6., 4891./30.,  -561./8.,  527./30.,-469./240.]
     9942    # 4th derivative
     9943    forwrdderiv[3,0,:]= [        1.,       -4.,        6.,       -4.,        1.,        0.,        0.,        0.,        0.]
     9944    forwrdderiv[3,1,:]= [        3.,      -14.,       26.,      -24.,       11.,       -2.,        0.,        0.,        0.]
     9945    forwrdderiv[3,2,:]= [    35./6.,      -31.,   137./2.,  -242./3.,   107./2.,      -19.,    17./6.,        0.,        0.]
     9946    forwrdderiv[3,3,:]= [    28./3.,  -111./2.,      142., -1219./6.,      176.,  -185./2.,    82./3.,    -7./2.,        0.]
     9947    forwrdderiv[3,4,:]= [ 1069./80.,-1316./15.,15289./60., -2144./5.,10993./24.,-4772./15., 2803./20., -536./15., 967./240.]
     9948
     9949    # Dictionary with the description of the operations
     9950    Lopn = {'add': '+', 'div': '/',  'inv': '^(-1)', 'mul': '*', 'pot': '^',         \
     9951      'sub': '-'}
     9952
     9953    opavail = ['add', 'centerderiv', 'div', 'forwrdderiv', 'inv', 'mul', 'pot', 'sub']
     9954
     9955    if opN[0:11] == 'centerderiv':
     9956        opn = 'centerderiv'
     9957        Nderiv = int(opN.split(',')[1])
     9958        order = int(opN.split(',')[2])
     9959        dim = int(opN.split(',')[3])
     9960        Lopn[opn]= str(Nderiv)+'-center_deriv_ord('+str(order)+ ')[' + str(dim) + ']'
     9961        if not centerorder.has_key(Nderiv-1):
     9962            print errormsg
     9963            print '  ' + fname + ': ', Nderiv,'-center derivative does not exist!!'
     9964            print '  existing ones:', list(centeroder.keys())+1
     9965            quit()
     9966        derivorders = centerorder[Nderiv-1]
     9967        if not searchInlist(derivorders,order-1):
     9968            print errormsg
     9969            print '  ' + fname + ': ', Nderiv,'-center derivative does not have ' +  \
     9970              'order:',  order,'!!'
     9971            print '  existing ones:', np.array(derivorders) + 1
     9972            quit()
     9973        if dim > len(matAshape)-1:
     9974            print errormsg
     9975            print '  '+fname+ ': wrong dimension:', dim, 'to compute the center-' +  \
     9976              'derivative!!'
     9977            print '  matrix has a smaller rank:', len(matAshape)
     9978            quit(-1)
     9979        derivcoeff = centerderiv[Nderiv-1,order-1,:]
     9980         
     9981    elif opN[0:11] == 'forwrdderiv':
     9982        opn = 'forwrdderiv'
     9983        Nderiv = int(opN.split(',')[1])
     9984        order = int(opN.split(',')[2])
     9985        dim = int(opN.split(',')[3])
     9986        Lopn[opn]= str(Nderiv)+'-forward_deriv_ord('+str(order)+ ')[' + str(dim) + ']'
     9987        if not forwrdorder.has_key(Nderiv-1):
     9988            print errormsg
     9989            print '  ' + fname + ': ', Nderiv,'-forward derivative does not exist!!'
     9990            print '  existing ones:', list(forwrdoder.keys())+1
     9991            quit()
     9992        derivorders = forwrdorder[Nderiv-1]
     9993        if not searchInlist(derivorders,order-1):
     9994            print errormsg
     9995            print '  ' + fname + ': ', Nderiv,'-forward derivative does not have ' + \
     9996              'order:',  order,'!!'
     9997            print '  existing ones:', np.array(derivorders) + 1
     9998            quit()
     9999        if dim > len(matAshape)-1:
     10000            print errormsg
     10001            print '  '+fname+ ': wrong dimension:', dim, 'to compute the forward-' + \
     10002              'derivative!!'
     10003            print '  matrix has a smaller rank:', len(matAshape)
     10004            quit(-1)
     10005        derivcoeff = forwrdderiv[Nderiv-1,order-1,:]
     10006    else:
     10007        opn = opN
     10008
     10009    valtype = mat.dtype
     10010
     10011    if opn == 'add':
     10012        newmat = matA + matB
     10013    elif opn == 'centerderiv':
     10014        sdim = matAshape[dim]
     10015        # Getting the different slices for each coefficient
     10016        # newmat = matA_4*coef_4 + matA_3*coef_3 + matA_2*coef_2 + matA_1*coef_1 +
     10017        #   matA*coef0 + matA1*coef1 + mat2*coef2 + matA3*coef3 + mat4*coef4
     10018        newmat = np.zeros(tuple(matAshape), dtype=valtype)
     10019        mat0 = []
     10020        mat1_ = []
     10021        mat1 = []
     10022        mat2_ = []
     10023        mat2 = []
     10024        mat3_ = []
     10025        mat3 = []
     10026        mat4_ = []
     10027        mat4 = []
     10028
     10029        Nterms = 2*int((Nderiv-1)/2)+order+1
     10030        for idim in range(len(matAshape)):
     10031            if idim == dim:
     10032                sd = (Nterms-1)/2
     10033                ed = sdim - (Nterms-1)/2
     10034                mat0.append(slice(sd,ed))
     10035                mat1_.append(slice(sd-1,ed-1))
     10036                mat1.append(slice(sd+1,ed+1))
     10037                mat2_.append(slice(sd-2,ed-2))
     10038                mat2.append(slice(sd+2,ed+2))
     10039                mat3_.append(slice(sd-3,ed-3))
     10040                mat3.append(slice(sd+3,ed+3))
     10041                mat4_.append(slice(sd-4,ed-4))
     10042                mat4.append(slice(sd+4,ed+4))
     10043            else:
     10044                mat0.append(slice(0,matAshape[idim]+1))
     10045                mat1_.append(slice(0,matAshape[idim]+1))
     10046                mat1.append(slice(0,matAshape[idim]+1))
     10047                mat2_.append(slice(0,matAshape[idim]+1))
     10048                mat2.append(slice(0,matAshape[idim]+1))
     10049                mat3_.append(slice(0,matAshape[idim]+1))
     10050                mat3.append(slice(0,matAshape[idim]+1))
     10051                mat4_.append(slice(0,matAshape[idim]+1))
     10052                mat4.append(slice(0,matAshape[idim]+1))
     10053        if Nterms == 3:
     10054            newmat[tuple(mat0)] = matA[tuple(mat1_)]*derivcoeff[3] +                 \
     10055              matA[tuple(mat0)]*derivcoeff[4] + matA[tuple(mat1)]*derivcoeff[5]
     10056        elif Nterms == 5:
     10057            newmat[tuple(mat0)] = matA[tuple(mat2_)]*derivcoeff[2] +                 \
     10058              matA[tuple(mat1_)]*derivcoeff[3] +                                     \
     10059              matA[tuple(mat0)]*derivcoeff[4] +                                      \
     10060              matA[tuple(mat1)]*derivcoeff[5] + matA[tuple(mat2)]*derivcoeff[6]
     10061        elif Nterms == 7:
     10062            newmat[tuple(mat0)] = matA[tuple(mat3_)]*derivcoeff[1] +                 \
     10063              matA[tuple(mat2_)]*derivcoeff[2] + matA[tuple(mat1_)]*derivcoeff[3] +  \
     10064              matA[tuple(mat0)]*derivcoeff[4] +                                      \
     10065              matA[tuple(mat1)]*derivcoeff[5] + matA[tuple(mat2)]*derivcoeff[6] +    \
     10066              matA[tuple(mat3)]*derivcoeff[7]
     10067        elif Nterms == 9:
     10068            newmat[tuple(mat0)] = matA[tuple(mat4_)]*derivcoeff[0] +                 \
     10069              matA[tuple(mat3_)]*derivcoeff[1] + matA[tuple(mat2_)]*derivcoeff[2] +  \
     10070              matA[tuple(mat1_)]*derivcoeff[3] +                                     \
     10071              matA[tuple(mat0)]*derivcoeff[4] +                                      \
     10072              matA[tuple(mat1)]*derivcoeff[5] + matA[tuple(mat2)]*derivcoeff[6] +    \
     10073              matA[tuple(mat3)]*derivcoeff[7] + matA[tuple(mat4)]*derivcoeff[8]
     10074    elif opn == 'div':
     10075        newmat = matA / matB
     10076    elif opn == 'forwrdderiv':
     10077        sdim = matAshape[dim]
     10078        # Getting the different slices for each coefficient
     10079        # newmat = matA0*coef0 + matA1*coef1 + matA2*coef2 + matA3*coef3 +
     10080        #   matA4*coef4 + matA5*coef5 + matA6*coef6 + matA7*coef7 + matA8*coef8
     10081        newmat = np.zeros(tuple(matAshape), dtype=valtype)
     10082        mat0 = []
     10083        mat1 = []
     10084        mat2 = []
     10085        mat3 = []
     10086        mat4 = []
     10087        mat5 = []
     10088        mat6 = []
     10089        mat7 = []
     10090        mat8 = []
     10091
     10092        Nterms = Nderiv + order
     10093        for idim in range(len(matAshape)):
     10094            if idim == dim:
     10095                sd = 0
     10096                ed = sdim-Nterms+1
     10097                mat0.append(slice(sd,ed))
     10098                mat1.append(slice(sd+1,ed+1))
     10099                mat2.append(slice(sd+2,ed+2))
     10100                mat3.append(slice(sd+3,ed+3))
     10101                mat4.append(slice(sd+4,ed+4))
     10102                mat5.append(slice(sd+5,ed+5))
     10103                mat6.append(slice(sd+6,ed+6))
     10104                mat7.append(slice(sd+7,ed+7))
     10105                mat8.append(slice(sd+8,ed+8))
     10106            else:
     10107                mat0.append(slice(0,matAshape[idim]+1))
     10108                mat1.append(slice(0,matAshape[idim]+1))
     10109                mat2.append(slice(0,matAshape[idim]+1))
     10110                mat3.append(slice(0,matAshape[idim]+1))
     10111                mat4.append(slice(0,matAshape[idim]+1))
     10112                mat5.append(slice(0,matAshape[idim]+1))
     10113                mat6.append(slice(0,matAshape[idim]+1))
     10114                mat7.append(slice(0,matAshape[idim]+1))
     10115                mat8.append(slice(0,matAshape[idim]+1))
     10116        if Nterms == 2:
     10117            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10118              matA[tuple(mat1)]*derivcoeff[1]
     10119        elif Nterms == 3:
     10120            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10121              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2]
     10122        elif Nterms == 4:
     10123            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10124              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2] +    \
     10125              matA[tuple(mat3)]*derivcoeff[3]
     10126        elif Nterms == 5:
     10127            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10128              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2] +    \
     10129              matA[tuple(mat3)]*derivcoeff[3] + matA[tuple(mat4)]*derivcoeff[4]
     10130        elif Nterms == 6:
     10131            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10132              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2] +    \
     10133              matA[tuple(mat3)]*derivcoeff[3] + matA[tuple(mat4)]*derivcoeff[4] +    \
     10134              matA[tuple(mat5)]*derivcoeff[5]
     10135        elif Nterms == 7:
     10136            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10137              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2] +    \
     10138              matA[tuple(mat3)]*derivcoeff[3] + matA[tuple(mat4)]*derivcoeff[4] +    \
     10139              matA[tuple(mat5)]*derivcoeff[5] + matA[tuple(mat6)]*derivcoeff[6]
     10140        elif Nterms == 8:
     10141            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10142              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2] +    \
     10143              matA[tuple(mat3)]*derivcoeff[3] + matA[tuple(mat4)]*derivcoeff[4] +    \
     10144              matA[tuple(mat5)]*derivcoeff[5] + matA[tuple(mat6)]*derivcoeff[6] +    \
     10145              matA[tuple(mat7)]*derivcoeff[7]
     10146        elif Nterms == 9:
     10147            newmat[tuple(mat0)] = matA[tuple(mat0)]*derivcoeff[0] +                  \
     10148              matA[tuple(mat1)]*derivcoeff[1] + matA[tuple(mat2)]*derivcoeff[2] +    \
     10149              matA[tuple(mat3)]*derivcoeff[3] + matA[tuple(mat4)]*derivcoeff[4] +    \
     10150              matA[tuple(mat5)]*derivcoeff[5] + matA[tuple(mat6)]*derivcoeff[6] +    \
     10151              matA[tuple(mat7)]*derivcoeff[7] + matA[tuple(mat8)]*derivcoeff[8]
     10152    elif opn == 'inv':
     10153        newmat = retype(1., valtype) / matA
     10154    elif opn == 'mul':
     10155        newmat = matA * matB
     10156    elif opn == 'pot':
     10157        newmat = matA ** matB
     10158    elif opn == 'sub':
     10159        newmat = matA - matB
     10160    else:
     10161        print '  ' + fname + ": operation '" + opn + "' does not exist !!"
     10162        print '  available ones:', opavail
     10163        quit(-1)
     10164
     10165    return newmat
     10166
    982310167#quit()
    982410168
Note: See TracChangeset for help on using the changeset viewer.