Changeset 11734


Ignore:
Timestamp:
03/19/12 20:11:39 (13 years ago)
Author:
Eric.Larour
Message:

Hook up new issm toolkit

Location:
issm/trunk-jpl/src/c/objects/Numerics
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp

    r11679 r11734  
    260260        this->CheckConsistency();
    261261
    262         #ifdef _HAVE_PETSC_
    263                 if(this->dofsymmetrical){
    264                         /*only use row dofs to add values into global matrices: */
    265                        
    266                         if(this->row_fsize){
    267                                 /*first, retrieve values that are in the f-set from the g-set values matrix: */
    268                                 localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
    269                                 for(i=0;i<this->row_fsize;i++){
    270                                         for(j=0;j<this->row_fsize;j++){
    271                                                 *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
    272                                         }
     262        if(this->dofsymmetrical){
     263                /*only use row dofs to add values into global matrices: */
     264
     265                if(this->row_fsize){
     266                        /*first, retrieve values that are in the f-set from the g-set values matrix: */
     267                        localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
     268                        for(i=0;i<this->row_fsize;i++){
     269                                for(j=0;j<this->row_fsize;j++){
     270                                        *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
    273271                                }
    274                                 /*add local values into global  matrix, using the fglobaldoflist: */
    275                                 MatSetValues(Kff->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES);
    276 
    277                                 /*Free ressources:*/
    278                                 xfree((void**)&localvalues);
    279                         }
    280 
    281 
    282                         if((this->row_ssize!=0) && (this->row_fsize!=0)){
    283                                 /*first, retrieve values that are in the f and s-set from the g-set values matrix: */
    284                                 localvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double));
    285                                 for(i=0;i<this->row_fsize;i++){
    286                                         for(j=0;j<this->row_ssize;j++){
    287                                                 *(localvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]);
    288                                         }
     272                        }
     273                        /*add local values into global  matrix, using the fglobaldoflist: */
     274                        Kff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
     275
     276                        /*Free ressources:*/
     277                        xfree((void**)&localvalues);
     278                }
     279
     280
     281                if((this->row_ssize!=0) && (this->row_fsize!=0)){
     282                        /*first, retrieve values that are in the f and s-set from the g-set values matrix: */
     283                        localvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double));
     284                        for(i=0;i<this->row_fsize;i++){
     285                                for(j=0;j<this->row_ssize;j++){
     286                                        *(localvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]);
    289287                                }
    290                                 /*add local values into global  matrix, using the fglobaldoflist: */
    291                                 MatSetValues(Kfs->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,(const double*)localvalues,ADD_VALUES);
    292 
    293                                 /*Free ressources:*/
    294                                 xfree((void**)&localvalues);
    295                         }
    296                 }
    297                 else{
    298                         _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
    299                 }
    300         #else
    301                 _error_("not supported yet!");
    302         #endif
     288                        }
     289                        /*add local values into global  matrix, using the fglobaldoflist: */
     290                        Kfs->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,localvalues,ADD_VAL);
     291
     292                        /*Free ressources:*/
     293                        xfree((void**)&localvalues);
     294                }
     295        }
     296        else{
     297                _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
     298        }
    303299
    304300}
    … …  
    316312        this->CheckConsistency();
    317313
    318         #ifdef _HAVE_PETSC_
    319                 if(this->dofsymmetrical){
    320                         /*only use row dofs to add values into global matrices: */
    321 
    322                         if(this->row_fsize){
    323                                 /*first, retrieve values that are in the f-set from the g-set values matrix: */
    324                                 localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
    325                                 for(i=0;i<this->row_fsize;i++){
    326                                         for(j=0;j<this->row_fsize;j++){
    327                                                 *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
    328                                         }
     314        if(this->dofsymmetrical){
     315                /*only use row dofs to add values into global matrices: */
     316
     317                if(this->row_fsize){
     318                        /*first, retrieve values that are in the f-set from the g-set values matrix: */
     319                        localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double));
     320                        for(i=0;i<this->row_fsize;i++){
     321                                for(j=0;j<this->row_fsize;j++){
     322                                        *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]);
    329323                                }
    330                                 /*add local values into global  matrix, using the fglobaldoflist: */
    331                                 MatSetValues(Jff->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES);
    332 
    333                                 /*Free ressources:*/
    334                                 xfree((void**)&localvalues);
    335                         }
    336 
    337                 }
    338                 else{
    339                         _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
    340                 }
    341         #else
    342                 _error_("not supported yet!");
    343         #endif
     324                        }
     325                        /*add local values into global  matrix, using the fglobaldoflist: */
     326                        Jff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
     327
     328                        /*Free ressources:*/
     329                        xfree((void**)&localvalues);
     330                }
     331
     332        }
     333        else{
     334                _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!");
     335        }
    344336
    345337}

  • issm/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp

    r11679 r11734  
    166166        double* localvalues=NULL;
    167167
    168         #ifdef _HAVE_PETSC_
    169                 if(this->fsize){
    170                         /*first, retrieve values that are in the f-set from the g-set values vector: */
    171                         localvalues=(double*)xmalloc(this->fsize*sizeof(double));
    172                         for(i=0;i<this->fsize;i++){
    173                                 localvalues[i]=this->values[this->flocaldoflist[i]];
    174                         }
    175                         /*add local values into global  vector, using the fglobaldoflist: */
    176                         VecSetValues(pf->vector,this->fsize,this->fglobaldoflist,(const double*)localvalues,ADD_VALUES);
    177 
    178                         /*Free ressources:*/
    179                         xfree((void**)&localvalues);
    180                 }
    181         #else
    182                 _error_("not supported yet!");
    183         #endif
     168        if(this->fsize){
     169                /*first, retrieve values that are in the f-set from the g-set values vector: */
     170                localvalues=(double*)xmalloc(this->fsize*sizeof(double));
     171                for(i=0;i<this->fsize;i++){
     172                        localvalues[i]=this->values[this->flocaldoflist[i]];
     173                }
     174                /*add local values into global  vector, using the fglobaldoflist: */
     175                pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL);
     176
     177                /*Free ressources:*/
     178                xfree((void**)&localvalues);
     179        }
    184180       
    185181}
    … …  
    191187        double* localvalues=NULL;
    192188
    193         #ifdef _HAVE_PETSC_
    194                 if(this->fsize){
    195                         /*first, retrieve values that are in the f-set from the g-set values vector: */
    196                         localvalues=(double*)xmalloc(this->fsize*sizeof(double));
    197                         for(i=0;i<this->fsize;i++){
    198                                 localvalues[i]=this->values[this->flocaldoflist[i]];
    199                         }
    200                         /*add local values into global  vector, using the fglobaldoflist: */
    201                         VecSetValues(pf->vector,this->fsize,this->fglobaldoflist,(const double*)localvalues,INSERT_VALUES);
    202 
    203                         /*Free ressources:*/
    204                         xfree((void**)&localvalues);
    205                 }
    206         #else
    207                 _error_("not supported yet!");
    208         #endif
     189        if(this->fsize){
     190                /*first, retrieve values that are in the f-set from the g-set values vector: */
     191                localvalues=(double*)xmalloc(this->fsize*sizeof(double));
     192                for(i=0;i<this->fsize;i++){
     193                        localvalues[i]=this->values[this->flocaldoflist[i]];
     194                }
     195                /*add local values into global  vector, using the fglobaldoflist: */
     196                pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,INS_VAL);
     197
     198                /*Free ressources:*/
     199                xfree((void**)&localvalues);
     200        }
    209201
    210202}

  • issm/trunk-jpl/src/c/objects/Numerics/Matrix.cpp

    r11713 r11734  
    4848        this->matrix=NewMat(pM,pN);
    4949        #else
    50         this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
     50        this->matrix=new SeqMat(pM,pN);
    5151        #endif
    5252        #ifdef _HAVE_ADOLC_
    … …  
    6464        this->matrix=NewMat(pM,pN,sparsity);
    6565        #else
    66         this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
     66        this->matrix=new SeqMat(pM,pN,sparsity);
    6767        #endif
    6868        #ifdef _HAVE_ADOLC_
    … …  
    9797        xfree((void**)&idxn);
    9898        #else
    99         this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
     99        this->matrix=new SeqMat(serial_mat,pM,pN,sparsity);
    100100        #endif
    101101        #ifdef _HAVE_ADOLC_
    … …  
    113113        this->matrix=NewMat(pM,pN,connectivity,numberofdofspernode);
    114114        #else
    115         this->matrix=(double*)xcalloc(pM*pN,sizeof(double));
     115        this->matrix=new SeqMat(pM,pN,connectivity,numberofdofspernode);
    116116        #endif
    117117        #ifdef _HAVE_ADOLC_
    … …  
    126126        MatFree(&this->matrix);
    127127        #else
    128         xfree((void**)&this->matrix);
     128        delete this->matrix;
    129129        #endif
    130130        #ifdef _HAVE_ADOLC_
    … …  
    143143        MatView(this->matrix,PETSC_VIEWER_STDOUT_WORLD);
    144144        #else
    145         printf("Matrix size: %i-%i\n",M,N);
    146         for(i=0;i<M;i++){
    147                 for(j=0;j<N;j++){
    148                         printf("%g ",*(matrix+N*i+j));
    149                 }
    150                 printf("\n");
    151         }
     145        this->matrix->Echo();
    152146        #endif
    153147
    … …  
    173167        PetscMatrixToMatlabMatrix(&dataref,this->matrix);
    174168        #else
    175         _error_("not implemented yet!");
     169        dataref=this->matrix->ToMatlabMatrix();
    176170        #endif
    177171        return dataref;
    178172
     173}
     174/*}}}*/
     175/*FUNCTION MatlabMatrixToMatrix{{{1*/
     176Matrix* MatlabMatrixToMatrix(const mxArray* mxmatrix){
     177
     178        int dummy;
     179        Matrix* matrix=NULL;
     180
     181        /*allocate matrix object: */
     182        matrix=new Matrix();
     183
     184        #ifdef _HAVE_PETSC_
     185        MatlabMatrixToPetscMatrix(&matrix->matrix,&matrix->M,&matrix->N,mxmatrix);
     186        #else
     187        matrix->matrix=MatlabMatrixToSeqMat(mxmatrix);
     188        #endif
     189       
     190        return matrix;
    179191}
    180192/*}}}*/
    … …  
    190202                #endif
    191203        #else
    192                 /*do nothing:*/
     204                this->matrix->Assemble();
    193205        #endif
    194206
    … …  
    203215                MatNorm(this->matrix,ISSMToPetscNormMode(norm_type),&norm);
    204216        #else
    205                 _error_("not implemented yet!");
     217                norm=this->matrix->Norm(norm_type);
    206218        #endif
    207219        return norm;
    … …  
    215227                MatGetSize(this->matrix,pM,pN);
    216228        #else
    217                 _error_("not implemented yet!");
     229                this->matrix->GetSize(pM,pN);
    218230        #endif
    219231}
    … …  
    226238                MatGetLocalSize(this->matrix,pM,pN);
    227239        #else
    228                 _error_("not implemented yet!");
     240                this->matrix->GetLocalSize(pM,pN);
    229241        #endif
    230242}
    … …  
    238250                MatMultPatch(this->matrix,X->vector,AX->vector);
    239251        #else
    240                 _error_("not implemented yet!");
     252                this->matrix->MatMult(X->vector,AX->vector);
    241253        #endif
    242254}
    … …  
    249261        output=new Matrix();
    250262
    251         #ifdef _HAVE_PETSC
     263        #ifdef _HAVE_PETSC_
    252264                _assert_(this->matrix);
    253265                MatDuplicate(this->matrix,MAT_COPY_VALUES,&output->matrix);
    254266        #else
    255                 _error_("not implemented yet!");
     267                output->matrix=this->matrix->Duplicate();
    256268        #endif
    257269}
    … …  
    262274        double* output=NULL;
    263275
    264         #ifdef _HAVE_PETSC
     276        #ifdef _HAVE_PETSC_
    265277                MatToSerial(&output,this->matrix);
    266278        #else
    267                 _error_("not implemented yet!");
     279                output=this->matrix->ToSerial();
    268280        #endif
    269281        return output;
    270282}
    271283/*}}}*/
     284/*FUNCTION Matrix::SetValues{{{1*/
     285void Matrix::SetValues(int m,int* idxm,int n,int* idxn,double* values,InsMode mode){
     286
     287        #ifdef _HAVE_PETSC_
     288                MatSetValues(this->matrix,m,idxm,n,idxn,values,ISSMToPetscInsertMode(mode));
     289        #else
     290                this->matrix->SetValues(m,idxm,n,idxn,values,mode);
     291        #endif
     292}
     293/*}}}*/

  • issm/trunk-jpl/src/c/objects/Numerics/Matrix.h

    r11695 r11734  
    3737                Mat matrix;
    3838                #else
    39                 double* matrix;
     39                SeqMat* matrix;
    4040                #endif
    4141                #ifdef _HAVE_ADOLC_
    … …  
    6363                Matrix* Duplicate(void);
    6464                double* ToSerial(void);
     65                void SetValues(int m,int* idxm,int n,int* idxn,double* values,InsMode mode);
    6566                /*}}}*/
    6667
    6768};
     69/*API: */
     70#ifdef _SERIAL_
     71Matrix* MatlabMatrixToMatrix(const mxArray* mxmatrix);
     72#endif
     73
    6874#endif //#ifndef _MATRIX_H_

  • issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp

    r11713 r11734  
    2525Vector::Vector(){
    2626
    27         this->M=0;
    2827        #ifdef _HAVE_PETSC_
    2928        this->vector=NULL;
    … …  
    4241        this->vector=NewVec(pM,fromlocalsize);
    4342        #else
    44         this->M=pM;
    45         this->vector=(double*)xcalloc(pM,sizeof(double));
     43        this->vector=new SeqVec(pM,fromlocalsize);
    4644        #endif
    4745        #ifdef _HAVE_ADOLC_
    … …  
    5149/*}}}*/
    5250/*FUNCTION Vector::Vector(double* serial_vec,int M){{{1*/
    53 Vector::Vector(double* serial_vec,int pM){
     51Vector::Vector(double* serial_vec,int M){
    5452
    5553        int i,j;
    … …  
    5856                int* idxm=NULL;
    5957
    60                 this->vector=NewVec(pM);
     58                this->vector=NewVec(M);
    6159               
    6260                idxm=(int*)xmalloc(M*sizeof(int));
    … …  
    7068
    7169        #else
    72                 this->M=pM;
    73                 this->vector=(double*)xcalloc(pM,sizeof(double));
    74         #endif
    75         #ifdef _HAVE_ADOLC_
    76                 this->avector=(adouble*)xmalloc(pM*sizeof(adouble));
     70                this->vector=new SeqVec(serial_vec,M);
     71        #endif
     72        #ifdef _HAVE_ADOLC_
     73                this->avector=(adouble*)xmalloc(M*sizeof(adouble));
    7774        #endif
    7875}
    … …  
    8279Vector::Vector(Vec petsc_vec){
    8380
    84         /*Get Vector size*/
    85         VecGetSize(petsc_vec,&this->M);
    86 
    8781        /*copy vector*/
    8882        VecDuplicate(petsc_vec,&this->vector);
    … …  
    9892        VecFree(&this->vector);
    9993        #else
    100         xfree((void**)&this->vector);
     94        delete this->vector;
    10195        #endif
    10296        #ifdef _HAVE_ADOLC_
    … …  
    115109        VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD);
    116110        #else
    117         printf("Vector size: %i\n",M);
    118         for(i=0;i<M;i++){
    119                 printf("%g\n ",*(vector+i));
    120         }
    121         #endif
    122 
    123         #ifdef _HAVE_ADOLC_
    124         /*Not sure about that one. Should we use the overloaded operator >>?*/
    125         printf("ADOLC Vector equivalent:" );
    126         for(i=0;i<M;i++){
    127                 printf("%g\n ",*(avector+i));
    128         }
     111        this->vector->Echo();
     112        #endif
     113
     114        #ifdef _HAVE_ADOLC_
     115        /*do nothing for now: */
    129116        #endif
    130117}
    … …  
    139126        PetscVectorToMatlabVector(&dataref,this->vector);
    140127        #else
    141         _error_("not implemented yet!");
     128        dataref=this->vector->ToMatlabVector();
    142129        #endif
    143130        return dataref;
    144131
     132}
     133/*}}}*/
     134/*FUNCTION MatlabVectorToVector{{{1*/
     135Vector* MatlabVectorToVector(const mxArray* mxvector){
     136
     137        int dummy;
     138        Vector* vector=NULL;
     139
     140        /*allocate vector object: */
     141        vector=new Vector();
     142
     143        #ifdef _HAVE_PETSC_
     144        MatlabVectorToPetscVector(&vector->vector,&dummy,mxvector);
     145        #else
     146        vector->vector=MatlabVectorToSeqVec(mxvector);
     147        #endif
     148       
     149        return vector;
    145150}
    146151/*}}}*/
    … …  
    154159                VecAssemblyEnd(this->vector);
    155160        #else
    156                 /*do nothing*/
     161                this->vector->Assemble();
    157162        #endif
    158163
    … …  
    167172                VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode));
    168173        #else
    169                 _error_("not implemented yet!");
     174                this->vector->SetValues(ssize,list,values,mode);
    170175        #endif
    171176
    … …  
    179184                VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode));
    180185        #else
    181                 _error_("not implemented yet!");
     186                this->vector->SetValue(dof,value,mode);
    182187        #endif
    183188
    … …  
    191196                VecGetValues(this->vector,1,&dof,pvalue);
    192197        #else
    193                 _error_("not implemented yet!");
     198        this->vector->GetValue(pvalue,dof);
    194199        #endif
    195200}
    … …  
    202207                VecGetSize(this->vector,pM);
    203208        #else
    204                 _error_("not implemented yet!");
     209                this->vector->GetSize(pM);
    205210        #endif
    206211
    … …  
    214219                VecGetLocalSize(this->vector,pM);
    215220        #else
    216                 _error_("not implemented yet!");
     221                this->vector->GetLocalSize(pM);
    217222        #endif
    218223
    … …  
    229234                VecDuplicate(this->vector,&vec_output);
    230235                output->vector=vec_output;
    231                 VecGetSize(output->vector,&output->M);
    232         #else
    233                 _error_("not implemented yet!");
     236        #else
     237                output->vector=this->vector->Duplicate();
    234238        #endif
    235239        return output;
    … …  
    244248                VecSet(this->vector,value);
    245249        #else
    246                 _error_("not implemented yet!");
     250                this->vector->Set(value);
    247251        #endif
    248252
    … …  
    256260                VecAXPY(this->vector,a,X->vector);
    257261        #else
    258                 _error_("not implemented yet!");
     262                this->vector->AXPY(X->vector,a);
    259263        #endif
    260264}
    … …  
    267271                VecAYPX(this->vector,a,X->vector);
    268272        #else
    269                 _error_("not implemented yet!");
     273                this->vector->AYPX(X->vector,a);
    270274        #endif
    271275
    … …  
    280284                VecToMPISerial(&vec_serial, this->vector);
    281285        #else
    282                 _error_("not implemented yet!");
     286                vec_serial=this->vector->ToMPISerial();
    283287        #endif
    284288
    … …  
    293297                VecCopy(this->vector,to->vector);
    294298        #else
    295                 _error_("not implemented yet!");
     299                this->vector->Copy(to->vector);
    296300        #endif
    297301
    … …  
    306310                VecNorm(this->vector,ISSMToPetscNormMode(norm_type),&norm);
    307311        #else
    308                 _error_("not implemented yet!");
     312                norm=this->vector->Norm(norm_type);
    309313        #endif
    310314        return norm;
    … …  
    318322                VecScale(this->vector,scale_factor);
    319323        #else
    320                 _error_("not implemented yet!");
     324                this->vector->Scale(scale_factor);
    321325        #endif
    322326}
    … …  
    330334                VecDot(this->vector,vector->vector,&dot);
    331335        #else
    332                 _error_("not implemented yet!");
     336                dot=this->vector->Dot(vector->vector);
    333337        #endif
    334338        return dot;
    … …  
    342346                VecPointwiseDivide(this->vector,x->vector,y->vector);
    343347        #else
    344                 _error_("not implemented yet!");
    345         #endif
    346 }
    347 /*}}}*/
     348                this->vector->PointwiseDivide(x->vector,y->vector);
     349        #endif
     350}
     351/*}}}*/

  • issm/trunk-jpl/src/c/objects/Numerics/Vector.h

    r11704 r11734  
    3131        public:
    3232       
    33                 int M;
    34 
    3533                #ifdef _HAVE_PETSC_
    3634                Vec vector;
    3735                #else
    38                 double* vector;
     36                SeqVec* vector;
    3937                #endif
    4038                #ifdef _HAVE_ADOLC_
    … …  
    7472                /*}}}*/
    7573};
     74
     75/*API: */
     76#ifdef _SERIAL_
     77Vector* MatlabVectorToVector(const mxArray* mxvector);
     78#endif
     79
    7680#endif //#ifndef _VECTOR_H_
Note: See TracChangeset for help on using the changeset viewer.