Changeset 14875 for issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h

Timestamp:

05/03/13 15:56:49 (12 years ago)

Author:

Eric.Larour

Message:

CHG: new NORM_FROBENIUS for matrices.
Also fixed a lot of bugs in the norm routines of vectors and matrices in our issm toolkit.
Took out the Assemble2, along with Bucket.h iSend and Recv routiens.

File:

• issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (modified) (3 diffs)

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

@@ -295 +295 @@
                 }
                 /*}}}*/
-                /*FUNCTION Assemble2{{{*/
-                void Assemble2(){
-
-                        int           i;
-                        int           j;
-                        int           k;
-                        int           my_rank;
-                        int           num_procs;
-                        int          *RowRank             = NULL;
-
-                        DataSet     **bucketspercpu       = NULL;
-                        int          *bucketspercpu_sizes = NULL;
-                        MPI_Request  *requests            = NULL;
-                        MPI_Status   *statuses            = NULL;
-                        MPI_Status    status;
-                        int           num_requests        = 0;
-                        DataSet      *mybuckets           = NULL;
-                        int           lower_row;
-                        int           upper_row;
-                        int           count               = 0;
-
-                        int           size;
-
-                        /*some communicator info: */
-                        num_procs=IssmComm::GetSize();
-                        my_rank=IssmComm::GetRank();
-                        MPI_Comm comm=IssmComm::GetComm();
-
-                        /*First, make a vector of size M, which for each row between 0 and M-1, tells which cpu this row belongs to: */
-                        RowRank=DetermineRowRankFromLocalSize(M,m,comm);
-
-                        /*Now, sort out our dataset of buckets according to cpu ownership of rows: */
-                        bucketspercpu=xNew<DataSet*>(num_procs);
-                        bucketspercpu_sizes=xNew<int>(num_procs);
-
-                        for(i=0;i<num_procs;i++){
-                                DataSet* bucketsofcpu_i=new DataSet();
-                                for (j=0;j<buckets->Size();j++){
-                                        Bucket<doubletype>* bucket=(Bucket<doubletype>*)buckets->GetObjectByOffset(j);
-                                        bucket->SpawnBucketsPerCpu(bucketsofcpu_i,i,RowRank);
-                                }
-                                bucketspercpu[i]=bucketsofcpu_i;
-                                bucketspercpu_sizes[i]=bucketsofcpu_i->Size();
-                        }
-
-                        /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this 
-                         * dataset owns correspond to rows that are owned by cpu i, not j!:*/
-
-                        /*First, figure out how many requests are going to be sent by MPI_Isend. Do this a little bit better? */
-                        for(i=0;i<num_procs;i++){
-                                if(i!=my_rank){
-                                        num_requests+=bucketspercpu[i]->Size()*VECTORBUCKETSIZEOFREQUESTS; //this is to take into account all the MPI_ISend calls in each bucket.
-                                        num_requests++; //this is to take into account on MPI_ISend in BucketsSend.
-                                }
-                        }
-
-                        /*Initialize array to track requests and statuses: */
-                        requests=new MPI_Request[num_requests];
-                        statuses=new MPI_Status[num_requests];
-
-                        /*Now, go through all our bucketspercpu datasets, and send them to the corresponding cpus. Do not send our own buckets though!: */
-                        count=0; //count requests
-                        for(i=0;i<num_procs;i++){
-                                if(my_rank==i){
-                                        for(j=0;j<num_procs;j++){
-                                                if(j!=i){//only send the buckets that this cpu does not own.
-
-                                                        /*Go through the buckets belonging to cpu j, and send them accordingly. */
-                                                        DataSet* buckets=bucketspercpu[j];
-                                                        MPI_Isend(bucketspercpu_sizes+j,1,MPI_INT,j,1,comm,requests+count); count++; //we use bucketspercpu_sizes because we need a permanent buffer for an asynchronous send
-                                                        for(k=0;k<buckets->Size();k++){
-                                                                Bucket<doubletype>* bucket=(Bucket<doubletype>*)buckets->GetObjectByOffset(k);
-                                                                bucket->Isend(j,requests,&count,comm);
-                                                        }
-                                                }
-                                        }
-                                }
-                                else{
-
-                                        /*Receive buckets from cpu i, and add them to my own my_rank bucket list: */
-                                        /*First, are we receiving anything from sender_rank? :*/
-                                        MPI_Recv(&size,1, MPI_INT,i,1, comm, &status);
-
-                                        /*If so, started receiving extra buckets and plug them into out buckets: */
-                                        if(size){
-                                                for(j=0;j<size;j++){
-                                                        Bucket<doubletype>* bucket=new Bucket<doubletype>();
-                                                        bucket->Recv(i,comm);
-                                                        bucketspercpu[my_rank]->AddObject(bucket);
-                                                }
-                                        }
-                                }
-                        }
-                        /*Wait for all requests to complete: */
-                        MPI_Waitall(num_requests,requests,statuses);
-
-                        /*Every cpu now has a dataset of buckets  in bucketspercpu[my_rank], which holds all the values 
-                         *local to this cpu that should be added to the global matrix. Just do that: */
-                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,comm);
-                        mybuckets=bucketspercpu[my_rank];
-
-                        for(i=0;i<mybuckets->Size();i++){
-                                Bucket<doubletype>* bucket=(Bucket<doubletype>*)mybuckets->GetObjectByOffset(i);
-                                bucket->SetLocalVectorValues(this->vector,lower_row);
-                        }
-
-                        /*Free ressources:{{{*/
-                        xDelete<int>(RowRank);
-                        for(i=0;i<num_procs;i++){
-                                DataSet* buckets=bucketspercpu[i];
-                                delete buckets;
-                        }
-                        xDelete<DataSet*>(bucketspercpu);
-                        xDelete<int>(bucketspercpu_sizes);
-                        xDelete<MPI_Request>(requests);
-                        /*}}}*/
-                }
-                /*}}}*/
                 /*FUNCTION SetValues{{{*/
                 void SetValues(int ssize, int* list, doubletype* values, InsMode mode){
@@ -571 +453 @@
                                         local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,this->vector[i]);
                                         MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
+                                        MPI_Bcast(&norm,1,MPI_DOUBLE,0,IssmComm::GetComm());
                                         return norm;
                                         break;
@@ -577 +460 @@
                                         for(i=0;i<this->m;i++)local_norm+=pow(this->vector[i],2);
                                         MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, 0, IssmComm::GetComm());
+                                        MPI_Bcast(&norm,1,MPI_DOUBLE,0,IssmComm::GetComm());
                                         return sqrt(norm);
                                         break;