Changeset 20729


Ignore:
Timestamp:
06/14/16 08:53:20 (9 years ago)
Author:
Mathieu Morlighem
Message:

CHG: cosmetics

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

    r20728 r20729  
    2222
    2323        int         idima,idimb,idimc,idimd;
    24         IssmDouble dtemp_static[600];
     24        IssmDouble  dtemp_static[600];
     25        IssmDouble* dtemp_dynamic = NULL;
     26        IssmDouble* dtemp         = NULL;
    2527        _assert_(a && b && c && d);
    2628
    … …  
    5153                if (ncolc != idimc) _error_("Matrix B and C inner vectors not equal size.");
    5254                idimd=nrowc;
     55        }
     56
     57        /*Depending on the size of incoming matrices, we might need to use a dynamic allocation*/
     58        if(idima*idimc>600){
     59                IssmDouble *dtemp_dynamic=xNew<IssmDouble>(idima*idimc);
     60                dtemp = dtemp_dynamic;
     61        }
     62        else{
     63                dtemp = &dtemp_static[0];
    5364        }
    5465       
    … …  
    6172        /*  multiply (a*b)*c+d  */
    6273        if (idima*idimc*(idimb+idimd) <= idimb*idimd*(idima+idimc)) {
    63                 if(idima*idimc>600) {
    64                         /*Use dynamic allocation instead of dtemp_static*/
    65                         IssmDouble *dtemp;
    66                         dtemp=xNew<IssmDouble>(idima*idimc);
    67 
    68                         MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,dtemp,0);
    69                         MatrixMultiply(dtemp,idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
    70                         xDelete<IssmDouble>(dtemp);
    71                         return 1;
    72                 }
    73                 MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,&dtemp_static[0],0);
    74                 MatrixMultiply(&dtemp_static[0],idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
     74                MatrixMultiply(a,nrowa,ncola,itrna,b,nrowb,ncolb,itrnb,dtemp,0);
     75                MatrixMultiply(dtemp,idima,idimc,0,c,nrowc,ncolc,itrnc,d,iaddd);
    7576        }
    7677
    7778        /*  multiply a*(b*c)+d  */
    7879        else{
    79                 if(idimb*idimd>600) {
    80                         /*Use dynamic allocation instead of dtemp_static*/
    81                         IssmDouble *dtemp;
    82                         dtemp=xNew<IssmDouble>(idimb*idimd);
    83 
    84                         MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0);
    85                         MatrixMultiply(a,nrowa,ncola,itrna,dtemp,idimb,idimd,0,d,iaddd);
    86                         xDelete<IssmDouble>(dtemp);
    87                         return 1;
    88                 }
    89                 MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,&dtemp_static[0],0);
    90                 MatrixMultiply(a,nrowa,ncola,itrna,&dtemp_static[0],idimb,idimd,0,d,iaddd);
    91         }
    92 
     80                MatrixMultiply(b,nrowb,ncolb,itrnb,c,nrowc,ncolc,itrnc,dtemp,0);
     81                MatrixMultiply(a,nrowa,ncola,itrna,dtemp,idimb,idimd,0,d,iaddd);
     82        }
     83
     84        /*Cleanup and return*/
     85        xDelete<IssmDouble>(dtemp);
    9386        return 1;
    9487}/*}}}*/
    … …  
    337330
    338331        /*Intermediaries*/
    339         IssmDouble det;
    340         IssmDouble det_reciprocal;
     332        IssmDouble det,det_reciprocal;
    341333
    342334        /*Compute determinant*/
    … …  
    345337       
    346338        /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
    347         det_reciprocal = 1/det; 
     339        det_reciprocal = 1./det; 
    348340
    349341        /*Compute invert*/
    … …  
    458450
    459451        /*Intermediaries*/
    460         IssmDouble det;
    461         IssmDouble det_reciprocal;
     452        IssmDouble det,det_reciprocal;
    462453
    463454        /*Compute determinant*/
    … …  
    466457
    467458        /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
    468         det_reciprocal = 1/det; 
     459        det_reciprocal = 1./det; 
    469460
    470461        /*Compute invert*/
    … …  
    563554
    564555        /*Intermediaries*/
    565         IssmDouble det;
    566         IssmDouble det_reciprocal;
     556        IssmDouble det,det_reciprocal;
    567557
    568558        /*Compute determinant*/
    569559        Matrix4x4Determinant(&det,A);
    570         //if(fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
     560        if(fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
    571561
    572562        /*Multiplication is faster than division, so we multiply by the reciprocal*/
    573         det_reciprocal = 1/det;
     563        det_reciprocal = 1./det;
    574564
    575565        /*Compute adjoint matrix*/
Note: See TracChangeset for help on using the changeset viewer.