group_baryon_operator_w.cc

Go to the documentation of this file.
00001 // $Id: group_baryon_operator_w.cc,v 1.27 2009/03/19 17:17:20 mcneile Exp $
00002 /*! \file
00003  *  \brief Construct group baryon operators
00004  */
00005 
00006 
00007 #include "chromabase.h"
00008 
00009 #include "meas/hadron/group_baryon_operator_w.h"
00010 #include "meas/hadron/baryon_operator_factory_w.h"
00011 #include "meas/hadron/baryon_operator_aggregate_w.h"
00012 
00013 #include "meas/smear/quark_smearing_factory.h"
00014 #include "meas/smear/quark_smearing_aggregate.h"
00015 
00016 #include "meas/smear/link_smearing_aggregate.h"
00017 #include "meas/smear/link_smearing_factory.h"
00018 
00019 #include "meas/sources/source_smearing_aggregate.h"
00020 #include "meas/sources/source_smearing_factory.h"
00021 #include "meas/sinks/sink_smearing_aggregate.h"
00022 #include "meas/sinks/sink_smearing_factory.h"
00023 
00024 #include "meas/sources/dilutezN_source_const.h"
00025 #include "meas/sources/zN_src.h"
00026 #include "meas/smear/quark_source_sink.h"
00027 #include "meas/smear/displacement.h"
00028 
00029 #include "util/ferm/diractodr.h"
00030 #include "util/ft/sftmom.h"
00031 
00032 #include <map>
00033 
00034 using std::map;
00035 
00036 namespace Chroma
00037 { 
00038 
00039   void displacementSub(const multi1d<LatticeColorMatrix>& u, 
00040                        LatticeFermion& chi, 
00041                        int length, int dir)
00042   {
00043     LatticeFermion tmp;
00044     if (length > 0)
00045       for(int n = 0; n < length; ++n)
00046         {
00047           tmp = shift(chi, FORWARD, dir);
00048           chi = u[dir] * tmp;
00049         }
00050     else // If length = or < 0.  If length == 0, does nothing.
00051       for(int n = 0; n > length; --n)
00052         {
00053           tmp = shift(adj(u[dir])*chi, BACKWARD, dir);
00054           chi = tmp;
00055         }
00056   }
00057   
00058 
00059 
00060           //! Baryon sequential sources
00061           /*! \ingroup hadron */
00062           namespace GroupBaryonOperatorEnv
00063           {                     
00064                         //      =======
00065                 //      Readers
00066                         //      =======
00067 
00068                         void read( XMLReader& xml, const string& path, 
00069                                    GroupBaryonOperatorEnv::Params::Qprop_t::Solutions_t& input )
00070                 { 
00071                                 XMLReader inputtop(xml, path);
00072                   read(inputtop, "Soln_file_names", input.soln_file_names);
00073                 }
00074                 void read( XMLReader& xml, const string& path, 
00075                                    GroupBaryonOperatorEnv::Params::Qprop_t& input )
00076                 { 
00077                                 XMLReader inputtop(xml, path);
00078                   read(inputtop, "Quarks", input.solns);
00079                 }
00080                 void read( XMLReader& xml, const string& path, 
00081                                    GroupBaryonOperatorEnv::Params::dilution_t& input )
00082                 { 
00083                                 XMLReader inputtop(xml, path);
00084                                 input.spatial_mask_size.resize(3);
00085                   read(inputtop, "spatial_mask_size", input.spatial_mask_size);
00086                                 input.spatial_mask.resize(3);
00087                   read(inputtop, "spatial_mask", input.spatial_mask);
00088                                 input.spin_mask.resize(4);
00089                   read(inputtop, "spin_mask", input.spin_mask);
00090                                 input.color_mask.resize(3);
00091                   read(inputtop, "color_mask", input.color_mask);
00092                 }
00093                 void read( XMLReader& xml, const string& path, 
00094                                    GroupBaryonOperatorEnv::Params& input )
00095                 { 
00096                                 XMLReader inputtop(xml, path);
00097                   read(inputtop, "Qprops", input.qprop);
00098                   read(inputtop, "Cfg", input.gaugestuff.cfg);
00099                   read(inputtop, "DilutionScheme", input.dilution);
00100                 }
00101                 // Writer
00102                 void write( XMLWriter& xml, const string& path, 
00103                                     const GroupBaryonOperatorEnv::Params& param )
00104                 {
00105                   param.writeXML( xml, path );
00106                 }
00107 
00108                 //      ============
00109                 //      Params stuff
00110                 //      ============
00111                 
00112     //! Initialize
00113     Params::Params()
00114     {}
00115 
00116     //! Read parameters
00117     Params::Params( XMLReader& xml, const std::string& path )
00118     {
00119       XMLReader paramtop( xml, path );
00120                         // ============================
00121                         // number of quarks fixed to 3
00122                         // j_decay dir is 3 (time)
00123                         // space-time is 4-dimensions
00124                         // ============================
00125                         int NumberofQuarks = 3;
00126                         j_decay = 3;
00127                         nrow.resize( 4 );
00128       int version;
00129       read( paramtop, "Param/Name", name );
00130       read( paramtop, "Param/version", version );
00131                         QDPIO::cout<< name << " Version " << version <<endl;
00132       switch ( version )
00133       {
00134         case 2:
00135           break;
00136         default:
00137           QDPIO::cerr << name << ": parameter version " << version << " unsupported." << endl;
00138           QDP_abort( 1 );
00139       }
00140                         read( paramtop, "Param/LattSize", nrow );
00141       read( paramtop, "Param/Baryon_type", param.baryon_operator );
00142 
00143       read( paramtop, "Param/SourceSmearing", source_smearing );
00144       read( paramtop, "Param/SinkSmearing", sink_smearing );
00145 
00146       gaugestuff.link_smearing = readXMLGroup( paramtop, "Param/LinkSmearing", "LinkSmearingType" );
00147                         read( paramtop, "Param/InputFileName", InputFileName );
00148                         QDPIO::cout<< "Main input file is " << InputFileName <<endl;
00149                         read( paramtop, "Cfg", gaugestuff.cfg );
00150 
00151                         dilution.resize( NumberofQuarks );
00152                         for(int i=0; i < NumberofQuarks; ++i) dilution[ i ].N = 4; // Z(4) noise
00153                         for(int i=0; i < NumberofQuarks; ++i) dilution[ i ].j_decay = 3; // time-direction
00154 
00155                 } // end Params::Params
00156 
00157 
00158     // Writer
00159     void Params::writeXML( XMLWriter& xml, const string& path ) const
00160     {
00161       push( xml, path );
00162       int version = 1;
00163       write( xml, "version", version );
00164       write( xml, "nrow", nrow );
00165       //write( xml, "BaryonOperatorType", GroupBaryonOperatorEnv::name );
00166                         //xml << source_smearing.xml;
00167       //xml << sink_smearing.xml;
00168       //xml << link_smearing.xml;
00169       write( xml, "InputFileName", InputFileName );
00170       write( xml, "mom2_max", mom2_max );
00171       write( xml, "j_decay",  j_decay );
00172       pop( xml );
00173     } // end void Params::writeXML
00174 
00175 
00176                 // Used for hybrid list indices, NH and noise indices 
00177                 //      for indices:  ( ord, i,   j,   k,   which )                     
00178                 //      for NHlimit:  ( ord, NH0, NH1, NH2, which )                     
00179                 //      for NoiseInd: ( ord, 0,   1,   2,   which )                     
00180                 int DilSwap( int ord, int i, int j, int k, int which )   
00181                 {  
00182                   multi1d<int> ijk( 3 );                                                                                                                           
00183                   switch ( ord )                                                                                                                                                           
00184                   {                                                                                                                                                                                                        
00185                     case 0:  // [012] ijk                                                                                                                  
00186                       ijk[ 0 ] = i;                                                                                                                                        
00187                       ijk[ 1 ] = j;                                                                                                                                        
00188                       ijk[ 2 ] = k;                                                                                                                                        
00189                       break;                                                                                                                                                                       
00190                     case 1:  // [210] kji                                                                                                                  
00191                       ijk[ 0 ] = k;                                                                                                                                        
00192                       ijk[ 1 ] = j;                                                                                                                                        
00193                       ijk[ 2 ] = i;                                                                                                                                        
00194                       break;                                                                                                                                                                       
00195                     case 2:  // [021] ikj                                                                                                                  
00196                       ijk[ 0 ] = i;                                                                                                                                        
00197                       ijk[ 1 ] = k;                                                                                                                                        
00198                       ijk[ 2 ] = j;                                                                                                                                        
00199                       break;                                                                                                                                                                       
00200                     case 3:  // [102] jik                                                                                                                  
00201                       ijk[ 0 ] = j;                                                                                                                                        
00202                       ijk[ 1 ] = i;                                                                                                                                        
00203                       ijk[ 2 ] = k;                                                                                                                                        
00204                       break;                                                                                                                                                                       
00205                     case 4:  // [120] jki                                                                                                                  
00206                       ijk[ 0 ] = j;                                                                                                                                        
00207                       ijk[ 1 ] = k;                                                                                                                                        
00208                       ijk[ 2 ] = i;
00209                       break;                                                                                                                                                                       
00210                     case 5:  // [201] kij                                                                                                                  
00211                       ijk[ 0 ] = k;                                                                                                                                        
00212                       ijk[ 1 ] = i;                                                                                                                                        
00213                       ijk[ 2 ] = j;                                                                                                                                        
00214                       break;
00215                                 default:
00216                                         cerr<<"say what? "<<ord<<" "<<i<<" "<<j<<" "<<k<<" "<<which<<endl;
00217                                         exit(1);
00218                   }                                                                                                                                                                                                        
00219                   return ( ijk[ which ] );                                                                                                                 
00220                 }
00221                 //
00222                 // index things back to [012] ordering
00223                 //
00224                 int DilSwapInv( int ord, int i, int j, int k, int which )   
00225                 {  
00226                   multi1d<int> ijk( 3 );                                                                                                                           
00227                   switch ( ord )                                                                                                                                                           
00228                   {                                                                                                                                                                                                        
00229                     case 0:  // [012] ijk                                                                                                                  
00230                       ijk[ 0 ] = i;                                                                                                                                        
00231                       ijk[ 1 ] = j;                                                                                                                                        
00232                       ijk[ 2 ] = k;                                                                                                                                        
00233                       break;                                                                                                                                                                       
00234                     case 1:  // [210] kji                                                                                                                  
00235                       ijk[ 0 ] = k;                                                                                                                                        
00236                       ijk[ 1 ] = j;                                                                                                                                        
00237                       ijk[ 2 ] = i;                                                                                                                                        
00238                       break;                                                                                                                                                                       
00239                     case 2:  // [021] ikj                                                                                                                  
00240                       ijk[ 0 ] = i;                                                                                                                                        
00241                       ijk[ 1 ] = k;                                                                                                                                        
00242                       ijk[ 2 ] = j;                                                                                                                                        
00243                       break;                                                                                                                                                                       
00244                     case 3:  // [102] jik                                                                                                                  
00245                       ijk[ 0 ] = j;                                                                                                                                        
00246                       ijk[ 1 ] = i;                                                                                                                                        
00247                       ijk[ 2 ] = k;                                                                                                                                        
00248                       break;                                                                                                                                                                       
00249                     case 4:  // [120] jki                                                                                                                  
00250                       ijk[ 0 ] = k;                                                                                                                                        
00251                       ijk[ 1 ] = i;                                                                                                                                        
00252                       ijk[ 2 ] = j; 
00253                       break;                                                                                                                                                                       
00254                     case 5:  // [201] kij                                                                                                                  
00255                       ijk[ 0 ] = j;                                                                                                                                        
00256                       ijk[ 1 ] = k;                                                                                                                                        
00257                       ijk[ 2 ] = i;                                                                                                                                        
00258                       break;
00259                                 default:
00260                                         cerr<<"Say what? "<<ord<<" "<<i<<" "<<j<<" "<<k<<" "<<which<<endl;
00261                                         exit(1);
00262                   }
00263                   return ( ijk[ which ] );                                                                                                                 
00264                 }
00265                 
00266                 //      ====================
00267                 //      GroupBaryonQQQ stuff
00268                 //      ====================
00269                 
00270     //! Constructor
00271     GroupBaryonQQQ::GroupBaryonQQQ()
00272     {
00273       // The spin basis matrix to goto Dirac
00274       //spin_rotate_mat = adj( DiracToDRMat() );
00275                 }
00276     void GroupBaryonQQQ::init( const Params& p_ )
00277     {
00278                         myparams = p_;
00279                         quark.resize(3);
00280       // The spin basis matrix to goto Dirac
00281       spin_rotate_mat = adj( DiracToDRMat() );
00282     } // GroupBaryonQQQ::init
00283 
00284     //! Full constructor
00285     GroupBaryonQQQ::GroupBaryonQQQ( const Params& p ) :
00286         myparams( p )
00287                 {
00288       // The spin basis matrix to goto Dirac
00289       spin_rotate_mat = adj( DiracToDRMat() );
00290       // Factory constructions
00291       try
00292       {
00293         // Create the source quark smearing object
00294         {
00295           std::istringstream xml_s( myparams.source_smearing.source.xml );
00296           XMLReader smeartop( xml_s );
00297           const string smear_path = "/SourceSmearing";
00298           QDPIO::cout << "Source quark smearing type = " 
00299                                                     << myparams.source_smearing.source.id << endl;
00300                       Handle< QuarkSourceSink<LatticeFermion> > sourceSmearing(
00301                                                         TheFermSourceSmearingFactory::Instance().createObject(
00302                                                myparams.source_smearing.source.id, 
00303                                                                                          smeartop,
00304                                                                                          myparams.source_smearing.source.path, 
00305                                                                                          myparams.gaugestuff.u )
00306                                                         );
00307         }
00308         // Create the sink quark smearing object
00309         {
00310           std::istringstream xml_s( myparams.sink_smearing.sink.xml );
00311           XMLReader smeartop( xml_s );
00312           const string smear_path = "/SinkSmearing";
00313           QDPIO::cout << "Sink quark smearing type = " 
00314                                                     << myparams.sink_smearing.sink.id << endl;
00315                                         Handle< QuarkSourceSink<LatticeFermion> > sinkSmearing(
00316                                                         TheFermSinkSmearingFactory::Instance().createObject(
00317                                                                                         myparams.sink_smearing.sink.id, 
00318                                                                                         smeartop,
00319                                                                                         myparams.sink_smearing.sink.path, 
00320                                                                                         myparams.gaugestuff.u )
00321                                                         );
00322         }
00323       }
00324       catch ( const std::string & e )
00325       {
00326         QDPIO::cerr << name << ": Caught Exception smearing: " << e << endl;
00327         QDP_abort( 1 );
00328       }
00329 
00330     } // GroupBaryonQQQ::GroupBaryonQQQ
00331 
00332     //! Reader
00333     void ReadTextInput( Params& params,
00334                         multi1d<GroupBaryonOp>& AB, multi1d<GroupBaryonOp>& CB,
00335                         multi1d<GroupBaryonQQQ>& AQQQ, multi1d<GroupBaryonQQQ>& CQQQ )
00336     {
00337       QDPIO::cout << "Reading input from TEXT file : " << params.InputFileName <<endl;
00338       TextFileReader reader( params.InputFileName );
00339 // for testing
00340 #ifdef REDUCETOTIMEDILUTION
00341                         params.NdilReduce = 12;
00342 #else
00343                         params.NdilReduce = 1;
00344 #endif
00345                         // Lattice sizes
00346       // ------------------------------------------------------------------------------------
00347       reader >> params.nrow[ 0 ] >> params.nrow[ 1 ] >> params.nrow[ 2 ] >> params.nrow[ 3 ];
00348                         // ------------------------------------------------------------------------------------
00349 
00350                         params.Nmomenta = 1; // need to change this in the future when we do multiple momenta
00351                         params.mom2_max = 0; // this obviously needs to be changed and read in ...
00352                         params.j_decay  = 3; // this is OK
00353 
00354       // Orderings are more or less hard wired for now
00355                         //  ... to be continued later when I have the time
00356                         
00357                         //read( xml, "NumSourcePerm", params.NsrcOrderings );
00358                         // ------------------------------------------------------------------------------------
00359       reader >> params.NsrcOrderings;
00360                         // ------------------------------------------------------------------------------------
00361                         // 012    factor in front is  2  2
00362       // 210    factor in front is -2  2
00363       // 021    factor in front is  1 -1
00364       // 102    factor in front is  1 -1
00365       // 120    factor in front is -1 -1
00366       // 201    factor in front is -1 -1
00367 
00368       //read( xml, "SrcQuarkIndices", params.SrcOrderings );
00369       params.SrcOrderings.resize( params.NsrcOrderings );
00370       for(int i=0; i < params.NsrcOrderings; ++i)
00371       {
00372         params.SrcOrderings[ i ].resize( 3 );
00373                                 // ------------------------------------------------------------------------------------
00374         reader >> params.SrcOrderings[ i ][ 0 ] >> params.SrcOrderings[ i ][ 1 ] >> params.SrcOrderings[ i ][ 2 ];
00375                                 // ------------------------------------------------------------------------------------
00376       }
00377       //
00378                         // This again, is hard-wired to 012 case for now ...
00379                         //
00380                         //read( xml, "NumSinkPerm", params.NsnkOrderings );
00381                         // ------------------------------------------------------------------------------------
00382       reader >> params.NsnkOrderings;
00383                         // ------------------------------------------------------------------------------------
00384       params.SnkOrderings.resize( params.NsnkOrderings );
00385       //read( xml, "SnkQuarkIndices", params.SnkOrderings );
00386       for(int i=0; i < params.NsnkOrderings; ++i)
00387       {
00388         params.SnkOrderings[ i ].resize( 3 );
00389                                 // ------------------------------------------------------------------------------------
00390         reader >> params.SnkOrderings[ i ][ 0 ] >> params.SnkOrderings[ i ][ 1 ] >> params.SnkOrderings[ i ][ 2 ];
00391                                 // ------------------------------------------------------------------------------------
00392       }
00393       // Hybrid list sizes
00394       params.NH.resize( params.NsrcOrderings ); 
00395       for(int i=0; i < params.NH.size(); ++i) params.NH[ i ].resize( 3 ); // L,M,R quarks
00396         // Set the hybrid list sizes
00397         //   { something like this would be better                  }
00398         //   { Ndil = params.named_obj.prop.op[n].soln_files.size() }                   
00399       //XMLReader& xml;
00400       //read( xml, "Ndilutions", params.NH[0] );
00401                         // This is the default 012 ordering
00402                         // ------------------------------------------------------------------------------------
00403       reader >> params.NH[ 0 ][ 0 ] >> params.NH[ 0 ][ 1 ] >> params.NH[ 0 ][ 2 ];
00404                         // ------------------------------------------------------------------------------------
00405       for(int i=1; i < params.NsrcOrderings; ++i)
00406                         {
00407         for(int j=0; j < 3; ++j)
00408                                         params.NH[ i ][ j ] = params.NH[ 0 ][ (params.SrcOrderings[ i ][ j ]) ];
00409                         }
00410       //read( xml, "NumOperators", params.Noperators );
00411       // ------------------------------------------------------------------------------------
00412       reader >> params.Noperators;
00413                         // ------------------------------------------------------------------------------------
00414       // The Baryon operator names ; G1g_L3_TDT_25 ...
00415       //read( xml, "OperatorNames", params.Names );
00416       params.Names.resize( params.Noperators );
00417                         int nameindex;
00418       for(int i=0; i < params.Names.size(); ++i)
00419       {
00420                                 // ------------------------------------------------------------------------------------
00421         reader >> nameindex >> params.Names[ i ];
00422                                 // ------------------------------------------------------------------------------------
00423       }
00424       //read( xml, "NumDistinctQQQ", params.NQQQs );
00425                         // ------------------------------------------------------------------------------------
00426       reader >> params.NQQQs;
00427                         // ------------------------------------------------------------------------------------
00428 #ifdef MAKE_SINK_OPERATORS
00429       AQQQ.resize( params.NQQQs );
00430       for(int i=0; i < params.NQQQs; ++i)
00431       {
00432                                 AQQQ[ i ].init( params );
00433                         }
00434 #endif
00435 #ifdef MAKE_SOURCE_OPERATORS
00436       CQQQ.resize( params.NQQQs );
00437       for(int i=0; i < params.NQQQs; ++i)
00438       {
00439                                 CQQQ[ i ].init( params );
00440                         }
00441 #endif
00442                         //
00443                         // GroupBaryonOperator structures
00444                         //
00445 #ifdef MAKE_SINK_OPERATORS
00446       AB.resize( params.Noperators );
00447 #endif
00448 #ifdef MAKE_SOURCE_OPERATORS
00449       CB.resize( params.Noperators );
00450 #endif
00451                         //
00452                         // Using the BaryonOperator_t structure
00453                         //
00454       for(int n=0; n < params.Noperators; ++n)
00455       {
00456 #ifdef MAKE_SINK_OPERATORS
00457         // This is 1 now because all the 
00458                                 // contractions are summed over
00459                                 // AB[ n ].baryonoperator.orderings.resize( params.NsnkOrderings );
00460         // for(int ord=0; ord < params.NsnkOrderings; ++ord)
00461                                 AB[ n ].baryonoperator.orderings.resize( 1 );
00462         for(int ord=0; ord < 1; ++ord)
00463                                 {
00464                 AB[ n ].baryonoperator.orderings[ ord ].op.resize( params.NH[ 0 ][ 0 ]/params.NdilReduce, params.NH[ 0 ][ 1 ]/params.NdilReduce, params.NH[ 0 ][ 2 ]/params.NdilReduce );
00465                 for(int i=0; i < (params.NH[ 0 ][ 0 ]/params.NdilReduce); ++i)
00466                   for(int j=0; j < (params.NH[ 0 ][ 1 ]/params.NdilReduce); ++j)
00467                     for(int k=0; k < (params.NH[ 0 ][ 2 ]/params.NdilReduce); ++k)
00468                     {
00469                       AB[ n ].baryonoperator.orderings[ ord ].op( i, j, k ).ind.resize(1);
00470                       AB[ n ].baryonoperator.orderings[ ord ].op( i, j, k ).ind[ 0 ].elem.resize( params.Nmomenta, params.nrow[3] );
00471                       for(int t=0; t<params.nrow[3]; ++t)
00472                                                                         AB[ n ].baryonoperator.orderings[ ord ].op( i, j, k ).ind[ 0 ].elem( 0, t ) = cmplx( Real(0), Real(0) );
00473                                                         }
00474         }
00475                                 AB[ n ].baryonoperator.version = 1;
00476                                 AB[ n ].baryonoperator.mom2_max = params.mom2_max;
00477                                 AB[ n ].baryonoperator.j_decay = params.j_decay;
00478                                 //AB[ n ].baryonoperator.quarkOrderings = params.SnkOrderings;
00479                                 AB[ n ].baryonoperator.quarkOrderings.resize(1);
00480                                 AB[ n ].baryonoperator.quarkOrderings[0].resize(3);
00481                                 AB[ n ].baryonoperator.quarkOrderings[0][0] = 0;
00482                                 AB[ n ].baryonoperator.quarkOrderings[0][1] = 1;
00483                                 AB[ n ].baryonoperator.quarkOrderings[0][2] = 2;
00484 #endif
00485 #ifdef MAKE_SOURCE_OPERATORS
00486                                 //CB[ n ].baryonoperator.orderings.resize( params.NsrcOrderings );
00487         //for(int ord=0; ord < params.NsrcOrderings; ++ord)
00488                                 CB[ n ].baryonoperator.orderings.resize( 1 );
00489         for(int ord=0; ord < 1; ++ord)
00490                                 {
00491         CB[ n ].baryonoperator.orderings[ ord ].op.resize( params.NH[ 0 ][ 0 ]/params.NdilReduce, params.NH[ 0 ][ 1 ]/params.NdilReduce, params.NH[ 0 ][ 2 ]/params.NdilReduce );
00492         for(int i=0; i < (params.NH[ 0 ][ 0 ]/params.NdilReduce); ++i)
00493           for(int j=0; j < (params.NH[ 0 ][ 1 ]/params.NdilReduce); ++j)
00494             for(int k=0; k < (params.NH[ 0 ][ 2 ]/params.NdilReduce); ++k)
00495             {
00496               CB[ n ].baryonoperator.orderings[ ord ].op( i, j, k ).ind.resize(1);
00497               CB[ n ].baryonoperator.orderings[ ord ].op( i, j, k ).ind[ 0 ].elem.resize( params.Nmomenta, params.nrow[3] );
00498               for(int t=0; t<params.nrow[3]; ++t) CB[ n ].baryonoperator.orderings[ ord ].op( i, j, k ).ind[ 0 ].elem( 0, t ) = cmplx( Real(0), Real(0) );
00499             }
00500                                 }
00501                                 CB[ n ].baryonoperator.version = 1;
00502                                 CB[ n ].baryonoperator.mom2_max = params.mom2_max;
00503                                 CB[ n ].baryonoperator.j_decay = params.j_decay;
00504                                 //CB[ n ].baryonoperator.quarkOrderings = params.SrcOrderings;
00505                                 CB[ n ].baryonoperator.quarkOrderings.resize(1);
00506                                 CB[ n ].baryonoperator.quarkOrderings[0].resize(3);
00507                                 CB[ n ].baryonoperator.quarkOrderings[0][0] = 0;
00508                                 CB[ n ].baryonoperator.quarkOrderings[0][1] = 1;
00509                                 CB[ n ].baryonoperator.quarkOrderings[0][2] = 2;
00510 #endif
00511       } // end n (operators)
00512       for(int i=0; i < params.NQQQs; ++i)
00513       {
00514                                 int dL, hash;
00515                                 int q0spin,q1spin,q2spin,q0disp,q1disp,q2disp,nbops;
00516                                 reader >> hash >> q0spin >> q1spin >> q2spin >> q0disp >> q1disp >> q2disp >> dL >> nbops;
00517 #ifdef MAKE_SINK_OPERATORS
00518                                 AQQQ[ i ].quark[ 0 ].spin = q0spin -1;
00519                                 AQQQ[ i ].quark[ 1 ].spin = q1spin -1;
00520                                 AQQQ[ i ].quark[ 2 ].spin = q2spin -1;
00521                                 AQQQ[ i ].quark[ 0 ].displacement = q0disp;
00522                                 AQQQ[ i ].quark[ 1 ].displacement = q1disp;
00523                                 AQQQ[ i ].quark[ 2 ].displacement = q2disp;
00524               AQQQ[ i ].NBaryonOps = nbops;
00525         for(int j=0; j < 3; ++j)
00526         {
00527           if ( AQQQ[ i ].quark[ j ].displacement >= 0 )
00528                                         {
00529                                                 AQQQ[ i ].quark[ j ].disp_len = dL;
00530                                         } 
00531                                         else 
00532                                         {
00533                                                 AQQQ[ i ].quark[ j ].disp_len = -dL;
00534                                         }
00535                                 }
00536 #endif
00537 #ifdef MAKE_SOURCE_OPERATORS
00538                                 CQQQ[ i ].quark[ 0 ].spin = q0spin -1;
00539                                 CQQQ[ i ].quark[ 1 ].spin = q1spin -1;
00540                                 CQQQ[ i ].quark[ 2 ].spin = q2spin -1;
00541                                 CQQQ[ i ].quark[ 0 ].displacement = q0disp;
00542                                 CQQQ[ i ].quark[ 1 ].displacement = q1disp;
00543                                 CQQQ[ i ].quark[ 2 ].displacement = q2disp;
00544               CQQQ[ i ].NBaryonOps = nbops;
00545         for(int j=0; j < 3; ++j)
00546         {
00547           if ( CQQQ[ i ].quark[ j ].displacement >= 0 )
00548                                         {
00549                                                 CQQQ[ i ].quark[ j ].disp_len = dL;
00550                                         } 
00551                                         else 
00552                                         {
00553                                                 CQQQ[ i ].quark[ j ].disp_len = -dL;
00554                                         }
00555                                 }
00556 #endif
00557 #ifdef MAKE_SINK_OPERATORS
00558         for(int j=0; j < 3; ++j)
00559         {
00560           if ( AQQQ[ i ].quark[ j ].displacement == 0 )
00561           {
00562             AQQQ[ i ].quark[ j ].disp_dir = 0;
00563             AQQQ[ i ].quark[ j ].disp_len = 0;
00564             AQQQ[ i ].quark[ j ].disp_ind = 0;
00565           }
00566           else if ( AQQQ[ i ].quark[ j ].displacement > 0 )
00567           {
00568             AQQQ[ i ].quark[ j ].disp_dir = AQQQ[ i ].quark[ j ].displacement - 1;
00569             AQQQ[ i ].quark[ j ].disp_ind = AQQQ[ i ].quark[ j ].displacement;
00570             AQQQ[ i ].quark[ j ].disp_len = dL;
00571           }
00572           else
00573           {
00574             AQQQ[ i ].quark[ j ].disp_dir = -AQQQ[ i ].quark[ j ].displacement - 1;
00575             AQQQ[ i ].quark[ j ].disp_ind = -AQQQ[ i ].quark[ j ].displacement + 3;
00576             AQQQ[ i ].quark[ j ].disp_len = -dL;
00577           }
00578         } // end j
00579         AQQQ[ i ].coef.resize( AQQQ[ i ].NBaryonOps );
00580         AQQQ[ i ].whichBaryonOps.resize( AQQQ[ i ].NBaryonOps );
00581         AQQQ[ i ].baryon.resize( AQQQ[ i ].NBaryonOps );
00582         //             name                                    coefficient
00583         // AQQQ[ <mu,nu,tau> ].whichBaryonOps[ n ]    ,   AQQQ[ <mu,nu,tau> ].coef[ n ]
00584         //
00585         for(int n=0; n < AQQQ[ i ].NBaryonOps; ++n)
00586         {
00587           Real re, im;
00588                                         // ------------------------------------------------------------------------------------
00589           reader >> AQQQ[ i ].whichBaryonOps[ n ];
00590           reader >> re >> im;
00591                                         // ------------------------------------------------------------------------------------
00592           AQQQ[ i ].coef[ n ] = cmplx( re, im );
00593                                         //
00594                                         // Storage will be in the GroupBaryonOp class
00595                                         //
00596                                         AQQQ[ i ].baryon[ n ] = &AB[ (AQQQ[ i ].whichBaryonOps[ n ]) ];
00597         }
00598 #endif
00599 #ifdef MAKE_SOURCE_OPERATORS
00600         for(int j=0; j < 3; ++j)
00601         {
00602           if ( CQQQ[ i ].quark[ j ].displacement == 0 )
00603           {
00604             CQQQ[ i ].quark[ j ].disp_dir = 0;
00605             CQQQ[ i ].quark[ j ].disp_len = 0;
00606             CQQQ[ i ].quark[ j ].disp_ind = 0;
00607           }
00608           else if ( CQQQ[ i ].quark[ j ].displacement > 0 )
00609           {
00610             CQQQ[ i ].quark[ j ].disp_dir = CQQQ[ i ].quark[ j ].displacement - 1;
00611             CQQQ[ i ].quark[ j ].disp_ind = CQQQ[ i ].quark[ j ].displacement;
00612             CQQQ[ i ].quark[ j ].disp_len = dL;
00613           }
00614           else
00615           {
00616             CQQQ[ i ].quark[ j ].disp_dir = -CQQQ[ i ].quark[ j ].displacement - 1;
00617             CQQQ[ i ].quark[ j ].disp_ind = -CQQQ[ i ].quark[ j ].displacement + 3;
00618             CQQQ[ i ].quark[ j ].disp_len = -dL;
00619           }
00620         } // end j
00621         CQQQ[ i ].coef.resize( CQQQ[ i ].NBaryonOps );
00622         CQQQ[ i ].whichBaryonOps.resize( CQQQ[ i ].NBaryonOps );
00623         CQQQ[ i ].baryon.resize( CQQQ[ i ].NBaryonOps );
00624         //             name                                    coefficient
00625         // CQQQ[ <mu,nu,tau> ].whichBaryonOps[ n ]    ,   CQQQ[ <mu,nu,tau> ].coef[ n ]
00626         //
00627         for(int n=0; n < CQQQ[ i ].NBaryonOps; ++n)
00628         {
00629           Real re, im;
00630                                         // ------------------------------------------------------------------------------------
00631           reader >> CQQQ[ i ].whichBaryonOps[ n ];
00632           reader >> re >> im;
00633                                         // ------------------------------------------------------------------------------------
00634           CQQQ[ i ].coef[ n ] = cmplx( re, im );
00635                                         //
00636                                         // Storage will be in the GroupBaryonOp class
00637                                         //
00638                                         CQQQ[ i ].baryon[ n ] = &CB[ (CQQQ[ i ].whichBaryonOps[ n ]) ];
00639         }
00640 #endif
00641       } // i loop
00642       //
00643       // Now for the creation operators
00644       //
00645       // rotate by gamma4's to get the barred coefficients
00646       //       index     0  1   2   3
00647       //      gamma4 = ( 1, 1, -1, -1 )
00648       // and take the complex conjugate
00649 #ifdef MAKE_SOURCE_OPERATORS
00650       for(int i=0; i < params.NQQQs; ++i)
00651       {
00652         int spinCount, flipsign;
00653         spinCount = 0;
00654         flipsign = 0;
00655         for(int j=0; j < 3; ++j)
00656         {
00657           if ( CQQQ[ i ].quark[ j ].spin > 1 ) spinCount++;
00658         }
00659         if ( spinCount % 2 ) flipsign = 1;
00660         //
00661                                 //             name                           coefficient
00662         // CQQQ[ <mu,nu,tau> ].whichBaryonOps[ n ] , CQQQ[ i ].coef[ n ]
00663         //
00664         for(int n=0; n < CQQQ[ i ].NBaryonOps; ++n)
00665         {
00666           if ( flipsign )
00667           {
00668             // cc and -sign
00669             // CQQQ[ i ].coef[ n ] = cmplx( -re, im );
00670             // flip sign (cc at the correlator level)
00671                                                 CQQQ[ i ].coef[ n ] *= cmplx( Real(-1), Real(0) );
00672           }
00673                                         //
00674                                         // Storage will be in the GroupBaryonOp class
00675                                         //
00676                                         CQQQ[ i ].baryon[ n ] = &CB[ CQQQ[ i ].whichBaryonOps[ n ] ];
00677         } // n : CQQQ[ i ].baryon[ n ] loop
00678       } // i : CQQQ[] loop
00679 #endif                  
00680                         // solution file names
00681                         params.qprop.solns.resize( 3 );
00682       for(int n=0; n < 3; ++n)
00683                         {
00684                                 params.qprop.solns[ n ].soln_file_names.resize( params.NH[ 0 ][ n ] );
00685         for(int i=0; i < params.NH[ 0 ][ n ]; ++i)
00686                                 {
00687                                         // ------------------------------------------------------------------------------------
00688                 reader >> params.qprop.solns[ n ].soln_file_names[ i ];
00689                                         // ------------------------------------------------------------------------------------
00690                                 }
00691                         }
00692       reader.close();
00693 
00694       QDPIO::cout << "Reading input from text file DONE " << endl;
00695     } // void ReadTextInput
00696                 
00697                 
00698     //! Construct the source
00699     LatticeFermion
00700     GroupBaryonQQQ::operator() ( Params::dilution_t diln ) const
00701     {
00702       //QDPIO::cout << "Diluted random complex ZN source" << endl;
00703       if ( diln.spatial_mask_size.size() != Nd - 1 )
00704       {
00705         QDPIO::cerr << name << ": spatial mask size incorrect 1" << endl;
00706         QDP_abort( 1 );
00707       }
00708       if ( diln.spatial_mask.size() == 0 )
00709       {
00710         QDPIO::cerr << name << ": spatial mask incorrect 2" << endl;
00711         QDP_abort( 1 );
00712       }
00713       multi1d<int> lookup_dir( Nd - 1 );
00714       int mu = 0;
00715       for(int j=0; j < diln.spatial_mask_size.size(); ++j, ++mu)
00716       {
00717         if ( j == diln.j_decay ) ++mu;  // bump up to next dir
00718         lookup_dir[ j ] = mu;
00719       }
00720       for(int j=0; j < diln.spatial_mask.size(); ++j)
00721       {
00722         if ( diln.spatial_mask[ j ].size() != Nd - 1 )
00723         {
00724           QDPIO::cerr << name << ": spatial mask incorrect 3" << endl;
00725           QDP_abort( 1 );
00726         }
00727       }
00728       for(int c=0; c < diln.color_mask.size(); ++c)
00729       {
00730         if ( diln.color_mask[ c ] < 0 || diln.color_mask[ c ] >= Nc )
00731         {
00732           QDPIO::cerr << name << ": color mask incorrect 6" << endl;
00733           QDP_abort( 1 );
00734         }
00735       }
00736       for(int s=0; s < diln.spin_mask.size(); ++s)
00737       {
00738         if ( diln.spin_mask[ s ] < 0 || diln.spin_mask[ s ] >= Ns )
00739         {
00740           QDPIO::cerr << name << ": spin mask incorrect 7" << endl;
00741           QDP_abort( 1 );
00742         }
00743       }
00744       // Save current seed
00745       Seed ran_seed;
00746       QDP::RNG::savern( ran_seed );
00747       // Set the seed to desired value
00748       QDP::RNG::setrn( diln.ran_seed );
00749       // Create the noisy quark source on the entire lattice
00750       LatticeFermion quark_noise;
00751       zN_src( quark_noise, diln.N );
00752       // This is the filtered noise source to return
00753       LatticeFermion quark_source = zero;
00754       // Filter over the color and spin indices first
00755       for(int s=0; s < diln.spin_mask.size(); ++s)
00756       {
00757         int spin_source = diln.spin_mask[ s ];
00758         LatticeColorVector colvec = peekSpin( quark_noise, spin_source );
00759         LatticeColorVector dest = zero;
00760         for(int c=0; c < diln.color_mask.size(); ++c)
00761         {
00762           int color_source = diln.color_mask[ c ];
00763           LatticeComplex comp = peekColor( colvec, color_source );
00764           pokeColor( dest, comp, color_source );
00765         }
00766         pokeSpin( quark_source, dest, spin_source );
00767       }
00768       quark_noise = quark_source;  // reset
00769       // Filter over the spatial sites
00770       LatticeBoolean mask = false;  // this is the starting mask
00771       for(int n=0; n < diln.spatial_mask.size(); ++n)
00772       {
00773         LatticeBoolean btmp = true;
00774         for(int j=0; j < diln.spatial_mask[ n ].size(); ++j)
00775           btmp &= ( Layout::latticeCoordinate( lookup_dir[ j ] ) % diln.spatial_mask_size[ j ] ) == diln.spatial_mask[ n ][ j ];
00776         mask |= btmp;
00777       }
00778       // Filter over the time slices
00779       mask &= Layout::latticeCoordinate( diln.j_decay ) == diln.t_source;
00780       // Zap the unused sites
00781       quark_source = where( mask, quark_noise, Fermion( zero ) );
00782       // Reset the seed
00783       QDP::RNG::setrn( ran_seed );
00784       return quark_source;
00785     };
00786                 
00787                 
00788     //! Construct array of maps of displacements
00789     void
00790     GroupBaryonQQQ::displaceQuarks( multi1d< map<int, LatticeFermion> >& disp_quarks,
00791                                     const multi1d<LatticeFermion>& q,
00792                                                                                                                                                 int* qindices                                   
00793                                   ) const
00794     {
00795       START_CODE();
00796       //QDPIO::cout << __PRETTY_FUNCTION__ << ": entering" << endl;
00797       disp_quarks.resize( 3 );
00798       for(int i=0; i < disp_quarks.size(); ++i)
00799       {
00800         // Make some shorthands to ease my brain
00801         map<int, LatticeFermion>& disp_q = disp_quarks[ i ];
00802                                 const QuarkTerm_t& term_q = quark[ qindices[ i ] ];
00803         // If no entry, then create a displaced version of the quark
00804         if ( disp_q.find( term_q.displacement ) == disp_q.end() )
00805         {
00806           //QDPIO::cout << __func__
00807           //               << ": i=" << i
00808           //               << " disp=" << quark[i].displacement << " disp=" << term_q.displacement
00809           //               << " len =" << quark[i].disp_len                     << " len=" << term_q.disp_len
00810           //               << " dir =" << quark[i].disp_dir                     << " dir=" << term_q.disp_dir
00811           //               << " spin=" << quark[i].spin                 << " spin=" << term_q.spin
00812           //               << endl;
00813           LatticeFermion qq = q[ i ];
00814           displacementSub( myparams.gaugestuff.u, qq, term_q.disp_len, term_q.disp_dir );
00815           disp_q.insert( std::make_pair( term_q.displacement, qq ) );
00816         }
00817       } // for i
00818       //QDPIO::cout << __PRETTY_FUNCTION__ << ": exiting" << endl;
00819       END_CODE();
00820     } // void GroupBaryonQQQ::displaceQuarks
00821     
00822 
00823     //! First smear then displace the quarks
00824     void
00825     GroupBaryonQQQ::smearDisplaceQuarks( multi1d< map<int, LatticeFermion> >& disp_quarks,
00826                                          const LatticeFermion& q1,
00827                                          const LatticeFermion& q2,
00828                                          const LatticeFermion& q3,
00829                                                                                                                                                                  int* qindices
00830                                        ) const
00831     { 
00832                         START_CODE();
00833       //QDPIO::cout << __PRETTY_FUNCTION__ << ": entering" << endl;
00834       multi1d<LatticeFermion> q( 3 );
00835       // Rotate the quarks up-front.
00836       // NOTE: this step assumes the spin rotation commutes with the
00837       // quark smearing and the displacement
00838       // However, we are careful about the ordering of the smearing and
00839       // the displacement
00840 //      q[ 0 ] = rotateMat() * q1;
00841 //      q[ 1 ] = rotateMat() * q2;
00842 //      q[ 2 ] = rotateMat() * q3;
00843       q[ 0 ] = q1;
00844       q[ 1 ] = q2;
00845       q[ 2 ] = q3;
00846       // Displace after the smearing
00847       displaceQuarks( disp_quarks, q, qindices );
00848       //QDPIO::cout << __PRETTY_FUNCTION__ << ": exiting" << endl;
00849       END_CODE();
00850     } // void GroupBaryonQQQ::smearDisplaceQuarks
00851 
00852                 
00853                 //! hack to avoid compiler errors
00854                 LatticeComplex GroupBaryonQQQ::operator() ( const LatticeFermion& q1,
00855                                                             const LatticeFermion& q2,
00856                                                                                                                                                                                                 const LatticeFermion& q3
00857                                                                                                                                                                                         ) const
00858                 { 
00859                         START_CODE();
00860                         LatticeComplex d;
00861                         return d;
00862                         END_CODE();
00863                 }               
00864                 multi1d<LatticeComplex> GroupBaryonQQQ::operator() ( const LatticeFermion& q1,
00865                                                                      const LatticeFermion& q2,
00866                                                                                                                                                                                                                                  const LatticeFermion& q3,
00867                                                                                                                                                                                                                                  enum PlusMinus isign ) const
00868                 { 
00869                         START_CODE();
00870                         multi1d<LatticeComplex> d(1);
00871                         return d;
00872                         END_CODE();
00873                 }
00874     
00875     
00876     // for the creation operator                
00877     //! Compute the operator
00878     LatticeComplex 
00879     GroupBaryonQQQ::operator() ( const LatticeFermion& q1,
00880                                  const LatticeFermion& q2,
00881                                  const LatticeFermion& q3,
00882                                                                                                                                  int* qindices 
00883                                ) const
00884     { 
00885                         START_CODE();
00886                         if ( Nc != 3 ){    /* Code is specific to Ns=4 and Nc=3. */
00887                           QDPIO::cerr<<" code only works for Nc=3 and Ns=4\n";
00888                           QDP_abort(111) ;
00889                         }
00890 #if QDP_NC == 3
00891 
00892 
00893       // The result of displace and smearing (in some unspecified order here)
00894       multi1d< map<int, LatticeFermion> > disp_quarks;
00895       // Depending on whether this is the sink or source, do the 
00896       // appropriate combination of smearing and displacing
00897       smearDisplaceQuarks( disp_quarks, q1, q2, q3, qindices );
00898       // The return
00899       LatticeComplex d;
00900 
00901       LatticeColorVector c0 = peekSpin( disp_quarks[ 0 ].find( quark[ qindices[ 0 ] ].displacement ) ->second,
00902                                                                       quark[ qindices[ 0 ] ].spin );
00903       LatticeColorVector c1 = peekSpin( disp_quarks[ 1 ].find( quark[ qindices[ 1 ] ].displacement ) ->second,
00904                                                                       quark[ qindices[ 1 ] ].spin );
00905       LatticeColorVector c2 = peekSpin( disp_quarks[ 2 ].find( quark[ qindices[ 2 ] ].displacement ) ->second,
00906                                                                       quark[ qindices[ 2 ] ].spin );
00907 /*
00908       multi1d<LatticeFermion> q( 3 );
00909       // Rotate the quarks up-front.
00910       // NOTE: this step assumes the spin rotation commutes with the
00911       // quark smearing and the displacement
00912       // However, we are careful about the ordering of the smearing and
00913       // the displacement
00914       q[ 0 ] = rotateMat() * q1;
00915       q[ 1 ] = rotateMat() * q2;
00916       q[ 2 ] = rotateMat() * q3;
00917 */
00918       // Contract over color indices with antisym tensors
00919       d = colorContract( c0, c1, c2 );
00920 
00921       //QDPIO::cout << __PRETTY_FUNCTION__ << ": exiting" << endl;
00922       END_CODE();
00923       return d;
00924 
00925 #else
00926       LatticeComplex d;
00927       d = zero ; 
00928       return d;
00929 #endif
00930 
00931     } // LatticeComplex GroupBaryonQQQ::operator()
00932 
00933     // for the annihilation operator
00934     multi1d<LatticeComplex>
00935     GroupBaryonQQQ::operator() ( const LatticeFermion& q1,
00936                                  const LatticeFermion& q2,
00937                                  const LatticeFermion& q3,
00938                                                                                                                                  int* qindices,
00939                                  enum PlusMinus isign ) const
00940     { 
00941       START_CODE();
00942       if ( Nc != 3 ){    /* Code is specific to Ns=4 and Nc=3. */
00943         QDPIO::cerr<<" code only works for Nc=3 and Ns=4\n";
00944         QDP_abort(111) ;
00945       }
00946 #if QDP_NC == 3
00947 
00948 
00949       //QDPIO::cout << __PRETTY_FUNCTION__ << ": entering" << endl;
00950       // The result of displace and smearing (in some unspecified order here)
00951       multi1d< map<int, LatticeFermion> > disp_quarks;
00952       // Depending on whether this is the sink or source, do the 
00953       // appropriate combination of smearing and displacing
00954       smearDisplaceQuarks( disp_quarks, q1, q2, q3, qindices );
00955       // The return
00956       multi1d<LatticeComplex> d( 1 ); //( 6 )
00957       // Contract over color indices with antisym tensors
00958 
00959       LatticeColorVector c0 = peekSpin( disp_quarks[ 0 ].find( quark[ qindices[ 0 ] ].displacement ) ->second,
00960                                                                       quark[ qindices[ 0 ] ].spin );
00961       LatticeColorVector c1 = peekSpin( disp_quarks[ 1 ].find( quark[ qindices[ 1 ] ].displacement ) ->second,
00962                                                                       quark[ qindices[ 1 ] ].spin );
00963       LatticeColorVector c2 = peekSpin( disp_quarks[ 2 ].find( quark[ qindices[ 2 ] ].displacement ) ->second,
00964                                                                       quark[ qindices[ 2 ] ].spin );
00965 /*
00966       multi1d<LatticeFermion> q( 3 );
00967       // Rotate the quarks up-front.
00968       // NOTE: this step assumes the spin rotation commutes with the
00969       // quark smearing and the displacement
00970       // However, we are careful about the ordering of the smearing and
00971       // the displacement
00972       q[ 0 ] = rotateMat() * q1;
00973       q[ 1 ] = rotateMat() * q2;
00974       q[ 2 ] = rotateMat() * q3;
00975 */
00976       // Contract over color indices with antisym tensors
00977       d[ 0 ] = colorContract( c0, c1, c2 );
00978 
00979       //QDPIO::cout << __PRETTY_FUNCTION__ << ": exiting" << endl;
00980       END_CODE();
00981       return d;
00982 #else
00983 
00984       multi1d<LatticeComplex> d( 1 ); 
00985       return d;
00986 #endif
00987 
00988     } // multi1d<LatticeComplex> GroupBaryonQQQ::operator()
00989 
00990 
00991     //  ===========================
00992     //  GroupBaryonOp stuff
00993                 //      now only used to store data
00994     //  ===========================
00995         //! Serialize baryon_operator object
00996         multi1d<Complex> 
00997                 BaryonOperator_t::serialize()
00998         {
00999           int orderings_size = orderings.size();
01000           int op_size3 = orderings[ 0 ].op.size3();
01001           int op_size2 = orderings[ 0 ].op.size2();
01002           int op_size1 = orderings[ 0 ].op.size1();
01003           int ind_size = orderings[ 0 ].op( 0, 0, 0 ).ind.size();
01004           int elem_size2 = orderings[ 0 ].op( 0, 0, 0 ).ind[ 0 ].elem.size2();
01005           int elem_size1 = orderings[ 0 ].op( 0, 0, 0 ).ind[ 0 ].elem.size1();
01006           Complex ord_sizes, op_sizes1, op_sizes2, elem_sizes;
01007           ord_sizes = cmplx( Real( orderings_size ), Real( zero ) );
01008           op_sizes1 = cmplx( Real( op_size2 ), Real( op_size1 ) );
01009           op_sizes2 = cmplx( Real( op_size3 ), Real( ind_size ) );
01010           elem_sizes = cmplx( Real( elem_size2 ), Real( elem_size1 ) );
01011       multi1d<Complex> baryop_1d( 4 + orderings_size * op_size3 * op_size2 * op_size1 * ind_size * elem_size2 * elem_size1 );
01012           //QDPIO::cout << "baryop_size=" << baryop_1d.size() << endl; 
01013           int cnt = 0;
01014           baryop_1d[ cnt++ ] = ord_sizes;
01015           baryop_1d[ cnt++ ] = op_sizes1;
01016           baryop_1d[ cnt++ ] = op_sizes2;
01017           baryop_1d[ cnt++ ] = elem_sizes;
01018           for(int s=0; s < orderings.size(); ++s)              // orderings
01019           {
01020             for(int i=0; i < orderings[ s ].op.size3(); ++i)              // op_l
01021               for(int j=0; j < orderings[ s ].op.size2(); ++j)            // op_m
01022                 for(int k=0; k < orderings[ s ].op.size1(); ++k)          // op_r
01023                   for(int l=0; l < orderings[ s ].op( i, j, k ).ind.size(); ++l)     // ind
01024                     for(int a=0; a < orderings[ s ].op( i, j, k ).ind[ l ].elem.size2(); ++a)     // elem_l
01025                       for(int b=0; b < orderings[ s ].op( i, j, k ).ind[ l ].elem.size1(); ++b)   // elem_r
01026                         baryop_1d[ cnt++ ] = orderings[ s ].op( i, j, k ).ind[ l ].elem( a, b );
01027           }
01028           if ( cnt != baryop_1d.size() )
01029           {
01030             QDPIO::cerr << ": size mismatch in BaryonOperator_t serialization " << cnt << endl;
01031             QDP_abort( 1 );
01032           }
01033           return baryop_1d;
01034         }
01035                 
01036     //! Anonymous namespace
01037     namespace
01038     {
01039 
01040       //-------------------- callback functions ---------------------------------------
01041 
01042       //! Call-back function
01043       BaryonOperator<LatticeFermion>* groupBaryon( XMLReader& xml_in,
01044                                                    const std::string& path,
01045                                                    const multi1d<LatticeColorMatrix>& u )
01046       {
01047         return new GroupBaryonQQQ( Params( xml_in, path ) );
01048       }
01049       //! Local registration flag
01050       bool registered = false;
01051                         
01052     }  // end anonymous namespace
01053 
01054     const std::string name = "GROUP_BARYON";
01055 
01056     //! Register all the factories
01057     bool registerAll()
01058     {
01059       bool success = true;
01060       if ( ! registered )
01061       {
01062         // Required stuff
01063         success &= LinkSmearingEnv::registerAll();
01064                                 success &= QuarkSourceSmearingEnv::registerAll();
01065                                 success &= QuarkSinkSmearingEnv::registerAll();
01066 
01067         //! Register all the factories
01068         success &= Chroma::TheWilsonBaryonOperatorFactory::Instance().registerObject(name, groupBaryon );
01069 
01070         registered = true;
01071       }
01072       return success;
01073     } // registerAll()
01074 
01075 
01076   } // namespace GroupBaryonOperatorEnv
01077 
01078 
01079 }  // namespace Chroma
01080
Generated on Sun Nov 22 04:30:09 2009 for CHROMA by  doxygen 1.4.7