00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "Field/BareField.h"
00028 #include "Field/BrickExpression.h"
00029 #include "FieldLayout/FieldLayout.h"
00030 #include "Message/Communicate.h"
00031 #include "Message/GlobalComm.h"
00032 #include "Message/Tags.h"
00033 #include "Utility/Inform.h"
00034 #include "Utility/Unique.h"
00035 #include "Utility/IpplInfo.h"
00036 #include "Utility/IpplStats.h"
00037 #include "Profile/Profiler.h"
00038
00039 #ifdef IPPL_STDSTL
00040 #include <map>
00041 #include <utility>
00042 using namespace std;
00043 #else
00044 #include <map.h>
00045 #include <multimap.h>
00046 #endif // IPPL_STDSTL
00047
00048 #include <stdlib.h>
00049
00050 #ifdef IPPL_NETCDF
00051 #include <netcdf.h>
00052 #endif
00053
00055
00056
00057
00058
00059 template< class T, unsigned Dim >
00060 BareField<T,Dim>::BareField(const BareField<T,Dim>& a)
00061 : Layout( a.Layout ),
00062 Gc( a.Gc ),
00063 compressible_m( a.compressible_m )
00064 {
00065 TAU_TYPE_STRING(taustr, "void (" + CT(a) + " )" );
00066 TAU_PROFILE("BareField::BareField()", taustr, TAU_FIELD);
00067
00068
00069 clearDirtyFlag();
00070
00071
00072 Locals_ac.reserve( a.size_if() );
00073
00074 for ( const_iterator_if a_i = a.begin_if(); a_i != a.end_if(); ++a_i )
00075 {
00076
00077 LField<T,Dim> *lf = new LField<T,Dim>( *((*a_i).second) );
00078
00079 #ifdef IPPL_SGI_CC_721_TYPENAME_BUG
00080 Locals_ac.insert( Locals_ac.end(),
00081 ac_id_larray::value_type((*a_i).first,lf) );
00082 #else
00083 Locals_ac.insert( Locals_ac.end(),
00084 typename ac_id_larray::value_type((*a_i).first,lf) );
00085 #endif // IPPL_SGI_CC_721_TYPENAME_BUG
00086 }
00087
00088 getLayout().checkin( *this , Gc );
00089
00090 }
00091
00092
00094
00095
00096
00097
00098 template< class T, unsigned Dim >
00099 BareField<T,Dim>::~BareField() {
00100 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
00101 TAU_PROFILE("BareField::~BareField()", taustr, TAU_FIELD);
00102
00103 if (Layout != 0) {
00104 Layout->checkout(*this);
00105 }
00106 }
00107
00108
00110
00111
00112
00113
00114 template< class T, unsigned Dim >
00115 void
00116 BareField<T,Dim>::initialize(Layout_t & l) {
00117 TAU_TYPE_STRING(taustr, "void (" + CT(l) + " )" );
00118 TAU_PROFILE("BareField::initialize()", taustr, TAU_FIELD);
00119
00120
00121 if (Layout == 0) {
00122 Layout = &l;
00123 setup();
00124 }
00125 }
00126
00127 template< class T, unsigned Dim >
00128 void
00129 BareField<T,Dim>::initialize(Layout_t & l,
00130 const GuardCellSizes<Dim>& gc) {
00131 TAU_TYPE_STRING(taustr, "void (" + CT(l) + ", " + CT(gc) + " )" );
00132 TAU_PROFILE("BareField::initialize()", taustr, TAU_FIELD);
00133
00134
00135 if (Layout == 0) {
00136 Layout = &l;
00137 Gc = gc;
00138 setup();
00139 }
00140 }
00141
00142
00144
00145
00146
00147
00148
00149 template< class T, unsigned Dim >
00150 void
00151 BareField<T,Dim>::setup()
00152 {
00153 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
00154 TAU_PROFILE("BareField::setup()", taustr, TAU_FIELD | TAU_ASSIGN);
00155
00156
00157
00158 if ( ! getLayout().fitsGuardCells(Gc)) {
00159 ERRORMSG("Cannot construct IPPL BareField:" << endl);
00160 ERRORMSG(" One or more vnodes is too small for the guard cells." << endl);
00161 ERRORMSG(" Guard cell sizes = " << Gc << endl);
00162 ERRORMSG(" FieldLayout = " << getLayout() << endl);
00163 Ippl::abort();
00164 }
00165
00166
00167 clearDirtyFlag();
00168
00169
00170 Locals_ac.reserve( getLayout().size_iv() );
00171
00172
00173 for (typename Layout_t::iterator_iv v_i=getLayout().begin_iv();
00174 v_i != getLayout().end_iv();
00175 ++v_i)
00176 {
00177
00178 const NDIndex<Dim> &owned = (*v_i).second->getDomain();
00179 NDIndex<Dim> guarded = AddGuardCells( owned , Gc );
00180
00181
00182 int vnode = (*v_i).second->getVnode();
00183
00184
00185 LField<T,Dim> *lf = new LField<T,Dim>( owned, guarded, vnode );
00186
00187
00188 #ifdef __MWERKS__
00189 Locals_ac.insert( end_if(),
00190 ac_id_larray::value_type((*v_i).first,lf));
00191 #else
00192 #ifdef IPPL_SGI_CC_721_TYPENAME_BUG
00193 Locals_ac.insert( end_if(),
00194 ac_id_larray::value_type((*v_i).first,lf));
00195 #else
00196 Locals_ac.insert( end_if(),
00197 typename ac_id_larray::value_type((*v_i).first,lf));
00198 #endif // IPPL_SGI_CC_721_TYPENAME_BUG
00199 #endif // __MWERKS__
00200 }
00201
00202 getLayout().checkin( *this , Gc );
00203
00204 }
00205
00207
00208
00209
00210
00211
00212 template< class T, unsigned Dim>
00213 void
00214 BareField<T,Dim>::write(ostream& out)
00215 {
00216 TAU_TYPE_STRING(taustr, CT(*this) + " void (ostream )" );
00217 TAU_PROFILE("BareField::write()", taustr, TAU_FIELD | TAU_IO);
00218
00219
00220
00221
00222
00223 int tag = Ippl::Comm->next_tag(F_WRITE_TAG, F_TAG_CYCLE);
00224 if (Ippl::myNode() != 0) {
00225 for ( iterator_if local = begin_if(); local != end_if() ; ++local) {
00226
00227 LField<T,Dim>& l = *((*local).second);
00228 NDIndex<Dim>& lo = (NDIndex<Dim>&) l.getOwned();
00229 typename LField<T,Dim>::iterator rhs(l.begin());
00230
00231
00232 if (Ippl::myNode() != 0) {
00233 Message *mess = new Message();
00234 lo.putMessage(*mess);
00235 rhs.putMessage(*mess);
00236
00237 Ippl::Comm->send(mess, 0, tag);
00238 }
00239 }
00240 } else {
00241
00242 LField<T,Dim> data(getDomain(), getDomain());
00243 data.Uncompress();
00244
00245
00246 for ( iterator_if local = begin_if(); local != end_if() ; ++local) {
00247
00248 LField<T,Dim>& l = *((*local).second);
00249 NDIndex<Dim>& lo = (NDIndex<Dim>&) l.getOwned();
00250 typename LField<T,Dim>::iterator rhs(l.begin());
00251
00252
00253
00254
00255 typename LField<T,Dim>::iterator putloc = data.begin(lo);
00256 typename LField<T,Dim>::iterator getloc = l.begin(lo);
00257 for ( ; getloc != l.end() ; ++putloc, ++getloc ) {
00258
00259
00260 *putloc = *getloc;
00261 }
00262 }
00263
00264
00265 int remaining = getLayout().size_rdv();
00266
00267
00268 for ( ; remaining > 0; --remaining) {
00269
00270 int any_node = COMM_ANY_NODE;
00271 Message *mess = Ippl::Comm->receive_block(any_node, tag);
00272
00273
00274 NDIndex<Dim> lo;
00275 T compressed_data;
00276 typename LField<T,Dim>::iterator rhs(compressed_data);
00277 lo.getMessage(*mess);
00278 rhs.getMessage(*mess);
00279
00280
00281
00282 typename LField<T,Dim>::iterator putloc = data.begin(lo);
00283 for (unsigned elems=lo.size(); elems > 0; ++putloc, ++rhs, --elems)
00284 *putloc = *rhs;
00285
00286
00287 delete mess;
00288 }
00289
00290
00291 out << data;
00292 }
00293 }
00294
00295
00297
00298
00299
00300
00301
00302 template< class T, unsigned Dim >
00303 void BareField<T,Dim>::fillGuardCells(bool reallyFill) const
00304 {
00305
00306 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
00307 TAU_PROFILE("BareField::fillGuardCells()", taustr, TAU_FIELD);
00308
00309 TAU_PROFILE_TIMER(sendtimer, " fillGuardCells-send",
00310 taustr, TAU_FIELD);
00311 TAU_PROFILE_TIMER(sendcommtimer, " fillGuardCells-send-comm",
00312 taustr, TAU_FIELD);
00313 TAU_PROFILE_TIMER(findrectimer, " fillGuardCells-findreceive",
00314 taustr, TAU_FIELD);
00315 TAU_PROFILE_TIMER(localstimer, " fillGuardCells-locals",
00316 taustr, TAU_FIELD);
00317 TAU_PROFILE_TIMER(rectimer, " fillGuardCells-receive",
00318 taustr, TAU_FIELD);
00319 TAU_PROFILE_TIMER(reccommtimer, " fillGuardCells-receive-comm",
00320 taustr, TAU_FIELD);
00321 TAU_PROFILE_TIMER(localsexprtimer, " fillGuardCells-locals-expression",
00322 taustr, TAU_FIELD);
00323
00324
00325
00326 BareField<T,Dim>& ncf = const_cast<BareField<T,Dim>&>(*this);
00327
00328
00329 if (!reallyFill)
00330 return;
00331
00332
00333 ncf.clearDirtyFlag();
00334
00335
00336 if (Gc == GuardCellSizes<Dim>())
00337 return;
00338
00339
00340
00341
00342
00343 iterator_if lf_i, lf_e = ncf.end_if();
00344
00345 #ifdef IPPL_PRINTDEBUG
00346 Inform msg("fillGuardCells", INFORM_ALL_NODES);
00347 #endif
00348
00349
00350
00351
00352
00353 int nprocs = Ippl::getNodes();
00354 if (nprocs > 1) {
00355 TAU_PROFILE_START(findrectimer);
00356
00357 typedef multimap< NDIndex<Dim> , LField<T,Dim>* , less<NDIndex<Dim> > > ac_recv_type;
00358 ac_recv_type recv_ac;
00359 bool* recvmsg = new bool[nprocs];
00360 TAU_PROFILE_STOP(findrectimer);
00361
00362 TAU_PROFILE_START(sendtimer);
00363
00364 Message** mess = new Message*[nprocs];
00365 #ifdef IPPL_PRINTDEBUG
00366 int* ndomains = new int[nprocs];
00367 #endif
00368 int iproc;
00369 for (iproc=0; iproc<nprocs; ++iproc) {
00370 mess[iproc] = NULL;
00371 recvmsg[iproc] = false;
00372 #ifdef IPPL_PRINTDEBUG
00373 ndomains[iproc] = 0;
00374 #endif
00375 }
00376 TAU_PROFILE_STOP(sendtimer);
00377
00378 for (lf_i = ncf.begin_if(); lf_i != lf_e; ++lf_i) {
00379 TAU_PROFILE_START(findrectimer);
00380
00381 LField<T,Dim> &lf = *((*lf_i).second);
00382 const NDIndex<Dim> &lf_domain = lf.getAllocated();
00383
00384 typename Layout_t::touch_range_dv
00385 rrange(ncf.getLayout().touch_range_rdv(lf_domain));
00386 typename Layout_t::touch_iterator_dv rv_i;
00387 for (rv_i = rrange.first; rv_i != rrange.second; ++rv_i) {
00388
00389 NDIndex<Dim> sub = lf_domain.intersect( (*rv_i).first );
00390 typedef typename ac_recv_type::value_type value_type;
00391 recv_ac.insert(value_type(sub,&lf));
00392
00393 int rnode = (*rv_i).second->getNode();
00394 recvmsg[rnode] = true;
00395 }
00396 TAU_PROFILE_STOP(findrectimer);
00397
00398 TAU_PROFILE_START(sendtimer);
00399 const NDIndex<Dim>& lo = lf.getOwned();
00400 #ifdef IPPL_PRINTDEBUG
00401 msg << "Finding send overlap regions for domain " << lo << endl;
00402 #endif
00403
00404
00405 typename Layout_t::touch_range_dv
00406 range( ncf.getLayout().touch_range_rdv(lo,ncf.getGC()) );
00407 typename Layout_t::touch_iterator_dv remote_i;
00408 for (remote_i = range.first; remote_i != range.second; ++remote_i) {
00409
00410 NDIndex<Dim> intersection = lo.intersect( (*remote_i).first );
00411
00412 int rnode = (*remote_i).second->getNode();
00413 #ifdef IPPL_PRINTDEBUG
00414 msg << " Overlap domain = " << (*remote_i).first << endl;
00415 msg << " Inters. domain = " << intersection;
00416 msg << " --> node " << rnode << endl;
00417 #endif
00418
00419
00420
00421 T compressed_value;
00422 LFI msgval = lf.begin(intersection, compressed_value);
00423 msgval.TryCompress();
00424
00425
00426 if (!mess[rnode]) mess[rnode] = new Message;
00427 PAssert(mess[rnode]);
00428 mess[rnode]->put(intersection);
00429 mess[rnode]->put(msgval);
00430 #ifdef IPPL_PRINTDEBUG
00431 ndomains[rnode]++;
00432 #endif
00433 }
00434 TAU_PROFILE_STOP(sendtimer);
00435 }
00436 TAU_PROFILE_START(findrectimer);
00437 int remaining = 0;
00438 for (iproc=0; iproc<nprocs; ++iproc)
00439 if (recvmsg[iproc]) ++remaining;
00440 delete [] recvmsg;
00441 TAU_PROFILE_STOP(findrectimer);
00442
00443 TAU_PROFILE_START(sendtimer);
00444
00445 int tag = Ippl::Comm->next_tag( F_GUARD_CELLS_TAG , F_TAG_CYCLE );
00446 TAU_PROFILE_START(sendcommtimer);
00447
00448 for (iproc=0; iproc<nprocs; ++iproc) {
00449 if (mess[iproc]) {
00450 #ifdef IPPL_PRINTDEBUG
00451 msg << "fillGuardCells: Sending message to node " << iproc << endl
00452 << " number of domains = " << ndomains[iproc] << endl
00453 << " number of MsgItems = " << mess[iproc]->size() << endl;
00454 #endif
00455 Ippl::Comm->send(mess[iproc], iproc, tag);
00456 }
00457 }
00458 TAU_PROFILE_STOP(sendcommtimer);
00459 delete [] mess;
00460 #ifdef IPPL_PRINTDEBUG
00461 delete [] ndomains;
00462 #endif
00463 TAU_PROFILE_STOP(sendtimer);
00464
00465
00466
00467
00468 TAU_PROFILE_START(localstimer);
00469 for (lf_i = ncf.begin_if(); lf_i != lf_e; ++lf_i)
00470 {
00471
00472 LField<T,Dim> &lf = *(*lf_i).second;
00473 const NDIndex<Dim>& la = lf.getAllocated();
00474
00475
00476
00477
00478
00479
00480
00481 if (!lf.OverlapCacheInitialized())
00482 {
00483 for (iterator_if rf_i = ncf.begin_if(); rf_i != ncf.end_if(); ++rf_i)
00484 if ( rf_i != lf_i )
00485 {
00486
00487 LField<T,Dim> &rf = *(*rf_i).second;
00488 const NDIndex<Dim>& ro = rf.getOwned();
00489
00490
00491 if ( la.touches(ro) )
00492 lf.AddToOverlapCache(&rf);
00493 }
00494 }
00495
00496
00497
00498
00499 for (typename LField<T,Dim>::OverlapIterator rf_i = lf.BeginOverlap();
00500 rf_i != lf.EndOverlap(); ++rf_i)
00501 {
00502 LField<T, Dim> &rf = *(*rf_i);
00503 const NDIndex<Dim>& ro = rf.getOwned();
00504
00505 bool c1 = lf.IsCompressed();
00506 bool c2 = rf.IsCompressed();
00507 bool c3 = *rf.begin() == *lf.begin();
00508
00509
00510 if ( !( c1 && c2 && c3 ) )
00511 {
00512 TAU_PROFILE_START(localsexprtimer);
00513
00514
00515 NDIndex<Dim> intersection = la.intersect(ro);
00516
00517
00518 LFI rhs = rf.begin(intersection);
00519
00520
00521 if ( !(c1 && rhs.CanCompress(*lf.begin())) )
00522 {
00523
00524 lf.Uncompress();
00525
00526
00527 LFI lhs = lf.begin(intersection);
00528 #ifdef IPPL_PRINTDEBUG
00529 msg << " Inters. domain=" << intersection << endl;
00530 #endif
00531
00532 BrickExpression<Dim,LFI,LFI,OpAssign>(lhs,rhs).apply();
00533 }
00534
00535 TAU_PROFILE_STOP(localsexprtimer);
00536 }
00537 #ifdef IPPL_PRINTDEBUG
00538 else {
00539 msg << " Both sides compressed and equal ... val = ";
00540 msg << *rf.begin() << endl;
00541 }
00542 #endif
00543 }
00544 }
00545 TAU_PROFILE_STOP(localstimer);
00546
00547
00548
00549 TAU_PROFILE_START(rectimer);
00550 while (remaining>0) {
00551
00552 int any_node = COMM_ANY_NODE;
00553 TAU_PROFILE_START(reccommtimer);
00554 Message* rmess = Ippl::Comm->receive_block(any_node,tag);
00555 PAssert(rmess);
00556 --remaining;
00557 TAU_PROFILE_STOP(reccommtimer);
00558
00559
00560 int ndoms = rmess->size() / (Dim + 3);
00561 #ifdef IPPL_PRINTDEBUG
00562 msg << "fillGuardCells: Message received from node " << any_node
00563 << ", number of domains = " << ndoms << endl;
00564 #endif
00565 for (int idom=0; idom<ndoms; ++idom) {
00566
00567 NDIndex<Dim> intersection;
00568 intersection.getMessage(*rmess);
00569
00570
00571 T compressed_value;
00572 LFI rhs(compressed_value);
00573 rhs.getMessage(*rmess);
00574 #ifdef IPPL_PRINTDEBUG
00575 msg << "Received remote overlap region = " << intersection << endl;
00576 #endif
00577
00578
00579 typename ac_recv_type::iterator hit = recv_ac.find( intersection );
00580 PAssert( hit != recv_ac.end() );
00581
00582 LField<T,Dim> &lf = *(*hit).second;
00583
00584 #ifdef IPPL_PRINTDEBUG
00585 msg << " LHS compressed ? " << lf.IsCompressed();
00586 msg << ", LHS value = " << *lf.begin() << endl;
00587 msg << " RHS compressed ? " << rhs.IsCompressed();
00588 msg << ", RHS value = " << *rhs << endl;
00589 msg << " *rhs == *lf.begin() ? " << (*rhs == *lf.begin()) << endl;
00590 #endif
00591 if ( !(rhs.IsCompressed() && lf.IsCompressed() &&
00592 (*rhs == *lf.begin())) )
00593 {
00594
00595 lf.Uncompress();
00596 LFI lhs = lf.begin(intersection);
00597
00598 #ifdef IPPL_PRINTDEBUG
00599 msg << " Doing copy." << endl;
00600 #endif
00601 BrickExpression<Dim,LFI,LFI,OpAssign>(lhs,rhs).apply();
00602 }
00603
00604
00605 recv_ac.erase( hit );
00606 }
00607 delete rmess;
00608 }
00609 TAU_PROFILE_STOP(rectimer);
00610 }
00611 else {
00612
00613
00614
00615 TAU_PROFILE_START(localstimer);
00616 for (lf_i = ncf.begin_if(); lf_i != lf_e; ++lf_i)
00617 {
00618
00619 LField<T,Dim> &lf = *(*lf_i).second;
00620 const NDIndex<Dim>& la = lf.getAllocated();
00621
00622
00623
00624
00625
00626
00627
00628 if (!lf.OverlapCacheInitialized())
00629 {
00630 for (iterator_if rf_i = ncf.begin_if(); rf_i != ncf.end_if(); ++rf_i)
00631 if ( rf_i != lf_i )
00632 {
00633
00634 LField<T,Dim> &rf = *(*rf_i).second;
00635 const NDIndex<Dim>& ro = rf.getOwned();
00636
00637
00638 if ( la.touches(ro) )
00639 lf.AddToOverlapCache(&rf);
00640 }
00641 }
00642
00643
00644
00645
00646 for (typename LField<T,Dim>::OverlapIterator rf_i = lf.BeginOverlap();
00647 rf_i != lf.EndOverlap(); ++rf_i)
00648 {
00649 LField<T, Dim> &rf = *(*rf_i);
00650 const NDIndex<Dim>& ro = rf.getOwned();
00651
00652 bool c1 = lf.IsCompressed();
00653 bool c2 = rf.IsCompressed();
00654 bool c3 = *rf.begin() == *lf.begin();
00655
00656
00657 if ( !( c1 && c2 && c3 ) )
00658 {
00659 TAU_PROFILE_START(localsexprtimer);
00660
00661
00662 NDIndex<Dim> intersection = la.intersect(ro);
00663
00664
00665 LFI rhs = rf.begin(intersection);
00666
00667
00668 if ( !(c1 && rhs.CanCompress(*lf.begin())) )
00669 {
00670
00671 lf.Uncompress();
00672
00673
00674 LFI lhs = lf.begin(intersection);
00675 #ifdef IPPL_PRINTDEBUG
00676 msg << " Inters. domain=" << intersection << endl;
00677 #endif
00678
00679 BrickExpression<Dim,LFI,LFI,OpAssign>(lhs,rhs).apply();
00680 }
00681
00682 TAU_PROFILE_STOP(localsexprtimer);
00683 }
00684 #ifdef IPPL_PRINTDEBUG
00685 else {
00686 msg << " Both sides compressed and equal ... val = ";
00687 msg << *rf.begin() << endl;
00688 }
00689 #endif
00690 }
00691 }
00692 TAU_PROFILE_STOP(localstimer);
00693 }
00694 return;
00695 }
00696
00698
00699
00700
00701
00702
00703
00704 template <class T, unsigned Dim>
00705 void BareField<T,Dim>::setGuardCells(const T& val) const
00706 {
00707
00708 TAU_TYPE_STRING(taustr, CT(*this) + " void (" + CT(val) + ")" );
00709 TAU_PROFILE("BareField::setGuardCells()", taustr, TAU_FIELD);
00710
00711
00712 if (Gc == GuardCellSizes<Dim>())
00713 return;
00714
00715
00716 BareField<T,Dim>& ncf = const_cast<BareField<T,Dim>&>(*this);
00717
00718
00719 iterator_if lf_i, lf_e = ncf.end_if();
00720 for (lf_i = ncf.begin_if(); lf_i != lf_e; ++lf_i) {
00721
00722 LField<T,Dim>& lf = *((*lf_i).second);
00723
00724
00725
00726
00727
00728 if (lf.IsCompressed() && lf.getCompressedData() == val)
00729 continue;
00730
00731
00732
00733 const NDIndex<Dim>& adom = lf.getAllocated();
00734 const NDIndex<Dim>& odom = lf.getOwned();
00735 NDIndex<Dim> dom = odom;
00736
00737
00738
00739 LField<T,Dim> rf(odom,adom);
00740 rf.Compress(val);
00741
00742
00743 lf.Uncompress();
00744 LFI lhs, rhs;
00745
00746
00747 for (int idim = 0; idim < Dim; ++idim) {
00748
00749 dom[idim] = Index(adom[idim].first(),
00750 adom[idim].first() + Gc.left(idim) - 1);
00751 if (!dom[idim].empty()) {
00752
00753 lhs = lf.begin(dom);
00754 rhs = rf.begin(dom);
00755 BrickExpression<Dim,LFI,LFI,OpAssign>(lhs,rhs).apply();
00756 }
00757
00758
00759 dom[idim] = Index(adom[idim].last() - Gc.right(idim) + 1,
00760 adom[idim].last());
00761 if (!dom[idim].empty()) {
00762
00763 lhs = lf.begin(dom);
00764 rhs = rf.begin(dom);
00765 BrickExpression<Dim,LFI,LFI,OpAssign>(lhs,rhs).apply();
00766 }
00767
00768
00769
00770 dom[idim] = adom[idim];
00771 }
00772 }
00773
00774
00775 ncf.setDirtyFlag();
00776
00777 return;
00778 }
00779
00781
00782
00783
00784
00785
00786 template <class T, unsigned Dim>
00787 void BareField<T,Dim>::accumGuardCells()
00788 {
00789
00790 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
00791 TAU_PROFILE("BareField::accumGuardCells()", taustr, TAU_FIELD);
00792
00793 TAU_PROFILE_TIMER(sendtimer, " accumGuardCells-send",
00794 taustr, TAU_FIELD);
00795 TAU_PROFILE_TIMER(sendcommtimer, " accumGuardCells-send-comm",
00796 taustr, TAU_FIELD);
00797 TAU_PROFILE_TIMER(findrectimer, " accumGuardCells-findreceive",
00798 taustr, TAU_FIELD);
00799 TAU_PROFILE_TIMER(localstimer, " accumGuardCells-locals",
00800 taustr, TAU_FIELD);
00801 TAU_PROFILE_TIMER(rectimer, " accumGuardCells-receive",
00802 taustr, TAU_FIELD);
00803 TAU_PROFILE_TIMER(reccommtimer, " accumGuardCells-receive-comm",
00804 taustr, TAU_FIELD);
00805 TAU_PROFILE_TIMER(localsexprtimer, " accumGuardCells-locals-expression",
00806 taustr, TAU_FIELD);
00807
00808
00809 if (Gc == GuardCellSizes<Dim>())
00810 return;
00811
00812
00813 iterator_if lf_i, lf_e = end_if();
00814
00815 #ifdef IPPL_PRINTDEBUG
00816 Inform msg("accumGuardCells", INFORM_ALL_NODES);
00817 #endif
00818
00819
00820
00821
00822
00823 int nprocs = Ippl::getNodes();
00824 if (nprocs > 1) {
00825 TAU_PROFILE_START(findrectimer);
00826
00827 typedef multimap< NDIndex<Dim> , LField<T,Dim>* , less<NDIndex<Dim> > > ac_recv_type;
00828 ac_recv_type recv_ac;
00829 bool* recvmsg = new bool[nprocs];
00830 TAU_PROFILE_STOP(findrectimer);
00831
00832 TAU_PROFILE_START(sendtimer);
00833
00834 Message** mess = new Message*[nprocs];
00835 #ifdef IPPL_PRINTDEBUG
00836 int* ndomains = new int[nprocs];
00837 #endif
00838 int iproc;
00839 for (iproc=0; iproc<nprocs; ++iproc) {
00840 mess[iproc] = NULL;
00841 recvmsg[iproc] = false;
00842 #ifdef IPPL_PRINTDEBUG
00843 ndomains[iproc] = 0;
00844 #endif
00845 }
00846 TAU_PROFILE_STOP(sendtimer);
00847
00848 for (lf_i = begin_if(); lf_i != lf_e; ++lf_i) {
00849 TAU_PROFILE_START(findrectimer);
00850
00851 LField<T,Dim> &lf = *((*lf_i).second);
00852 const NDIndex<Dim>& lo = lf.getOwned();
00853
00854 typename Layout_t::touch_range_dv
00855 rrange( getLayout().touch_range_rdv(lo,Gc) );
00856 typename Layout_t::touch_iterator_dv rv_i;
00857 for (rv_i = rrange.first; rv_i != rrange.second; ++rv_i) {
00858
00859 NDIndex<Dim> sub = lo.intersect( (*rv_i).first );
00860 recv_ac.insert(ac_recv_type::value_type(sub,&lf));
00861
00862 int rnode = (*rv_i).second->getNode();
00863 recvmsg[rnode] = true;
00864 }
00865 TAU_PROFILE_STOP(findrectimer);
00866
00867 TAU_PROFILE_START(sendtimer);
00868 const NDIndex<Dim> &lf_domain = lf.getAllocated();
00869 #ifdef IPPL_PRINTDEBUG
00870 msg << "Finding send overlap regions for domain " << lf_domain << endl;
00871 #endif
00872
00873
00874 typename Layout_t::touch_range_dv
00875 range(getLayout().touch_range_rdv(lf_domain));
00876 typename Layout_t::touch_iterator_dv remote_i;
00877 for (remote_i = range.first; remote_i != range.second; ++remote_i) {
00878
00879 NDIndex<Dim> intersection = lf_domain.intersect( (*remote_i).first );
00880
00881 int rnode = (*remote_i).second->getNode();
00882 #ifdef IPPL_PRINTDEBUG
00883 msg << " Overlap domain = " << (*remote_i).first << endl;
00884 msg << " Inters. domain = " << intersection;
00885 msg << " --> node " << rnode << endl;
00886 #endif
00887
00888
00889
00890 T compressed_value;
00891 LFI msgval = lf.begin(intersection, compressed_value);
00892 msgval.TryCompress();
00893
00894
00895 if (!mess[rnode]) mess[rnode] = new Message;
00896 PAssert(mess[rnode]);
00897 mess[rnode]->put(intersection);
00898 mess[rnode]->put(msgval);
00899 #ifdef IPPL_PRINTDEBUG
00900 ndomains[rnode]++;
00901 #endif
00902 }
00903 TAU_PROFILE_STOP(sendtimer);
00904 }
00905 TAU_PROFILE_START(findrectimer);
00906 int remaining = 0;
00907 for (iproc=0; iproc<nprocs; ++iproc)
00908 if (recvmsg[iproc]) ++remaining;
00909 delete [] recvmsg;
00910 TAU_PROFILE_STOP(findrectimer);
00911
00912 TAU_PROFILE_START(sendtimer);
00913
00914 int tag = Ippl::Comm->next_tag( F_GUARD_CELLS_TAG , F_TAG_CYCLE );
00915 TAU_PROFILE_START(sendcommtimer);
00916
00917 for (iproc=0; iproc<nprocs; ++iproc) {
00918 if (mess[iproc]) {
00919 #ifdef IPPL_PRINTDEBUG
00920 msg << "accumGuardCells: Sending message to node " << iproc << endl
00921 << " number of domains = " << ndomains[iproc] << endl
00922 << " number of MsgItems = " << mess[iproc]->size() << endl;
00923 #endif
00924 Ippl::Comm->send(mess[iproc], iproc, tag);
00925 }
00926 }
00927 TAU_PROFILE_STOP(sendcommtimer);
00928 delete [] mess;
00929 #ifdef IPPL_PRINTDEBUG
00930 delete [] ndomains;
00931 #endif
00932 TAU_PROFILE_STOP(sendtimer);
00933
00934
00935
00936
00937 TAU_PROFILE_START(localstimer);
00938 for (lf_i=begin_if(); lf_i != lf_e; ++lf_i)
00939 {
00940
00941 LField<T,Dim> &lf = *(*lf_i).second;
00942 const NDIndex<Dim>& la = lf.getAllocated();
00943
00944
00945
00946
00947
00948
00949
00950 if (!lf.OverlapCacheInitialized())
00951 {
00952 for (iterator_if rf_i = begin_if(); rf_i != end_if(); ++rf_i)
00953 if ( rf_i != lf_i )
00954 {
00955
00956 LField<T,Dim> &rf = *(*rf_i).second;
00957 const NDIndex<Dim>& ro = rf.getOwned();
00958
00959
00960 if ( la.touches(ro) )
00961 lf.AddToOverlapCache(&rf);
00962 }
00963 }
00964
00965
00966
00967
00968 for (typename LField<T,Dim>::OverlapIterator rf_i = lf.BeginOverlap();
00969 rf_i != lf.EndOverlap(); ++rf_i)
00970 {
00971 LField<T, Dim> &rf = *(*rf_i);
00972 const NDIndex<Dim>& ro = rf.getOwned();
00973
00974 TAU_PROFILE_START(localsexprtimer);
00975
00976
00977 NDIndex<Dim> intersection = la.intersect(ro);
00978
00979
00980 LFI lhs = lf.begin(intersection);
00981
00982
00983 if ( !lhs.CanCompress(T()) ) {
00984
00985
00986 rf.Uncompress();
00987
00988
00989 LFI rhs = rf.begin(intersection);
00990
00991 #ifdef IPPL_PRINTDEBUG
00992 msg << " Inters. domain=" << intersection << endl;
00993 #endif
00994
00995 BrickExpression<Dim,LFI,LFI,OpAddAssign>(rhs,lhs).apply();
00996 }
00997 TAU_PROFILE_STOP(localsexprtimer);
00998 }
00999 }
01000 TAU_PROFILE_STOP(localstimer);
01001
01002
01003
01004 TAU_PROFILE_START(rectimer);
01005 while (remaining>0) {
01006
01007 int any_node = COMM_ANY_NODE;
01008 TAU_PROFILE_START(reccommtimer);
01009 Message* rmess = Ippl::Comm->receive_block(any_node,tag);
01010 PAssert(rmess);
01011 --remaining;
01012 TAU_PROFILE_STOP(reccommtimer);
01013
01014
01015 int ndoms = rmess->size() / (Dim + 3);
01016 #ifdef IPPL_PRINTDEBUG
01017 msg << "accumGuardCells: Message received from node " << any_node
01018 << ", number of domains = " << ndoms << endl;
01019 #endif
01020 for (int idom=0; idom<ndoms; ++idom) {
01021
01022 NDIndex<Dim> intersection;
01023 intersection.getMessage(*rmess);
01024
01025
01026 T compressed_value;
01027 LFI rhs(compressed_value);
01028 rhs.getMessage(*rmess);
01029 #ifdef IPPL_PRINTDEBUG
01030 msg << "Received remote overlap region = " << intersection << endl;
01031 #endif
01032
01033
01034 typename ac_recv_type::iterator hit = recv_ac.find( intersection );
01035 PAssert( hit != recv_ac.end() );
01036
01037
01038 LField<T,Dim> &lf = *(*hit).second;
01039
01040
01041 lf.Uncompress();
01042 LFI lhs = lf.begin(intersection);
01043
01044
01045 BrickExpression<Dim,LFI,LFI,OpAddAssign>(lhs,rhs).apply();
01046
01047
01048 recv_ac.erase( hit );
01049 }
01050 delete rmess;
01051 }
01052 TAU_PROFILE_STOP(rectimer);
01053 }
01054 else {
01055
01056
01057
01058 TAU_PROFILE_START(localstimer);
01059 for (lf_i=begin_if(); lf_i != lf_e; ++lf_i)
01060 {
01061
01062 LField<T,Dim> &lf = *(*lf_i).second;
01063 const NDIndex<Dim>& la = lf.getAllocated();
01064
01065
01066
01067
01068
01069
01070
01071 if (!lf.OverlapCacheInitialized())
01072 {
01073 for (iterator_if rf_i = begin_if(); rf_i != end_if(); ++rf_i)
01074 if ( rf_i != lf_i )
01075 {
01076
01077 LField<T,Dim> &rf = *(*rf_i).second;
01078 const NDIndex<Dim>& ro = rf.getOwned();
01079
01080
01081 if ( la.touches(ro) )
01082 lf.AddToOverlapCache(&rf);
01083 }
01084 }
01085
01086
01087
01088
01089 for (typename LField<T,Dim>::OverlapIterator rf_i = lf.BeginOverlap();
01090 rf_i != lf.EndOverlap(); ++rf_i)
01091 {
01092 LField<T, Dim> &rf = *(*rf_i);
01093 const NDIndex<Dim>& ro = rf.getOwned();
01094
01095 TAU_PROFILE_START(localsexprtimer);
01096
01097
01098 NDIndex<Dim> intersection = la.intersect(ro);
01099
01100
01101 LFI lhs = lf.begin(intersection);
01102
01103
01104 if ( !lhs.CanCompress(T()) ) {
01105
01106 rf.Uncompress();
01107
01108
01109 LFI rhs = rf.begin(intersection);
01110
01111 #ifdef IPPL_PRINTDEBUG
01112 msg << " Inters. domain=" << intersection << endl;
01113 #endif
01114
01115 BrickExpression<Dim,LFI,LFI,OpAddAssign>(rhs,lhs).apply();
01116 }
01117 TAU_PROFILE_STOP(localsexprtimer);
01118 }
01119 }
01120 TAU_PROFILE_STOP(localstimer);
01121 }
01122
01123
01124
01125 setDirtyFlag();
01126
01127 fillGuardCellsIfNotDirty();
01128
01129 Compress();
01130
01131 return;
01132 }
01133
01135
01136
01137
01138
01139
01140 template< class T, unsigned Dim >
01141 void BareField<T,Dim>::write(char* fname) const
01142 {
01143 TAU_TYPE_STRING(taustr, CT(*this) + " void (char * )" );
01144 TAU_PROFILE("BareField::write()", taustr, TAU_FIELD | TAU_IO);
01145
01146 #ifdef IPPL_NETCDF
01147
01148 int ncid;
01149 int varVID;
01150
01151
01152 ncid = nccreate(fname, NC_CLOBBER);
01153
01154
01155 int fillmode = ncsetfill(ncid,NC_NOFILL);
01156
01157
01158
01159
01160 const NDIndex<Dim> &domain = getDomain();
01161
01162
01163
01164
01165
01166 int dimids[Dim];
01167 if ( Dim==1 )
01168 {
01169 dimids[0] = ncdimdef(ncid, "nx", domain[0].length());
01170 }
01171 else if ( Dim==2 )
01172 {
01173 dimids[1] = ncdimdef(ncid, "nx", domain[1].length());
01174 dimids[0] = ncdimdef(ncid, "ny", domain[0].length());
01175 }
01176 else if ( Dim==3 )
01177 {
01178 dimids[2] = ncdimdef(ncid, "nx", domain[0].length());
01179 dimids[1] = ncdimdef(ncid, "ny", domain[1].length());
01180 dimids[0] = ncdimdef(ncid, "nz", domain[2].length());
01181 }
01182 else if ( Dim>3 )
01183 {
01184 ERRORMSG("BareField: Can't write more than 3 dimensions." << endl);
01185 }
01186 varVID = ncvardef(ncid, fname, NC_DOUBLE, Dim, dimids);
01187
01188 int errRet = ncendef(ncid);
01189
01190
01191 long startIndex[Dim];
01192 long countIndex[Dim];
01193 int rcode;
01194
01195 for (const_iterator_if l_i = begin_if(); l_i != end_if(); ++l_i)
01196 {
01197
01198 LField<T,Dim> *ldf = (*l_i).second.get();
01199 const NDIndex<Dim> &Owned = ldf->getOwned();
01200
01201 int size = 1;
01202 for( int i = 0 ; i < Dim ; i++ )
01203 {
01204 startIndex[Dim - i - 1] = Owned[i].first();
01205 countIndex[Dim - i - 1] = Owned[i].length();
01206 size *= countIndex[Dim - i - 1];
01207 }
01208 double* buffer = new double[size];
01209
01210 int icount = 0;
01211 typename LField<T,Dim>::iterator lf_bi = ldf->begin();
01212 if ( Dim==1 )
01213 {
01214 int n0 = ldf->size(0);
01215 for (int i0=0; i0<n0; ++i0)
01216 buffer[icount++] = lf_bi.offset(i0);
01217 }
01218 else if ( Dim==2 )
01219 {
01220 int n0 = ldf->size(0);
01221 int n1 = ldf->size(1);
01222 for (int i1=0; i1<n1; ++i1)
01223 for (int i0=0; i0<n0; ++i0)
01224 buffer[icount++] = lf_bi.offset(i0,i1);
01225 }
01226 else if ( Dim==3 )
01227 {
01228 int n0 = ldf->size(0);
01229 int n1 = ldf->size(1);
01230 int n2 = ldf->size(2);
01231 for (int i2=0; i2<n2; ++i2)
01232 for (int i1=0; i1<n1; ++i1)
01233 for (int i0=0; i0<n0; ++i0)
01234 buffer[icount++] = lf_bi.offset(i0,i1,i2);
01235 }
01236 rcode = ncvarput(ncid, varVID, startIndex,countIndex, buffer);
01237 if ( rcode != 0)
01238 {
01239 ERRORMSG("BareField: ncvarput() error, rcode=" << rcode << endl);
01240 }
01241 delete [] buffer;
01242 }
01243 rcode = ncclose(ncid);
01244 #else
01245 ERRORMSG("[Bare]Field::write(\"filename\") requires that you run 'conf' "
01246 << endl << " "
01247 << "with the NETCDF option when building libippl.a; you haven't."
01248 << endl);
01249 #endif
01250 }
01251
01253
01254
01255
01256
01257
01258
01259
01260 template< class T, unsigned Dim >
01261 void BareField<T,Dim>::writeb(char* fname) const
01262 {
01263 TAU_TYPE_STRING(taustr, CT(*this) + " void (char * )" );
01264 TAU_PROFILE("BareField::writeb()", taustr, TAU_FIELD | TAU_IO);
01265 Inform outC(0, fname, Inform::OVERWRITE);
01266
01267 int icount = 0;
01268 for (const_iterator_if l_i = begin_if(); l_i != end_if(); ++l_i)
01269 {
01270
01271 LField<T,Dim> *ldf = (*l_i).second.get();
01272 const NDIndex<Dim> &Owned = ldf->getOwned();
01273
01274 outC << "vnode = " << icount++ << endl;
01275 typename LField<T,Dim>::iterator lf_bi = ldf->begin();
01276 if ( Dim==1 )
01277 {
01278 int n0 = ldf->size(0);
01279 int l0 = -Gc.left(0);
01280 int r0 = n0 + Gc.right(0);
01281 for (int i0=l0; i0<r0; ++i0)
01282 {
01283 outC << " [" << i0;
01284 outC << "]=" << lf_bi.offset(i0);
01285 }
01286 }
01287 else if ( Dim==2 )
01288 {
01289 int n0 = ldf->size(0);
01290 int n1 = ldf->size(1);
01291 int l0 = -Gc.left(0);
01292 int l1 = -Gc.left(1);
01293 int r0 = n0 + Gc.right(0);
01294 int r1 = n1 + Gc.right(1);
01295 for (int i1=l1; i1<r1; ++i1)
01296 {
01297 for (int i0=l0; i0<r0; ++i0)
01298 {
01299 outC << " [" << i0;
01300 outC << "][" << i1;
01301 outC << "]=" << lf_bi.offset(i0,i1);
01302 }
01303 outC << endl;
01304 }
01305 }
01306 else if ( Dim==3 )
01307 {
01308 int n0 = ldf->size(0);
01309 int n1 = ldf->size(1);
01310 int n2 = ldf->size(2);
01311 int l0 = -Gc.left(0);
01312 int l1 = -Gc.left(1);
01313 int l2 = -Gc.left(2);
01314 int r0 = n0 + Gc.right(0);
01315 int r1 = n1 + Gc.right(1);
01316 int r2 = n2 + Gc.right(2);
01317 for (int i2=l2; i2<r2; ++i2)
01318 {
01319 for (int i1=l1; i1<r1; ++i1)
01320 for (int i0=l0; i0<r0; ++i0)
01321 {
01322 outC << " [" << i0;
01323 outC << "][" << i1;
01324 outC << "][" << i2;
01325 outC << "]=" << lf_bi.offset(i0,i1,i2);
01326 }
01327 outC << endl;
01328 }
01329 }
01330 else
01331 {
01332 ERRORMSG(" can not write for larger than three dimensions " << endl);
01333 }
01334 outC << "------------------ " << endl;
01335 }
01336 }
01337
01339
01340
01341
01342
01343
01344 template<class T, unsigned Dim>
01345 void BareField<T,Dim>::Compress() const
01346 {
01347 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
01348 TAU_PROFILE("BareField::Compress()", taustr, TAU_FIELD);
01349
01350 if (!compressible_m) return;
01351
01352
01353 BareField<T,Dim>& ncf = const_cast<BareField<T,Dim>&>(*this);
01354 for (iterator_if lf_i = ncf.begin_if(); lf_i != ncf.end_if(); ++lf_i)
01355 (*lf_i).second->TryCompress(isDirty());
01356 }
01357
01358 template<class T, unsigned Dim>
01359 void BareField<T,Dim>::Uncompress() const
01360 {
01361 TAU_TYPE_STRING(taustr, CT(*this) + " void ()" );
01362 TAU_PROFILE("BareField::Uncompress()", taustr, TAU_FIELD);
01363
01364
01365 BareField<T,Dim>& ncf = const_cast<BareField<T,Dim>&>(*this);
01366 for (iterator_if lf_i = ncf.begin_if(); lf_i != ncf.end_if(); ++lf_i)
01367 (*lf_i).second->Uncompress();
01368 }
01369
01370
01371
01372
01373
01374 template<class T, unsigned Dim>
01375 double BareField<T,Dim>::CompressedFraction() const
01376 {
01377 TAU_TYPE_STRING(taustr, CT(*this) + " double ()" );
01378 TAU_PROFILE("BareField::CompressedFraction()", taustr, TAU_FIELD);
01379
01380
01381
01382 double elements[2] = { 0.0, 0.0};
01383
01384 for (const_iterator_if lf_i=begin_if(); lf_i != end_if(); ++lf_i)
01385 {
01386
01387 LField<T,Dim> &lf = *(*lf_i).second;
01388
01389 double s = lf.getOwned().size();
01390
01391 elements[0] += s;
01392
01393 if ( lf.IsCompressed() )
01394
01395 elements[1] += s;
01396 }
01397
01398 double reduced_elements[2] = { 0.0, 0.0};
01399
01400 reduce(elements,elements+2,reduced_elements,OpAddAssign());
01401
01402 return reduced_elements[1]/reduced_elements[0];
01403 }
01404
01405
01407 template<class T, unsigned Dim>
01408 void BareField<T,Dim>::Repartition(UserList* userlist)
01409 {
01410 TAU_TYPE_STRING(taustr, CT(*this) + " void (UserList * )" );
01411 TAU_PROFILE("BareField::Repartition()", taustr, TAU_FIELD);
01412
01413
01414 Layout_t *newLayout = (Layout_t *)( userlist );
01415
01416
01417 BareField<T,Dim> tempField( *newLayout, Gc );
01418
01419
01420 tempField = *this;
01421
01422
01423 Locals_ac = tempField.Locals_ac;
01424
01425
01426 }
01427
01428
01430
01431
01432 template<class T, unsigned Dim>
01433 void BareField<T,Dim>::notifyUserOfDelete(UserList* userlist)
01434 {
01435 TAU_TYPE_STRING(taustr, CT(*this) + " void (UserList * )" );
01436 TAU_PROFILE("BareField::notifyUserOfDelete()", taustr, TAU_FIELD);
01437
01438
01439
01440 if (Layout != 0 && Layout->get_Id() == userlist->getUserListID()) {
01441
01442
01443
01444 Layout = 0;
01445 } else {
01446
01447
01448
01449 WARNMSG("BareField: notified of unknown UserList being deleted.");
01450 WARNMSG(endl);
01451 }
01452 }
01453
01454
01455
01456 template<unsigned Dim, bool exists>
01457 class LFieldDimTag {
01458 #ifdef IPPL_PURIFY
01459
01460 public:
01461 LFieldDimTag() {}
01462 LFieldDimTag(const LFieldDimTag<Dim,exists> &) {}
01463 LFieldDimTag<Dim,exists>&
01464 operator=(const LFieldDimTag<Dim,exists> &) { return *this; }
01465 #endif
01466 };
01467
01468
01469 template <class T>
01470 inline
01471 T* PtrOffset(T* ptr, const NDIndex<1U>& pos, const NDIndex<1U>& alloc,
01472 LFieldDimTag<1U,true>) {
01473 T* newptr = ptr + pos[0].first() - alloc[0].first();
01474 return newptr;
01475 }
01476
01477 template <class T>
01478 inline
01479 T* PtrOffset(T* ptr, const NDIndex<2U>& pos, const NDIndex<2U>& alloc,
01480 LFieldDimTag<2U,true>) {
01481 T* newptr = ptr + pos[0].first() - alloc[0].first() +
01482 alloc[0].length() * ( pos[1].first() - alloc[1].first() );
01483 return newptr;
01484 }
01485
01486 template <class T>
01487 inline
01488 T* PtrOffset(T* ptr, const NDIndex<3U>& pos, const NDIndex<3U>& alloc,
01489 LFieldDimTag<3U,true>) {
01490 T* newptr = ptr + pos[0].first() - alloc[0].first() +
01491 alloc[0].length() * ( pos[1].first() - alloc[1].first() +
01492 alloc[1].length() * ( pos[2].first() - alloc[2].first() ) );
01493 return newptr;
01494 }
01495
01496 template <class T, unsigned Dim>
01497 inline
01498 T* PtrOffset(T* ptr, const NDIndex<Dim>& pos, const NDIndex<Dim>& alloc,
01499 LFieldDimTag<Dim,false>) {
01500 T* newptr = ptr;
01501 int n=1;
01502 for (unsigned int d=0; d<Dim; ++d) {
01503 newptr += n * (pos[d].first() - alloc[d].first());
01504 n *= alloc[d].length();
01505 }
01506 return newptr;
01507 }
01508
01509
01511
01512
01513
01514
01515 template<class T, unsigned Dim>
01516 T& BareField<T,Dim>::localElement(const NDIndex<Dim>& Indexes) const
01517 {
01518 TAU_PROFILE_STMT(T tauT);
01519 TAU_TYPE_STRING(taustr, "BareField<" + CT(tauT) + ",Dim> (" + CT(Indexes) + " )" );
01520 TAU_PROFILE("localElement()", taustr, TAU_FIELD);
01521
01522
01523
01524
01525
01526
01527 const_iterator_if lf_i = begin_if();
01528 const_iterator_if lf_end = end_if();
01529 for ( ; lf_i != lf_end ; ++lf_i ) {
01530 LField<T,Dim>& lf(*(*lf_i).second);
01531
01532 if ( lf.getOwned().contains( Indexes ) ) {
01533
01534
01535 lf.Uncompress();
01536
01537
01538 NDIndex<Dim> alloc = lf.getAllocated();
01539 T* pdata = PtrOffset(lf.getP(), Indexes, alloc,
01540 LFieldDimTag<Dim,(Dim<=3)>());
01541 return *pdata;
01542 }
01543 }
01544
01545
01546 ERRORMSG("BareField::localElement: attempt to access non-local index ");
01547 ERRORMSG(Indexes << " on node " << Ippl::myNode() << endl);
01548 ERRORMSG("Occurred in a BareField with layout = " << getLayout() << endl);
01549 ERRORMSG("Calling abort ..." << endl);
01550 Ippl::abort();
01551 return *((*((*(begin_if())).second)).begin());
01552 }
01553
01554
01556
01557
01558
01559
01561 template <class T, unsigned int Dim>
01562 void BareField<T,Dim>::getsingle(const NDIndex<Dim>& Indexes, T& r) const
01563 {
01564 TAU_TYPE_STRING(taustr, "void (" + CT(Indexes) + ", " + CT(r) + " )" );
01565 TAU_PROFILE("BareField::getsingle()", taustr, TAU_FIELD);
01566
01567
01568
01569
01570
01571
01572
01573
01574 if ( (!Indexes.touches(Layout->getDomain())) &&
01575 Indexes.touches(AddGuardCells(Layout->getDomain(),Gc)) ) {
01576 getsingle_bc(Indexes,r);
01577 return;
01578 }
01579
01580
01581 int tag = Ippl::Comm->next_tag(F_GETSINGLE_TAG, F_TAG_CYCLE);
01582
01583
01584 const_iterator_if lf_end = end_if();
01585 for (const_iterator_if lf_i = begin_if(); lf_i != lf_end ; ++lf_i) {
01586 if ( (*lf_i).second->getOwned().contains( Indexes ) ) {
01587
01588
01589 LField<T,Dim>& lf( *(*lf_i).second );
01590 typename LField<T,Dim>::iterator lp = lf.begin(Indexes);
01591
01592
01593 r = *lp;
01594
01595
01596 if (Ippl::getNodes() > 1) {
01597 Message *mess = new Message;
01598 ::putMessage(*mess, r);
01599 Ippl::Comm->broadcast_others(mess, tag);
01600 }
01601
01602 return;
01603 }
01604 }
01605
01606
01607 if (Ippl::getNodes() > 1) {
01608 int any_node = COMM_ANY_NODE;
01609 Message *mess = Ippl::Comm->receive_block(any_node,tag);
01610 ::getMessage(*mess, r);
01611 delete mess;
01612 }
01613 else {
01614
01615
01616 ERRORMSG("Illegal single-value index " << Indexes);
01617 ERRORMSG(" in BareField::getsingle()" << endl);
01618 }
01619 }
01620
01622
01623
01624
01625
01626
01627
01628
01629
01630 template<class T, unsigned D>
01631 void
01632 BareField<T,D>::getsingle_bc(const NDIndex<D>& Indexes, T& r) const
01633 {
01634
01635
01636
01637
01638 int authority_proc = -1;
01639
01640
01641
01642 const_iterator_if lf_end = end_if();
01643 for (const_iterator_if lf_i = begin_if(); lf_i != lf_end ; ++lf_i) {
01644
01645 if ( (*lf_i).second->getAllocated().contains( Indexes ) ) {
01646
01647 authority_proc = Ippl::myNode();
01648
01649
01650 LField<T,D>& lf( *(*lf_i).second );
01651 typename LField<T,D>::iterator lp = lf.begin(Indexes);
01652 r = *lp;
01653
01654
01655 if (Ippl::getNodes() > 1)
01656 break;
01657
01658
01659 return;
01660 }
01661 }
01662
01663
01664
01665 typename FieldLayout<D>::touch_range_dv range =
01666 Layout->touch_range_rdv(Indexes,Gc);
01667 for (typename FieldLayout<D>::touch_iterator_dv p=range.first;
01668 p!=range.second;++p)
01669
01670 if ( authority_proc > (*p).second->getNode() )
01671
01672 authority_proc = (*p).second->getNode();
01673
01674
01675 int tag = Ippl::Comm->next_tag(F_GETSINGLE_TAG, F_TAG_CYCLE);
01676
01677 if ( authority_proc == Ippl::myNode() ) {
01678
01679 Message *mess = new Message;
01680 ::putMessage(*mess, r);
01681 Ippl::Comm->broadcast_others(mess, tag);
01682 }
01683 else {
01684
01685 Message *mess = Ippl::Comm->receive_block(authority_proc,tag);
01686 ::getMessage(*mess, r);
01687 delete mess;
01688 }
01689
01690 }
01691
01692
01693
01694
01695
01696