lemur::api::DOCID_T indri::index::MemoryIndex::addDocument( indri::api::ParsedDocument& document ) {
  indri::thread::ScopedLock sl( _writeLock );
  
  unsigned int position = 0;
  unsigned int extentIndex = 0;
  indri::utility::greedy_vector<indri::parse::TagExtent *> openTags;
  indri::utility::greedy_vector<indri::parse::TagExtent *> indexedTags;
  unsigned int indexedTerms = 0;
  indri::utility::greedy_vector<char*>& words = document.terms;
  term_entry* entries = 0;

  // assign a document ID
  lemur::api::DOCID_T documentID = _corpusStatistics.maximumDocument;
  _corpusStatistics.totalDocuments++;
  _corpusStatistics.maximumDocument++;
  
  _termList.clear();

  // move words into inverted lists, recording model statistics as we go
  for( position = 0; position < words.size(); position++ ) {
    const char* word = words[position];
    
    if( !word || *word == 0 ) {
      _termList.addTerm(0);
      continue;
    }

    int wordLength = strlen(word);

    if( wordLength >= lemur::file::Keyfile::MAX_KEY_LENGTH-1 ) {
      _termList.addTerm(0);     
      continue;
    }

    // fetch everything we know about this word so far
    term_entry* entry = _lookupTerm( word );

    // store information about this term location
    indri::index::TermData* termData = entry->termData;

    // store this term in the direct list
    _termList.addTerm( entry->termID );

    assert( entry->list.termFrequency() == termData->corpus.totalCount );

    // store this term in the inverted list
    if( !entry->marked() )
      entry->list.startDocument( documentID );
    entry->list.addLocation( position ); 
    termData->corpus.totalCount++;

    assert( entry->list.termFrequency() == termData->corpus.totalCount );

    // link this term_entry onto a list of ones we've seen
    if( entries == 0 ) {
      entry->mark();
      entries = entry;
    } else if( !entry->marked() ) {
      entry->next = entries;
      entries = entry;
    }

    // update our open tag knowledge
    _addOpenTags( indexedTags, openTags, document.tags, extentIndex, position );
    _removeClosedTags( openTags, position );

    // for every open tag, we want to record that we've seen the 
    for( indri::utility::greedy_vector<indri::parse::TagExtent *>::iterator tag = openTags.begin(); tag != openTags.end(); tag++ ) {
      int id = _fieldID( (*tag)->name );
      indri::index::TermFieldStatistics* termField = &entry->termData->fields[id - 1];
      termField->addOccurrence( documentID );

      indri::index::FieldStatistics* field = &_fieldData[id - 1];
      field->addOccurrence( documentID );
    }

    indexedTerms++;
  }

  _corpusStatistics.totalTerms += words.size();

  // need to add any tags that contain no text at the end of a document
  _addOpenTags( indexedTags, openTags, document.tags, extentIndex, position );
  _removeClosedTags( openTags, position );

  // go through the list of terms we've seen and update doc length counts
  term_entry* entry = entries;
  int uniqueTerms = 0;

  while( entry ) {
    indri::index::TermData* termData = entry->termData;
    term_entry* old = entry;

    termData->maxDocumentLength = lemur_compat::max<int>( termData->maxDocumentLength, words.size() );
    termData->minDocumentLength = lemur_compat::min<int>( termData->minDocumentLength, words.size() );
    termData->corpus.documentCount++;

    entry->list.endDocument();
    entry = entry->hasNext() ? entry->next : 0;
    old->clearMark();
//.........这里部分代码省略.........

void VertexBuffer::addBuffer(VertexBuffer *vb)
{
    ScopedLock sl(&_mutex);  
    findLast()->next = vb;
}

bool Track::IsMuted() const
{
   ScopedLock sl(fMutex);
   return fMuted;
}

 void RangeDeleterMockEnv::waitForNthGetCursor(uint64_t nthCall) {
     scoped_lock sl(_envStatMutex);
     while (_getCursorsCallCount < nthCall) {
         _cursorsCallCountUpdatedCV.wait(sl.boost());
     }
 }

  void SessionFactory::CreateSession(Node *node, const Id &session_id,
      SessionType type, AuthFactory::AuthType auth_type,
      const QSharedPointer<KeyShare> &public_keys)
  {
    CreateRound cr;
    switch(type) {
      case NULL_ROUND:
        cr = &TCreateRound<NullRound>;
        break;
      case SHUFFLE:
        cr = &TCreateRound<ShuffleRound>;
        break;
      case BULK:
        cr = &TCreateRound<BulkRound>;
        break;
      case REPEATING_BULK:
        cr = &TCreateRound<RepeatingBulkRound>;
        break;
      case CSBULK:
        cr = &TCreateBulkRound<CSBulkRound, NeffKeyShuffleRound>;
        break;
      case NEFF_SHUFFLE:
        cr = &TCreateRound<NeffShuffleRound>;
      case BLOG_DROP_REACTIVE:
        cr = TCreateBlogDropRound<Parameters::ParameterType_CppECHashingProduction,
           NeffShuffleRound, false>;
        break;
      case BLOG_DROP_PROACTIVE:
        cr = TCreateBlogDropRound<Parameters::ParameterType_CppECHashingProduction,
           NeffShuffleRound, true>;
        break;
      default:
        qFatal("Invalid session type");
    }

    QSharedPointer<IAuthenticate> authe(AuthFactory::CreateAuthenticate(
          node, auth_type, public_keys));

    Session *session = new Session(node->GetGroupHolder(), authe, session_id,
        node->GetNetwork(), cr);

    QObject::connect(node->GetOverlay().data(), SIGNAL(Disconnecting()),
        session, SLOT(CallStop()));

    QSharedPointer<Session> psession(session);
    session->SetSharedPointer(psession);
    node->GetSessionManager().AddSession(psession);

    psession->SetSink(node->GetSink().data());
    if(node->GetPrivateIdentity().GetLocalId() ==
        node->GetGroupHolder()->GetGroup().GetLeader())
    {
      QSharedPointer<IAuthenticator> autho(AuthFactory::CreateAuthenticator(
            node, auth_type, public_keys));
      QSharedPointer<SessionLeader> sl(new SessionLeader(
            node->GetGroupHolder()->GetGroup(), node->GetPrivateIdentity(),
            node->GetNetwork(), psession, autho));
      node->GetSessionManager().AddSessionLeader(sl);
      sl->Start();
    } else {
      psession->Start();
    }
  }

bool
SharedPortEndpoint::InitRemoteAddress()
{
		// Why do we read SharedPortServer's address from a file rather
		// than getting it passed down to us via the environment or
		// having a (configurable) fixed port?  Because the
		// SharedPortServer daemon may be listening via CCB, and its CCB
		// contact info may not be known as soon as it is started up
		// or may even change over time.

		// Why don't we just use a daemon client object to find the
		// address of the SharedPortServer daemon?  Because daemon
		// client assumes we want the best address for _us_ to connect
		// to.  That's not necessarily the public address that we want
		// to advertise for others to connect to.

	MyString shared_port_server_ad_file;
	if( !param(shared_port_server_ad_file,"SHARED_PORT_DAEMON_AD_FILE") ) {
		EXCEPT("SHARED_PORT_DAEMON_AD_FILE must be defined");
	}

	FILE *fp = safe_fopen_wrapper_follow(shared_port_server_ad_file.Value(),"r");
	if( !fp ) {
		dprintf(D_ALWAYS,"SharedPortEndpoint: failed to open %s: %s\n",
				shared_port_server_ad_file.Value(), strerror(errno));
		return false;
	}

	int adIsEOF = 0, errorReadingAd = 0, adEmpty = 0;
	ClassAd *ad = new ClassAd(fp, "[classad-delimiter]", adIsEOF, errorReadingAd, adEmpty);
	ASSERT(ad);
	fclose( fp );

		// avoid leaking ad when returning from this function
	counted_ptr<ClassAd> smart_ad_ptr(ad);

	if( errorReadingAd ) {
		dprintf(D_ALWAYS,"SharedPortEndpoint: failed to read ad from %s.\n",
				shared_port_server_ad_file.Value());
		return false;
	}

	MyString public_addr;
	if( !ad->LookupString(ATTR_MY_ADDRESS,public_addr) ) {
		dprintf(D_ALWAYS,
				"SharedPortEndpoint: failed to find %s in ad from %s.\n",
				ATTR_MY_ADDRESS, shared_port_server_ad_file.Value());
		return false;
	}

	Sinful sinful(public_addr.Value());
	sinful.setSharedPortID( m_local_id.Value() );

		// if there is a private address, set the shared port id on that too
	char const *private_addr = sinful.getPrivateAddr();
	if( private_addr ) {
		Sinful private_sinful( private_addr );
		private_sinful.setSharedPortID( m_local_id.Value() );
		sinful.setPrivateAddr( private_sinful.getSinful() );
	}

	// Next, look for alternate command strings
	std::string commandStrings;
	if (ad->EvaluateAttrString(ATTR_SHARED_PORT_COMMAND_SINFULS, commandStrings))
	{
		m_remote_addrs.clear();
		StringList sl(commandStrings.c_str());
		sl.rewind();
		const char *commandSinfulStr;
		while ((commandSinfulStr = sl.next()))
		{
			Sinful altsinful(commandSinfulStr);
			altsinful.setSharedPortID(m_local_id.Value());
			char const *private_addr = sinful.getPrivateAddr();
			if (private_addr)
			{
				Sinful private_sinful(private_addr);
				private_sinful.setSharedPortID(m_local_id.Value());
				altsinful.setPrivateAddr(private_sinful.getSinful());
			}
			m_remote_addrs.push_back(altsinful);
		}
	}

	m_remote_addr = sinful.getSinful();

	return true;
}

m64p_error open_rom(const unsigned char* romimage, unsigned int size)
{
    md5_state_t state;
    md5_byte_t digest[16];
    romdatabase_entry* entry;
    char buffer[256];
    unsigned char imagetype;
    int i;

    /* check input requirements */
    if (rom != NULL)
    {
        DebugMessage(M64MSG_ERROR, "open_rom(): previous ROM image was not freed");
        return M64ERR_INTERNAL;
    }
    if (romimage == NULL || !is_valid_rom(romimage))
    {
        DebugMessage(M64MSG_ERROR, "open_rom(): not a valid ROM image");
        return M64ERR_INPUT_INVALID;
    }

    /* Clear Byte-swapped flag, since ROM is now deleted. */
    g_MemHasBeenBSwapped = 0;
    /* allocate new buffer for ROM and copy into this buffer */
    rom_size = size;
    rom = (unsigned char *) malloc(size);
    if (rom == NULL)
        return M64ERR_NO_MEMORY;
    memcpy(rom, romimage, size);
    swap_rom(rom, &imagetype, rom_size);

    memcpy(&ROM_HEADER, rom, sizeof(m64p_rom_header));

    /* Calculate MD5 hash  */
    md5_init(&state);
    md5_append(&state, (const md5_byte_t*)rom, rom_size);
    md5_finish(&state, digest);
    for ( i = 0; i < 16; ++i )
        sprintf(buffer+i*2, "%02X", digest[i]);
    buffer[32] = '\0';
    strcpy(ROM_SETTINGS.MD5, buffer);

    /* add some useful properties to ROM_PARAMS */
    ROM_PARAMS.systemtype = rom_country_code_to_system_type(ROM_HEADER.Country_code);
    ROM_PARAMS.vilimit = rom_system_type_to_vi_limit(ROM_PARAMS.systemtype);
    ROM_PARAMS.aidacrate = rom_system_type_to_ai_dac_rate(ROM_PARAMS.systemtype);

    memcpy(ROM_PARAMS.headername, ROM_HEADER.Name, 20);
    ROM_PARAMS.headername[20] = '\0';
    trim(ROM_PARAMS.headername); /* Remove trailing whitespace from ROM name. */

    /* Look up this ROM in the .ini file and fill in goodname, etc */
    if ((entry=ini_search_by_md5(digest)) != NULL ||
        (entry=ini_search_by_crc(sl(ROM_HEADER.CRC1),sl(ROM_HEADER.CRC2))) != NULL)
    {
        strncpy(ROM_SETTINGS.goodname, entry->goodname, 255);
        ROM_SETTINGS.goodname[255] = '\0';
        ROM_SETTINGS.savetype = entry->savetype;
        ROM_SETTINGS.status = entry->status;
        ROM_SETTINGS.players = entry->players;
        ROM_SETTINGS.rumble = entry->rumble;
    }
    else
    {
        strcpy(ROM_SETTINGS.goodname, ROM_PARAMS.headername);
        strcat(ROM_SETTINGS.goodname, " (unknown rom)");
        ROM_SETTINGS.savetype = NONE;
        ROM_SETTINGS.status = 0;
        ROM_SETTINGS.players = 0;
        ROM_SETTINGS.rumble = 0;
    }

    /* print out a bunch of info about the ROM */
    DebugMessage(M64MSG_INFO, "Goodname: %s", ROM_SETTINGS.goodname);
    DebugMessage(M64MSG_INFO, "Name: %s", ROM_HEADER.Name);
    imagestring(imagetype, buffer);
    DebugMessage(M64MSG_INFO, "MD5: %s", ROM_SETTINGS.MD5);
    DebugMessage(M64MSG_INFO, "CRC: %x %x", sl(ROM_HEADER.CRC1), sl(ROM_HEADER.CRC2));
    DebugMessage(M64MSG_INFO, "Imagetype: %s", buffer);
    DebugMessage(M64MSG_INFO, "Rom size: %d bytes (or %d Mb or %d Megabits)", rom_size, rom_size/1024/1024, rom_size/1024/1024*8);
    DebugMessage(M64MSG_VERBOSE, "ClockRate = %x", sl(ROM_HEADER.ClockRate));
    DebugMessage(M64MSG_INFO, "Version: %x", sl(ROM_HEADER.Release));
    if(sl(ROM_HEADER.Manufacturer_ID) == 'N')
        DebugMessage(M64MSG_INFO, "Manufacturer: Nintendo");
    else
        DebugMessage(M64MSG_INFO, "Manufacturer: %x", sl(ROM_HEADER.Manufacturer_ID));
    DebugMessage(M64MSG_VERBOSE, "Cartridge_ID: %x", ROM_HEADER.Cartridge_ID);
    countrycodestring(ROM_HEADER.Country_code, buffer);
    DebugMessage(M64MSG_INFO, "Country: %s", buffer);
    DebugMessage(M64MSG_VERBOSE, "PC = %x", sl((unsigned int)ROM_HEADER.PC));
    DebugMessage(M64MSG_VERBOSE, "Save type: %d", ROM_SETTINGS.savetype);

    //Prepare Hack for GOLDENEYE
    isGoldeneyeRom = 0;
    if(strcmp(ROM_PARAMS.headername, "GOLDENEYE") == 0)
       isGoldeneyeRom = 1;

    return M64ERR_SUCCESS;
}

void CheatStart(eCheatMode CheatMode, char *CheatNumList)
{
    /* if cheat codes are disabled, then we don't have to do anything */
    if (CheatMode == CHEAT_DISABLE || (CheatMode == CHEAT_LIST && strlen(CheatNumList) == 0))
    {
        DebugMessage(M64MSG_STATUS, "Cheat codes disabled.");
        return;
    }

    /* get the ROM header for the currently loaded ROM image from the core */
    if ((*CoreDoCommand)(M64CMD_ROM_GET_HEADER, sizeof(l_RomHeader), &l_RomHeader) != M64ERR_SUCCESS)
    {
        DebugMessage(M64MSG_WARNING, "couldn't get ROM header information from core library");
        return;
    }

    /* generate section name from ROM's CRC and country code */
    char RomSection[24];
    sprintf(RomSection, "%08X-%08X-C:%X", sl(l_RomHeader.CRC1), sl(l_RomHeader.CRC2), l_RomHeader.Country_code & 0xff);

    /* parse through the cheat INI file and load up any cheat codes found for this ROM */
    ReadCheats(RomSection);
    if (!l_RomFound || l_CheatCodesFound == 0)
    {
        DebugMessage(M64MSG_WARNING, "no cheat codes found for ROM image '%.20s'", l_RomHeader.Name);
        CheatFreeAll();
        return;
    }

    /* handle the list command */
    if (CheatMode == CHEAT_SHOW_LIST)
    {
        DebugMessage(M64MSG_INFO, "%i cheat code(s) found for ROM '%s'", l_CheatCodesFound, l_CheatGameName);
        sCheatInfo *pCur = l_CheatList;
        while (pCur != NULL)
        {
            if (pCur->Description == NULL)
                DebugMessage(M64MSG_INFO, "   %i: %s", pCur->Number, pCur->Name);
            else
                DebugMessage(M64MSG_INFO, "   %i: %s (%s)", pCur->Number, pCur->Name, pCur->Description);
            if(pCur->VariableLine != -1)
            {
                int i;
                for (i = 0; i < pCur->Codes[pCur->VariableLine].var_count; i++)
                    DebugMessage(M64MSG_INFO, "      %i: %s", i, pCur->Codes[pCur->VariableLine].variable_names[i]);
            }
            pCur = pCur->Next;
        }
        CheatFreeAll();
        return;
    }

    /* handle all cheats enabled mode */
    if (CheatMode == CHEAT_ALL)
    {
        sCheatInfo *pCur = l_CheatList;
        while (pCur != NULL)
        {
            CheatActivate(pCur);
            pCur = pCur->Next;
        }
        CheatFreeAll();
        return;
    }

    /* handle list of cheats enabled mode */
    if (CheatMode == CHEAT_LIST)
    {
        int option, number;
        char *cheat_next;
        sCheatInfo *pCheat;
        while(CheatNumList != NULL && *CheatNumList)
        {
            if ((cheat_next = strchr(CheatNumList, ',')) != NULL)
            {
                *cheat_next = 0;
                cheat_next ++;
            }

            if (strchr(CheatNumList, '-') != NULL)
            {
                sscanf(CheatNumList, "%i-%i", &number, &option); /* option */
            }
            else
            {
                option=0;
                sscanf(CheatNumList, "%i", &number);
            }

            pCheat = CheatFindCode(number);
            if (pCheat == NULL)
                DebugMessage(M64MSG_WARNING, "invalid cheat code number %i", number);
            else
            {
                if (pCheat->VariableLine != -1 && pCheat->Count > pCheat->VariableLine && option < pCheat->Codes[pCheat->VariableLine].var_count)
                    pCheat->Codes[pCheat->VariableLine].var_to_use = option;
                CheatActivate(pCheat);
            }

            CheatNumList = cheat_next;
//.........这里部分代码省略.........

static void
sl (void)
{
  FILE *f = fopen (pidfile, "w");
  if (f == NULL)
    exit (1);

  fprintf (f, "%lld\n", (long long) getpid ());
  fflush (f);

  struct flock fl =
    {
      .l_type = F_WRLCK,
      .l_start = 0,
      .l_whence = SEEK_SET,
      .l_len = 1
    };
  if (fcntl (fileno (f), F_SETLK, &fl) != 0)
    exit (1);

  sigset_t ss;
  sigfillset (&ss);
  sigsuspend (&ss);
  exit (0);
}


static void
do_prepare (int argc, char *argv[])
{
  if (command == NULL)
    command = argv[0];

  if (pidfile)
    sl ();

  int fd = mkstemp (pidfilename);
  if (fd == -1)
    {
      puts ("mkstemp failed");
      exit (1);
    }

  write (fd, " ", 1);
  close (fd);
}


static int
do_test (void)
{
  pthread_t th;
  if (pthread_create (&th, NULL, tf, NULL) != 0)
    {
      puts ("pthread_create failed");
      return 1;
    }

  do
    sleep (1);
  while (access (pidfilename, R_OK) != 0);

  if (pthread_cancel (th) != 0)
    {
      puts ("pthread_cancel failed");
      return 1;
    }

  void *r;
  if (pthread_join (th, &r) != 0)
    {
      puts ("pthread_join failed");
      return 1;
    }

  sleep (1);

  FILE *f = fopen (pidfilename, "r+");
  if (f == NULL)
    {
      puts ("no pidfile");
      return 1;
    }

  long long ll;
  if (fscanf (f, "%lld\n", &ll) != 1)
    {
      puts ("could not read pid");
      unlink (pidfilename);
      return 1;
    }

  struct flock fl =
    {
      .l_type = F_WRLCK,
      .l_start = 0,
      .l_whence = SEEK_SET,
      .l_len = 1
    };
  if (fcntl (fileno (f), F_GETLK, &fl) != 0)
//.........这里部分代码省略.........

bool SmugMug::isUploading()
{
	ScopedLock sl(lock);
	return uploadThreads.size() > 0 || dupeThreads.size() > 0;
}

// This is what does the work in the background thread
UINT CHexEditDoc::RunAerialThread()
{
    // Keep looping until we are told to die
    for (;;)
    {
        // Signal that we are waiting then wait for start_aerial_event_ to be pulsed
        {
            CSingleLock sl(&docdata_, TRUE);
            aerial_state_ = WAITING;
        }
        TRACE1("+++ BGAerial: waiting for %p\n", this);
        DWORD wait_status = ::WaitForSingleObject(HANDLE(start_aerial_event_), INFINITE);
        docdata_.Lock();
        aerial_state_ = SCANNING;
        docdata_.Unlock();
        start_aerial_event_.ResetEvent();
        ASSERT(wait_status == WAIT_OBJECT_0);
        TRACE1("+++ BGAerial: got event for %p\n", this);

        if (AerialProcessStop())
            continue;

        // Reset for new scan
        docdata_.Lock();
        aerial_fin_ = false;
        aerial_addr_ = 0;
        FILE_ADDRESS file_len = length_;
        int file_bpe = bpe_;
        unsigned char *file_dib = FreeImage_GetBits(dib_);
        unsigned dib_size = FreeImage_GetDIBSize(dib_);
        docdata_.Unlock();
        TRACE("+++ BGAerial: using bitmap at %p\n", file_dib);

        // Get the file buffer
        size_t buf_len = (size_t)min(file_len, 65536);
        ASSERT(aerial_buf_ == NULL);
        aerial_buf_ = new unsigned char[buf_len];

        for (;;)
        {
            // First check if we need to stop
            if (AerialProcessStop())
                break;   // stop processing and go back to WAITING state

            // Check if we have finished scanning the file
            if (aerial_addr_ >= file_len)
            {
                TRACE2("+++ BGAerial: finished scan for %p at address %p\n", this, file_dib + 3*size_t(aerial_addr_/file_bpe));
                CSingleLock sl(&docdata_, TRUE); // Protect shared data access

                aerial_fin_ = true;
                break;                          // falls out to end_scan
            }

            // Get the next buffer full from the file and scan it
            size_t got = GetData(aerial_buf_, buf_len, aerial_addr_, 3);
            ASSERT(got <= buf_len);

            unsigned char *pbm = file_dib + 3*size_t(aerial_addr_/file_bpe);    // where we write to bitmap
            unsigned char *pbuf;                                        // where we read from the file buffer
            for (pbuf = aerial_buf_; pbuf < aerial_buf_ + got; pbuf += file_bpe, pbm += 3)
            {
                int r, g, b;
                r = g = b = 0;
                for (unsigned char *pp = pbuf; pp < pbuf + file_bpe; ++pp)
                {
                    r += GetRValue(kala_[*pp]);
                    g += GetGValue(kala_[*pp]);
                    b += GetBValue(kala_[*pp]);
                }
                *pbm     = unsigned char(b/file_bpe);
                *(pbm+1) = unsigned char(g/file_bpe);
                *(pbm+2) = unsigned char(r/file_bpe);
            }
            aerial_addr_ += got;
        }

        delete[] aerial_buf_;
        aerial_buf_ = NULL;
    }
    return 0;   // never reached
}

int main(int argc, char** argv)
{

    int ppw=10;     // Point per wavelength
    std::string filename="kitenormcond10.txt";


    std::vector<double> freqs;
    freqs.push_back(5);
    freqs.push_back(10);
    freqs.push_back(20);
    freqs.push_back(40);
    freqs.push_back(80);
    freqs.push_back(160);
    freqs.push_back(320);
    freqs.push_back(640);



    std::vector<double> norm_sl(freqs.size());
    std::vector<double> norm_dl(freqs.size());
    std::vector<double> norm_combined1(freqs.size());
    std::vector<double> norm_combined2(freqs.size());
    std::vector<double> cond_sl(freqs.size());
    std::vector<double> cond_dl(freqs.size());
    std::vector<double> cond_combined1(freqs.size());
    std::vector<double> cond_combined2(freqs.size());



    clock_t start, finish;
    double time;
    start=clock();


#ifdef BEM2DMPI
    MPI_Init(&argc, &argv);


    int nprow=4; // Number of rows in process grid
    int npcol=2; // Number of columns in process grid
    int mb=24;  // Row Block size
    int nb=24;  // Column Block size
    bem2d::BlacsSystem* b=bem2d::BlacsSystem::Initialize(nprow,npcol,mb,nb);

    // Exit if Context could not be created or process does not belong to context

    if (!b) {
        std::cout <<  "Could not create Blacs context" << std::endl;
        MPI_Finalize();
        exit(1);
    }
    if ((b->get_myrow()==-1)&&(b->get_mycol()==-1)) {
        MPI_Finalize();
        exit(0);
    }
#endif

    for (int j=0; j<freqs.size(); j++) {

        bem2d::freqtype k= {(double)freqs[j],0};
        double eta1=k.re; // Coupling between conj. double and single layer pot.
        double eta2=cbrt(k.re*k.re);
        bem2d::pCurve kobj(new bem2d::Kite);
        int n=(int)(kobj->Length()*k.re*ppw/2.0/bem2d::PI);
        bem2d::AnalyticCurve kite(n,kobj);
        bem2d::pGeometry pgeom=kite.GetGeometry();

        bem2d::PolBasis::AddBasis(0,pgeom); // Add constant basis functions


        // Discretize the single and double layer potential

        bem2d::SingleLayer sl(k);
        bem2d::ConjDoubleLayer cdl(k);
        bem2d::DoubleLayer dl(k);

        bem2d::QuadOption quadopts;

        quadopts.L=3;
        quadopts.N=5;
        quadopts.sigma=0.15;

#ifdef BEM2DMPI
        if (b->IsRoot()) {
            std::cout << "Discretize Kernels with n=" << n << std::endl;
        }
#else
        std::cout << "Discretize Kernels with n=" << n << std::endl;
#endif



        bem2d::Matrix dsl=*(DiscreteKernel(*pgeom,quadopts,sl));
        bem2d::Matrix ddl=(*DiscreteKernel(*pgeom,quadopts,dl));
        bem2d::Matrix dcdl=*(DiscreteKernel(*pgeom,quadopts,cdl));
        bem2d::Matrix Id=*(EvalIdent(*pgeom, quadopts));
        bem2d::Matrix combined1=Id+2.0*dcdl-bem2d::complex(0,2.0)*eta1*dsl;
        bem2d::Matrix combined2=Id+2.0*dcdl-bem2d::complex(0,2.0)*eta2*dsl;

//.........这里部分代码省略.........