Changeset 14834


Ignore:
Timestamp:
05/01/13 15:56:42 (12 years ago)
Author:
Eric.Larour
Message:

CHG: done coding Assemble routines for vectors and matrices. Hook up to Solver not yet implemented. Will work on MUMPS solver next.
Assemble routine is not yet debugged. No segfault or mem-leak, but something is off in the results.

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/objects/Bucket.h

    r14833 r14834  
    187187                        int        *modes_forcpu       = NULL;
    188188
     189
    189190                        /*initialize buffers: */
    190191                        row_indices_forcpu=*prow_indices_forcpu;
    … …  
    216217                };
    217218                /*}}}*/
     219                void Marshall(int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu){ /*{{{*/
     220
     221                        /*intermediary: */
     222                        int         i;
     223                        int         j;
     224
     225                        /*buffers: */
     226                        int        *row_indices_forcpu = NULL;
     227                        doubletype *values_forcpu      = NULL;
     228                        int        *modes_forcpu       = NULL;
     229
     230
     231                        /*initialize buffers: */
     232                        row_indices_forcpu=*prow_indices_forcpu;
     233                        values_forcpu=*pvalues_forcpu;
     234                        modes_forcpu=*pmodes_forcpu;
     235
     236                        /*fill buffers with out values and indices and modes: */
     237                        for(i=0;i<m;i++){
     238                                row_indices_forcpu[i]=idxm[i];
     239                                values_forcpu[i]=values[i*n+j];
     240                                modes_forcpu[i]=mode;
     241                        }
     242
     243                        /*increment buffer for next Bucket who will marshall his data: */
     244                        row_indices_forcpu+=m;
     245                        values_forcpu+=m;
     246                        modes_forcpu+=m;
     247
     248                        /*output modified buffers: */
     249                        *prow_indices_forcpu=row_indices_forcpu;
     250                        *pvalues_forcpu=values_forcpu;
     251                        *pmodes_forcpu=modes_forcpu;
     252                };
     253                /*}}}*/
    218254                int MarshallSize(void){ /*{{{*/
    219255
    220                         if(type=MATRIX_BUCKET){
     256                        if(type==MATRIX_BUCKET){
    221257                                return m*n;
    222258                        }

  • issm/trunk-jpl/src/c/toolkits/issm/IssmAbsMat.h

    r14822 r14834  
    3434*/
    3535
     36template <class doubletype> class IssmAbsVec;
     37class Parameters;
     38
    3639template <class doubletype>
    3740class IssmAbsMat{
    … …  
    5356                virtual void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode)=0;
    5457                virtual void Convert(MatrixType type)=0;
     58                virtual IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pf, Parameters* parameters)=0;
    5559};
    5660

  • issm/trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h

    r14671 r14834  
    4444
    4545                /*IssmDenseMat constructors, destructors*/
    46                 /*FUNCTION IssmDenseMat(){{{*/
     46                /*IssmDenseMat(){{{*/
    4747                IssmDenseMat(){
    4848
    … …  
    5252                }
    5353                /*}}}*/
    54                 /*FUNCTION IssmDenseMat(int M,int N){{{*/
     54                /*IssmDenseMat(int M,int N){{{*/
    5555                IssmDenseMat(int pM,int pN){
    5656
    … …  
    6161                }
    6262                /*}}}*/
    63                 /*FUNCTION IssmDenseMat(int M,int N, doubletype sparsity){{{*/
     63                /*IssmDenseMat(int M,int N, doubletype sparsity){{{*/
    6464                IssmDenseMat(int pM,int pN, doubletype sparsity){
    6565
    … …  
    7070                }
    7171                /*}}}*/
    72                 /*FUNCTION IssmDenseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
     72                /*IssmDenseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
    7373                IssmDenseMat(int m,int n,int pM,int pN,int* d_nnz,int* o_nnz){
    7474
    … …  
    7979                }
    8080                /*}}}*/
    81                 /*FUNCTION IssmDenseMat(doubletype* serial_mat,int M,int N,doubletype sparsity){{{*/
     81                /*IssmDenseMat(doubletype* serial_mat,int M,int N,doubletype sparsity){{{*/
    8282                IssmDenseMat(doubletype* serial_mat,int pM,int pN,doubletype sparsity){
    8383
    … …  
    9292                }
    9393                /*}}}*/
    94                 /*FUNCTION IssmDenseMat(int M,int N, int connectivity, int numberofdofspernode){{{*/
     94                /*IssmDenseMat(int M,int N, int connectivity, int numberofdofspernode){{{*/
    9595                IssmDenseMat(int pM,int pN, int connectivity,int numberofdofspernode){
    9696
    … …  
    101101                }
    102102                /*}}}*/
    103                 /*FUNCTION ~IssmDenseMat(){{{*/
     103                /*~IssmDenseMat(){{{*/
    104104                ~IssmDenseMat(){
    105105
    … …  
    110110                /*}}}*/
    111111
    112                 /*IssmDenseMat specific routines */
    113                 /*FUNCTION Echo{{{*/
     112                /*IssmAbsMat virtual functions*/
     113                /*Echo{{{*/
    114114                void Echo(void){
    115115
    … …  
    124124                }
    125125                /*}}}*/
    126                 /*FUNCTION Assemble{{{*/
     126                /*Assemble{{{*/
    127127                void Assemble(void){
    128128
    … …  
    131131                }
    132132                /*}}}*/
    133                 /*FUNCTION Norm{{{*/
     133                /*Norm{{{*/
    134134                doubletype Norm(NormMode mode){
    135135
    … …  
    156156                }
    157157                /*}}}*/
    158                 /*FUNCTION GetSize{{{*/
     158                /*GetSize{{{*/
    159159                void GetSize(int* pM,int* pN){
    160160
    … …  
    164164                }
    165165                /*}}}*/
    166                 /*FUNCTION GetLocalSize{{{*/
     166                /*GetLocalSize{{{*/
    167167                void GetLocalSize(int* pM,int* pN){
    168168
    … …  
    172172                }
    173173                /*}}}*/
    174                 /*FUNCTION MatMult{{{*/
     174                /*MatMult{{{*/
    175175                void MatMult(IssmAbsVec<doubletype>* Xin,IssmAbsVec<doubletype>* AXin){
    176176
    … …  
    203203                }
    204204                /*}}}*/
    205                 /*FUNCTION Duplicate{{{*/
     205                /*Duplicate{{{*/
    206206                IssmDenseMat<doubletype>* Duplicate(void){
    207207
    … …  
    212212                }
    213213                /*}}}*/
    214                 /*FUNCTION ToSerial{{{*/
     214                /*ToSerial{{{*/
    215215                doubletype* ToSerial(void){
    216216
    … …  
    225225                }
    226226                /*}}}*/
    227                 /*FUNCTION SetValues{{{*/
     227                /*SetValues{{{*/
    228228                void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode){
    229229
    … …  
    243243                }
    244244                /*}}}*/
    245                 /*FUNCTION Convert{{{*/
     245                /*Convert{{{*/
    246246                void Convert(MatrixType type){
    247247
    … …  
    250250                }
    251251                /*}}}*/         
     252                /*Solve{{{*/
     253                IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pfin, Parameters* parameters){
     254
     255                        /*First off, we assume that the type of IssmAbsVec is IssmSeqVec. So downcast: */
     256                        IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin;
     257
     258                        #ifdef _HAVE_GSL_
     259                                /*Intermediary: */
     260                                int M,N,N2;
     261                                IssmSeqVec<IssmDouble> *uf = NULL;
     262
     263                                Kff->GetSize(&M,&N);
     264                                pf->GetSize(&N2);
     265
     266                                if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
     267                                if(M!=N)_error_("Stiffness matrix should be square!");
     268                                #ifdef _HAVE_ADOLC_
     269                                        ensureContiguousLocations(N);
     270                                #endif
     271                                IssmDouble *x  = xNew<IssmDouble>(N);
     272                                #ifdef _HAVE_ADOLC_
     273                                        SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
     274                                #els
     275                                        SolverxSeq(x,Kff->matrix,pf->vector,N);
     276                                #endif
     277
     278                                uf=new IssmSeqVec<IssmDouble>(x,N);
     279                                xDelete(x);
     280
     281                                /*Assign output pointers:*/
     282                                IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();
     283                                out->vector=(IssmAbsVec<IssmDouble>*)uf;
     284                                *pout=out;
     285
     286                        #else
     287                                _error_("GSL support not compiled in!");
     288                        #endif
     289
     290                }/*}}}*/
    252291
    253292};

  • issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h

    r14822 r14834  
    3434template <class doubletype> class IssmDenseMat;
    3535template <class doubletype> class IssmMpiDenseMat;
     36class Parameters;
    3637
    3738template <class doubletype>
    … …  
    199200                        matrix->convert(type);
    200201                }/*}}}*/
     202                IssmVec<doubletype>* Solve(IssmVec<doubletype>* pf, Parameters* parameters){ /*{{{*/
     203                       
     204                        IssmVec<doubletype>* outvector=NULL;
     205
     206                        outvector=new IssmVec<doubletype>();
     207
     208                        outvector->vector=this->matrix->Solve(pf->vector,parameters);
     209
     210                }/*}}}*/
    201211
    202212};

  • issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h

    r14833 r14834  
    142142                                                        printf("%g ",this->matrix[j*this->N+k]);
    143143                                                }
     144                                                printf("\n");
    144145                                        }
    145146                                }
    … …  
    302303                        RowRank=DetermineRowRankFromLocalSize(M,m,comm);
    303304
    304                         /*Now, sort out our dataset of buckets according to cpu ownership of rows: */
     305                        /*Now, sort out our dataset of buckets according to cpu ownership of rows: {{{*/
    305306                        bucketsforcpu=xNew<DataSet*>(num_procs);
    306307
    … …  
    313314                                bucketsforcpu[i]=bucketsofcpu_i;
    314315                        }
    315 
    316                         /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this
     316                        /*}}}*/
     317
     318                        /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this  {{{
    317319                         * dataset owns correspond to rows that are owned by cpu i, not j!. Out of all the buckets we own, make row,col,value,insert_mode
    318320                         * vectors that will be shipped around the cluster: */
    319321                        this->BucketsBuildScatterBuffers(&numvalues_forcpu,&row_indices_forcpu,&col_indices_forcpu,&values_forcpu,&modes_forcpu,bucketsforcpu,num_procs);
    320 
    321                         /*Now, we need to allocate on each cpu arrays to receive data from all the other cpus. To know what we need to allocate, we need
     322                        /*}}}*/
     323
     324                        /*Now, we need to allocate on each cpu arrays to receive data from all the other cpus. To know what we need to allocate, we need  {{{
    322325                         *some scatter calls: */
    323326                        numvalues_fromcpu   = xNew<int>(num_procs);
    324327                        for(i=0;i<num_procs;i++){
    325                                 MPI_Scatter(numvalues_forcpu,num_procs,MPI_INT,numvalues_fromcpu+i,1,MPI_INT,i,comm);
    326                         }
    327                         for(i=0;i<num_procs;i++){
    328                                 row_indices_fromcpu[i]=xNew<int>(numvalues_fromcpu[i]);
    329                                 col_indices_fromcpu[i]=xNew<int>(numvalues_fromcpu[i]);
    330                                 values_fromcpu[i]=xNew<doubletype>(numvalues_fromcpu[i]);
    331                                 modes_fromcpu[i]=xNew<int>(numvalues_fromcpu[i]);
    332                         }
    333 
     328                                MPI_Scatter(numvalues_forcpu,1,MPI_INT,numvalues_fromcpu+i,1,MPI_INT,i,comm);
     329                        }
     330                       
     331                        row_indices_fromcpu=xNew<int*>(num_procs);
     332                        col_indices_fromcpu=xNew<int*>(num_procs);
     333                        values_fromcpu=xNew<doubletype*>(num_procs);
     334                        modes_fromcpu=xNew<int*>(num_procs);
     335                        for(i=0;i<num_procs;i++){
     336                                int size=numvalues_fromcpu[i];
     337                                if(size){
     338                                        row_indices_fromcpu[i]=xNew<int>(size);
     339                                        col_indices_fromcpu[i]=xNew<int>(size);
     340                                        values_fromcpu[i]=xNew<doubletype>(size);
     341                                        modes_fromcpu[i]=xNew<int>(size);
     342                                }
     343                                else{
     344                                        row_indices_fromcpu[i]=NULL;
     345                                        col_indices_fromcpu[i]=NULL;
     346                                        values_fromcpu[i]=NULL;
     347                                        modes_fromcpu[i]=NULL;
     348                                }
     349                        }
     350                        /*}}}*/
     351
     352                        /*Scatter values around: {{{*/
    334353                        /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given
    335354                         * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the MPI_Scatterv prototype:
    … …  
    344363                        }
    345364
    346                         /*Start the scattering: */
    347365                        for(i=0;i<num_procs;i++){
    348366                                MPI_Scatterv( row_indices_forcpu, sendcnts, displs, MPI_INT, row_indices_fromcpu[i], numvalues_fromcpu[i], MPI_INT, i, comm);
    … …  
    351369                                MPI_Scatterv( modes_forcpu, sendcnts, displs, MPI_INT, modes_fromcpu[i], numvalues_fromcpu[i], MPI_INT, i, comm);
    352370                        }
    353                        
    354                         /*Plug values into global matrix: */
     371                        /*}}}*/
     372
     373                        /*Plug values into global matrix: {{{*/
    355374                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,comm);
    356375                        for(i=0;i<num_procs;i++){
    … …  
    366385                                }
    367386                        }
    368                        
     387                        /*}}}*/
     388
    369389                        /*Free ressources:{{{*/
    370390                        xDelete<int>(RowRank);
    … …  
    523543
    524544                        /*figure out size of buffers per cpu: */
     545
     546                        numvalues_forcpu=xNew<int>(num_procs);
    525547                        for(i=0;i<num_procs;i++){
    526548                                DataSet    *buckets            = bucketsforcpu[i];
    … …  
    565587                        }
    566588
     589                        /*sanity check: */
     590                        if (temp_row_indices_forcpu!=row_indices_forcpu+total_size)_error_("problem with marshalling of buckets");
     591                        if (temp_col_indices_forcpu!=col_indices_forcpu+total_size)_error_("problem with marshalling of buckets");
     592                        if (temp_values_forcpu!=values_forcpu+total_size)_error_("problem with marshalling of buckets");
     593                        if (temp_modes_forcpu!=modes_forcpu+total_size)_error_("problem with marshalling of buckets");
     594
    567595                        /*output buffers: */
    568                         *pnumvalues_forcpu   = row_indices_forcpu;
     596                        *pnumvalues_forcpu   = numvalues_forcpu;
    569597                        *prow_indices_forcpu = row_indices_forcpu;
    570598                        *pcol_indices_forcpu = col_indices_forcpu;
    … …  
    573601                }
    574602                /*}}}*/         
     603                /*Solve{{{*/
     604                IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pfin, Parameters* parameters){
     605
     606                        printf("IssmMpiDenseMat solve not implemented yet!");
     607                        exit(1);
     608
     609                }/*}}}*/
    575610};
    576611

  • issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h

    r14822 r14834  
    126126                /*FUNCTION Assemble{{{*/
    127127                void Assemble(){
     128
     129
     130                        int           i,j;
     131
     132                        int         *RowRank            = NULL;
     133                        int           num_procs;
     134
     135                        int        *row_indices_forcpu = NULL;
     136                        int        *col_indices_forcpu = NULL;
     137                        int        *modes_forcpu       = NULL;
     138                        doubletype *values_forcpu      = NULL;
     139                        int         *numvalues_forcpu   = NULL;
     140                        DataSet     **bucketsforcpu       = NULL;
     141
     142                        int        **row_indices_fromcpu = NULL;
     143                        int        **col_indices_fromcpu = NULL;
     144                        int        **modes_fromcpu       = NULL;
     145                        doubletype **values_fromcpu      = NULL;
     146                        int         *numvalues_fromcpu   = NULL;
     147
     148                        int           lower_row;
     149                        int           upper_row;
     150                        int*          sendcnts            = NULL;
     151                        int*          displs              = NULL;
     152                        int           count               = 0;
     153
     154                        /*some communicator info: */
     155                        num_procs=IssmComm::GetSize();
     156                        MPI_Comm comm=IssmComm::GetComm();
     157
     158
     159                        /*First, make a vector of size M, which for each row between 0 and M-1, tells which cpu this row belongs to: */
     160                        RowRank=DetermineRowRankFromLocalSize(M,m,comm);
     161
     162
     163                        /*Now, sort out our dataset of buckets according to cpu ownership of rows: {{{*/
     164                        bucketsforcpu=xNew<DataSet*>(num_procs);
     165
     166                        for(i=0;i<num_procs;i++){
     167                                DataSet* bucketsofcpu_i=new DataSet();
     168                                for (j=0;j<buckets->Size();j++){
     169                                        Bucket<doubletype>* bucket=(Bucket<doubletype>*)buckets->GetObjectByOffset(j);
     170                                        bucket->SpawnBucketsPerCpu(bucketsofcpu_i,i,RowRank);
     171                                }
     172                                bucketsforcpu[i]=bucketsofcpu_i;
     173                        }
     174                        /*}}}*/
     175
     176                        /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this  {{{
     177                         * dataset owns correspond to rows that are owned by cpu i, not j!. Out of all the buckets we own, make row,col,value,insert_mode
     178                         * vectors that will be shipped around the cluster: */
     179                        this->BucketsBuildScatterBuffers(&numvalues_forcpu,&row_indices_forcpu,&values_forcpu,&modes_forcpu,bucketsforcpu,num_procs);
     180                        /*}}}*/
     181
     182                        /*Now, we need to allocate on each cpu arrays to receive data from all the other cpus. To know what we need to allocate, we need  {{{
     183                         *some scatter calls: */
     184                        numvalues_fromcpu   = xNew<int>(num_procs);
     185                        for(i=0;i<num_procs;i++){
     186                                MPI_Scatter(numvalues_forcpu,1,MPI_INT,numvalues_fromcpu+i,1,MPI_INT,i,comm);
     187                        }
     188                       
     189                        row_indices_fromcpu=xNew<int*>(num_procs);
     190                        values_fromcpu=xNew<doubletype*>(num_procs);
     191                        modes_fromcpu=xNew<int*>(num_procs);
     192                        for(i=0;i<num_procs;i++){
     193                                int size=numvalues_fromcpu[i];
     194                                if(size){
     195                                        row_indices_fromcpu[i]=xNew<int>(size);
     196                                        values_fromcpu[i]=xNew<doubletype>(size);
     197                                        modes_fromcpu[i]=xNew<int>(size);
     198                                }
     199                                else{
     200                                        row_indices_fromcpu[i]=NULL;
     201                                        values_fromcpu[i]=NULL;
     202                                        modes_fromcpu[i]=NULL;
     203                                }
     204                        }
     205                        /*}}}*/
     206
     207                        /*Scatter values around: {{{*/
     208                        /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given
     209                         * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the MPI_Scatterv prototype:
     210                         * int MPI_Scatterv( void *sendbuf, int *sendcnts, int *displs, MPI_Datatype sendtype, void *recvbuf, int recvcnt, MPI_Datatype recvtype, int root, MPI_Comm comm) :*/
     211                        sendcnts=xNew<int>(num_procs);
     212                        displs=xNew<int>(num_procs);
     213                        count=0;
     214                        for(i=0;i<num_procs;i++){
     215                                sendcnts[i]=numvalues_forcpu[i];
     216                                displs[i]=count;
     217                                count+=numvalues_forcpu[i];
     218                        }
     219
     220                        for(i=0;i<num_procs;i++){
     221                                MPI_Scatterv( row_indices_forcpu, sendcnts, displs, MPI_INT, row_indices_fromcpu[i], numvalues_fromcpu[i], MPI_INT, i, comm);
     222                                MPI_Scatterv( values_forcpu, sendcnts, displs, MPI_DOUBLE, values_fromcpu[i], numvalues_fromcpu[i], MPI_DOUBLE, i, comm);
     223                                MPI_Scatterv( modes_forcpu, sendcnts, displs, MPI_INT, modes_fromcpu[i], numvalues_fromcpu[i], MPI_INT, i, comm);
     224                        }
     225                        /*}}}*/
     226
     227                        /*Plug values into global vector: {{{*/
     228                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,comm);
     229                        for(i=0;i<num_procs;i++){
     230                                int  numvalues=numvalues_fromcpu[i];
     231                                int* rows=row_indices_fromcpu[i];
     232                                doubletype* values=values_fromcpu[i];
     233                                int* mods=modes_fromcpu[i];
     234
     235                                for(j=0;j<numvalues;j++){
     236                                        if(mods[j]==ADD_VAL) *(vector+(rows[j]-lower_row))+=values[j];
     237                                        else *(vector+(rows[j]-lower_row))=values[j];
     238                                }
     239                        }
     240                        /*}}}*/
     241
     242                        /*Free ressources:{{{*/
     243                        xDelete<int>(RowRank);
     244                        xDelete<int>(row_indices_forcpu);
     245                        xDelete<int>(modes_forcpu);
     246                        xDelete<doubletype>(values_forcpu);
     247                        xDelete<int>(numvalues_forcpu);
     248                       
     249                        for(i=0;i<num_procs;i++){
     250                                DataSet* buckets=bucketsforcpu[i];
     251                                delete buckets;
     252                        }
     253                        xDelete<DataSet*>(bucketsforcpu);
     254
     255                        for(i=0;i<num_procs;i++){
     256                                int* rows=row_indices_fromcpu[i];
     257                                int* modes=modes_fromcpu[i];
     258                                doubletype* values=values_fromcpu[i];
     259
     260                                xDelete<int>(rows);
     261                                xDelete<int>(modes);
     262                                xDelete<doubletype>(values);
     263                        }
     264                        xDelete<int*>(row_indices_fromcpu);
     265                        xDelete<int*>(modes_fromcpu);
     266                        xDelete<doubletype*>(values_fromcpu);
     267                        xDelete<int>(numvalues_fromcpu);
     268                       
     269                        xDelete<int>(sendcnts);
     270                        xDelete<int>(displs);
     271                        /*}}}*/
     272
     273
     274                }
     275                /*}}}*/
     276                /*FUNCTION Assemble2{{{*/
     277                void Assemble2(){
    128278
    129279                        int           i;
    … …  
    461611                }
    462612                /*}}}*/
     613                /*FUNCTION BucketsBuildScatterBuffers{{{*/
     614                void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){
     615
     616
     617                        /*intermediary: */
     618                        int         i,j;
     619                        int         count                   = 0;
     620                        int         total_size              = 0;
     621                        int        *temp_row_indices_forcpu = NULL;
     622                        doubletype *temp_values_forcpu      = NULL;
     623                        int        *temp_modes_forcpu       = NULL;
     624
     625                        /*output: */
     626                        int        *numvalues_forcpu        = NULL;
     627                        int        *row_indices_forcpu      = NULL;
     628                        doubletype *values_forcpu           = NULL;
     629                        int        *modes_forcpu            = NULL;
     630
     631                        /*figure out size of buffers per cpu: */
     632
     633                        numvalues_forcpu=xNew<int>(num_procs);
     634                        for(i=0;i<num_procs;i++){
     635                                DataSet    *buckets            = bucketsforcpu[i];
     636                               
     637                                count=0;
     638                                for(j=0;j<buckets->Size();j++){
     639                                        Bucket<doubletype>* bucket =(Bucket<doubletype>*)buckets->GetObjectByOffset(j);
     640                                        count+=bucket->MarshallSize();
     641                                }
     642
     643                                numvalues_forcpu[i]=count;
     644                        }
     645
     646                        /*now, figure out size of  total buffers (for all cpus!): */
     647                        count=0;
     648                        for(i=0;i<num_procs;i++){
     649                                count+=numvalues_forcpu[i];
     650                        }
     651                        total_size=count;
     652
     653                        /*Allocate buffers: */
     654                        row_indices_forcpu = xNew<int>(total_size);
     655                        values_forcpu = xNew<doubletype>(total_size);
     656                        modes_forcpu = xNew<int>(total_size);
     657
     658                        /*we are going to march through the buffers, and marshall data onto them, so in order to not
     659                         *lose track of where these buffers are located in memory, we are going to work using copies
     660                         of them: */
     661                        temp_row_indices_forcpu=row_indices_forcpu;
     662                        temp_values_forcpu=values_forcpu;
     663                        temp_modes_forcpu=modes_forcpu;
     664
     665                        /*Fill buffers: */
     666                        for(i=0;i<num_procs;i++){
     667                                DataSet    *buckets            = bucketsforcpu[i];
     668                                for(j=0;j<buckets->Size();j++){
     669                                        Bucket<doubletype>* bucket =(Bucket<doubletype>*)buckets->GetObjectByOffset(j);
     670                                        bucket->Marshall(&temp_row_indices_forcpu,&temp_values_forcpu,&temp_modes_forcpu); //pass in the address of the buffers, so as to have the Marshall routine increment them.
     671                                }
     672                        }
     673
     674                        /*sanity check: */
     675                        if (temp_row_indices_forcpu!=row_indices_forcpu+total_size)_error_("problem with marshalling of buckets");
     676                        if (temp_values_forcpu!=values_forcpu+total_size)_error_("problem with marshalling of buckets");
     677                        if (temp_modes_forcpu!=modes_forcpu+total_size)_error_("problem with marshalling of buckets");
     678
     679                        /*output buffers: */
     680                        *pnumvalues_forcpu   = numvalues_forcpu;
     681                        *prow_indices_forcpu = row_indices_forcpu;
     682                        *pvalues_forcpu      = values_forcpu;
     683                        *pmodes_forcpu       = modes_forcpu;
     684                }
     685                /*}}}*/         
    463686};
    464687#endif //#ifndef _ISSM_MPI_VEC_H_       

  • issm/trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp

    r14792 r14834  
    1 /*!\file SolverxSeq
    2  * \brief implementation of sequential solver using the GSL librarie
     1/*!\file IssmSolver
     2 * \brief implementation of solver
    33 */
    44
    … …  
    2020#endif
    2121
    22 void IssmSolve(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kffin, IssmVec<IssmDouble>* pfin, Parameters* parameters){/*{{{*/
    23 
    24         /*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */
    25         IssmSeqVec<IssmDouble>* pf=(IssmSeqVec<IssmDouble>*)pfin->vector;
    26         IssmDenseMat<IssmDouble>* Kff=(IssmDenseMat<IssmDouble>*)Kffin->matrix;
    27 
    28 #ifdef _HAVE_GSL_
    29         /*Intermediary: */
    30         int M,N,N2;
    31         IssmSeqVec<IssmDouble> *uf = NULL;
    32 
    33         Kff->GetSize(&M,&N);
    34         pf->GetSize(&N2);
    35 
    36         if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !");
    37         if(M!=N)_error_("Stiffness matrix should be square!");
    38 #ifdef _HAVE_ADOLC_
    39         ensureContiguousLocations(N);
    40 #endif
    41         IssmDouble *x  = xNew<IssmDouble>(N);
    42 #ifdef _HAVE_ADOLC_
    43         SolverxSeq(x,Kff->matrix,pf->vector,N,parameters);
    44 #else
    45         SolverxSeq(x,Kff->matrix,pf->vector,N);
    46 #endif
    47 
    48         uf=new IssmSeqVec<IssmDouble>(x,N);
    49         xDelete(x);
    50 
    51         /*Assign output pointers:*/
    52         IssmVec<IssmDouble>* out=new IssmVec<IssmDouble>();
    53         out->vector=(IssmAbsVec<IssmDouble>*)uf;
    54         *pout=out;
    55 
    56 #else
    57         _error_("GSL support not compiled in!");
    58 #endif
    59 
    60 }/*}}}*/
     22void IssmSolve(IssmVec<IssmDouble>** pout,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf, Parameters* parameters){/*{{{*/
     23
     24        /*Let matrix decide, to retain object orientation: */
     25        IssmVec<IssmDouble>* outvector=NULL;
     26
     27        outvector=Kff->Solve(pf,parameters);
     28
     29        *pout=outvector;
     30}
     31/*}}}*/
    6132void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
    6233
Note: See TracChangeset for help on using the changeset viewer.