Changeset 13196


Ignore:
Timestamp:
08/30/12 12:54:28 (13 years ago)
Author:
utke
Message:

CHG: streamline the allocation of the temporary vector for the solution (single owner allocates & deletes)

Location:
issm/trunk-jpl/src/c/modules/Solverx
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/Solverx/Solverx.h

    r13195 r13196  
    2424
    2525void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf,Parameters* parameters);
    26 void SolverxSeq(IssmPDouble** pX,IssmPDouble* A,IssmPDouble* B,int n);
     26void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
    2727
    2828#ifdef _HAVE_ADOLC_
    29 void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n, Parameters* parameters);
     29void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
    3030ADOLC_ext_fct EDF_for_solverx;
    3131#endif

  • issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp

    r13195 r13196  
    2929        int M,N,N2,s;
    3030        SeqVec *uf = NULL;
    31         IssmDouble *x  = NULL;
    3231
    3332        Kff->GetSize(&M,&N);
    … …  
    3635        if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
    3736        if(M!=N)_error_("Stiffness matrix should be square!");
    38 
     37        IssmDouble *x  = xNew<IssmDouble>(N);
    3938#ifdef _HAVE_ADOLC_
    40         SolverxSeq(&x,Kff->matrix,pf->vector,N,parameters);
     39        SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
    4140#else
    42         SolverxSeq(&x,Kff->matrix,pf->vector,N);
     41        SolverxSeq(x,Kff->matrix,pf->vector,N);
    4342#endif
    44         uf=new SeqVec(x,N);
     43        uf=new SeqVec(x,N);     
     44        xDelete(x);
    4545
    4646        /*Assign output pointers:*/
    … …  
    5656int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y) {
    5757    if(m*(m+1)!=n)_error_("Stiffness matrix should be square!");
    58     SolverxSeq(&y,x, x+m*m, m);
     58    SolverxSeq(y,x, x+m*m, m);
    5959    return 0;
    6060}
    6161
    62 void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n, Parameters* parameters){//{{{
     62void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){//{{{
    6363        // pack inputs to conform to the EDF-prescribed interface
    64         IssmDouble*  adoubleEDF_X=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
    65         for(int i=0; i<n*n;i++)adoubleEDF_X[i]    =A[i]; // pack matrix
    66         for(int i=0; i<n;  i++)adoubleEDF_X[i+n*n]=B[i]; // pack the right hand side
    67         IssmPDouble* pdoubleEDF_X=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
    68         IssmPDouble* pdoubleEDF_Y=xNew<IssmPDouble>(n);       // provide space to transfer outputs during call_ext_fct
     64        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
     65        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i]; // pack matrix
     66        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
     67        IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
     68        IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);       // provide space to transfer outputs during call_ext_fct
    6969        // call the wrapped solver through the registry entry we retrieve from parameters
    7070        call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
    71                      n*(n+1), pdoubleEDF_X, adoubleEDF_X,
    72                      n, pdoubleEDF_Y, B);
    73         xDelete(adoubleEDF_X);
    74         xDelete(pdoubleEDF_X);
    75         xDelete(pdoubleEDF_Y);
     71                     n*(n+1), pdoubleEDFin, adoubleEDFin,
     72                     n, pdoubleEDFout,X);
     73        xDelete(adoubleEDFin);
     74        xDelete(pdoubleEDFin);
     75        xDelete(pdoubleEDFout);
    7676}
    7777/*}}}*/
    7878#endif
    79 void SolverxSeq(IssmPDouble** pX,IssmPDouble* A,IssmPDouble* B,int n){ //{{{
     79void SolverxSeq(IssmPDouble * X, IssmPDouble * A, IssmPDouble * B,int n){ //{{{
    8080        #ifdef _HAVE_GSL_
    8181        /*GSL Matrices and vectors: */
    … …  
    8787        /*A will be modified by LU decomposition. Use copy*/
    8888        double* Acopy = xNew<double>(n*n);
    89         xMemCpy<double>(Acopy,A,n*n);
     89        xMemCpy(Acopy,A,n*n);
    9090
    9191        /*Initialize gsl matrices and vectors: */
    … …  
    103103
    104104        /*Copy result*/
    105         double* X = xNew<double>(n);
    106         memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
     105        xMemCpy(X,gsl_vector_ptr(x,0),n);
    107106
    108107        /*Clean up and assign output pointer*/
    109         xDelete<double>(Acopy);
     108        xDelete(Acopy);
    110109        gsl_permutation_free(p);
    111110        gsl_vector_free(x);
    112         *pX=X;
    113111        #endif
    114112}
Note: See TracChangeset for help on using the changeset viewer.