自学教程：C++ z函数代码示例

double leekesler::hdep(){    double tr = T/Tcr;    return tr*tr*I() + (1.0 - z())*tr;}

GLdouble GLVector::operator *(const GLVector &v )const{    return  x() * v.x() + y() * v.y() + z() * v.z();}

bool GLVector::operator == (const GLVector &v) const{ return ((v.x() == x()) & (v.y() == y()) & (v.z() == z()));}

bool testSymbolSpace() {    bool passed = true;    symbol::Space<int> point;    symbol::Symbol x("x");    symbol::Symbol y("y");    symbol::Symbol z("z"); // reserved for "misses"    // nothing to get from empty space    passed &= (point.get(x) == NULL);    passed &= (point.get(y) == NULL);    passed &= (point.get(z) == NULL);    point.del(x);    point.del(y);    point.del(z);    // one key    point.set(x, 42);    passed &= (*point.get(x) == 42);    passed &= (point.get(z) == NULL);    point.del(z);     point.set(x, 99);    passed &= (*point.get(x) == 99);    point.del(x);    passed &= (point.get(x) == NULL);    point.del(z);     // two keys    point.set(x, 1);    point.set(y, 2); // append    passed &= (*point.get(x) == 1);    passed &= (*point.get(y) == 2);    passed &= (point.get(z) == NULL);    point.del(z);     point.del(x);    passed &= (point.get(x) == NULL);    passed &= (*point.get(y) == 2);    passed &= (point.get(z) == NULL);    point.del(z);     point.del(y);    passed &= (point.get(x) == NULL);    passed &= (point.get(y) == NULL);    passed &= (point.get(z) == NULL);    // two keys, again, different order    point.set(y, 2);    point.set(x, 1); // insert before head    passed &= (*point.get(x) == 1);    passed &= (*point.get(y) == 2);    passed &= (point.get(z) == NULL);    point.del(z);     point.del(x);    passed &= (point.get(x) == NULL);    passed &= (*point.get(y) == 2);    passed &= (point.get(z) == NULL);    point.del(z);     point.del(y);    passed &= (point.get(x) == NULL);    passed &= (point.get(y) == NULL);    passed &= (point.get(z) == NULL);    // two keys, with overwritting this time    point.set(y, 2);    point.set(x, 1);    passed &= (*point.get(x) == 1);    passed &= (*point.get(y) == 2);    point.set(x, 3);    passed &= (*point.get(x) == 3);    passed &= (*point.get(y) == 2);    point.set(y, 4);    passed &= (*point.get(x) == 3);    passed &= (*point.get(y) == 4);    point.del(x);    point.del(y);    passed &= (point.get(x) == NULL);    passed &= (point.get(y) == NULL);    passed &= (point.get(z) == NULL);    // three keys    point.set(x, 1);    point.set(z, 3);    point.set(y, 2); // insert in the middle!    passed &= (*point.get(x) == 1);    passed &= (*point.get(y) == 2);    passed &= (*point.get(z) == 3);    // overwrite at head, middle, and end    point.set(x, 4);    point.set(y, 5);     point.set(z, 6);    passed &= (*point.get(x) == 4);    passed &= (*point.get(y) == 5);    passed &= (*point.get(z) == 6);    // del from middle    point.del(y);    passed &= (*point.get(x) == 4);    passed &= (point.get(y) == NULL);    passed &= (*point.get(z) == 6);    symbol::Space<int> letters;//.........这里部分代码省略.........

const GLVector GLVector::vectorMult( const GLVector& v)const{ return GLVector(   y() * v.z() - z() * v.y(),   z() * v.x() - x() * v.z(),   x() * v.y() - y() * v.x());}

void MVertex::writePLY2(FILE *fp){  if(_index < 0) return; // negative index vertices are never saved  fprintf(fp, "%.16g %.16g %.16g/n", x(), y(), z());}

  void CollocationIntegrator::setupFG() {    // Interpolation order    deg_ = getOption("interpolation_order");    // All collocation time points    std::vector<long double> tau_root = collocationPointsL(deg_, getOption("collocation_scheme"));    // Coefficients of the collocation equation    vector<vector<double> > C(deg_+1, vector<double>(deg_+1, 0));    // Coefficients of the continuity equation    vector<double> D(deg_+1, 0);    // Coefficients of the quadratures    vector<double> B(deg_+1, 0);    // For all collocation points    for (int j=0; j<deg_+1; ++j) {      // Construct Lagrange polynomials to get the polynomial basis at the collocation point      Polynomial p = 1;      for (int r=0; r<deg_+1; ++r) {        if (r!=j) {          p *= Polynomial(-tau_root[r], 1)/(tau_root[j]-tau_root[r]);        }      }      // Evaluate the polynomial at the final time to get the      // coefficients of the continuity equation      D[j] = zeroIfSmall(p(1.0L));      // Evaluate the time derivative of the polynomial at all collocation points to      // get the coefficients of the continuity equation      Polynomial dp = p.derivative();      for (int r=0; r<deg_+1; ++r) {        C[j][r] = zeroIfSmall(dp(tau_root[r]));      }      // Integrate polynomial to get the coefficients of the quadratures      Polynomial ip = p.anti_derivative();      B[j] = zeroIfSmall(ip(1.0L));    }    // Symbolic inputs    MX x0 = MX::sym("x0", f_.input(DAE_X).sparsity());    MX p = MX::sym("p", f_.input(DAE_P).sparsity());    MX t = MX::sym("t", f_.input(DAE_T).sparsity());    // Implicitly defined variables (z and x)    MX v = MX::sym("v", deg_*(nx_+nz_));    vector<int> v_offset(1, 0);    for (int d=0; d<deg_; ++d) {      v_offset.push_back(v_offset.back()+nx_);      v_offset.push_back(v_offset.back()+nz_);    }    vector<MX> vv = vertsplit(v, v_offset);    vector<MX>::const_iterator vv_it = vv.begin();    // Collocated states    vector<MX> x(deg_+1), z(deg_+1);    for (int d=1; d<=deg_; ++d) {      x[d] = reshape(*vv_it++, this->x0().shape());      z[d] = reshape(*vv_it++, this->z0().shape());    }    casadi_assert(vv_it==vv.end());    // Collocation time points    vector<MX> tt(deg_+1);    for (int d=0; d<=deg_; ++d) {      tt[d] = t + h_*tau_root[d];    }    // Equations that implicitly define v    vector<MX> eq;    // Quadratures    MX qf = MX::zeros(f_.output(DAE_QUAD).sparsity());    // End state    MX xf = D[0]*x0;    // For all collocation points    for (int j=1; j<deg_+1; ++j) {      //for (int j=deg_; j>=1; --j) {      // Evaluate the DAE      vector<MX> f_arg(DAE_NUM_IN);      f_arg[DAE_T] = tt[j];      f_arg[DAE_P] = p;      f_arg[DAE_X] = x[j];      f_arg[DAE_Z] = z[j];      vector<MX> f_res = f_.call(f_arg);      // Get an expression for the state derivative at the collocation point      MX xp_j = C[0][j] * x0;      for (int r=1; r<deg_+1; ++r) {        xp_j += C[r][j] * x[r];      }//.........这里部分代码省略.........

LOCAL_C void testFix(CArrayFix<TBuf<0x10> >& aFix)//// Test all methods//	{	test.Next(_L("Test all methods"));	test(aFix.Count()==0);	test(aFix.Length()==sizeof(TBuf<0x10>));	aFix.Compress();	test(TRUE);	aFix.Reset();	test(TRUE);	TKeyArrayFix kk(0,ECmpNormal,0x10);	test(TRUE);	aFix.Sort(kk);	test(TRUE);	TBuf<0x10> aa(_L("aaaaa"));	aFix.InsertL(0,aa);	test(TRUE);	aFix[0].Fill(' ');	test(TRUE);	TBuf<0x10> z(aFix[0]);    z.Length();	test(TRUE);	aFix[0].Fill('a');	test(TRUE);	TInt pp;	test(aFix.Find(aa,kk,pp)==0);	test(pp==0);	aFix.Delete(0);	TBuf<0x10> bb(_L("bbbbb"));	aFix.AppendL(bb);	test(aFix.Count()==1);	test(aFix.InsertIsqAllowDuplicatesL(aa,kk)==0);	test(aFix.InsertIsqAllowDuplicatesL(bb,kk)==2);	test(aFix.FindIsq(aa,kk,pp)==0);	test(pp==0);	aFix.Reset();	for(TInt index=0;index<KTestGranularity*7/2;index++)		aFix.AppendL(aa);	const TBuf<0x10> *end=NULL;	const TBuf<0x10> *ptr=NULL;	for(TInt index2=0;index2<KTestGranularity*7/2;index2++)		{		if (end==ptr)			{			end=aFix.End(index2);			ptr=&aFix[index2];			TInt seglen=end-ptr;			test(seglen==KTestGranularity || seglen==(aFix.Count()-index2));			}		test(&aFix[index2]==ptr++);		}	const TBuf<0x10> *bak=NULL;	ptr=NULL;	for(TInt index3=KTestGranularity*7/2;index3>0;index3--)		{		if (bak==ptr)			{			bak=aFix.Back(index3);			ptr=&aFix[index3-1]+1;			TInt seglen=ptr-bak;			test(seglen==KTestGranularity || seglen==index3 || seglen==index3%KTestGranularity);			}		test(&aFix[index3-1]==--ptr);		}		//Test ExpandL	//Expand array in slot 1	TBuf16<0x10> exp;	exp=_L("abc AbC");	aFix.InsertL(0,exp);	aFix.InsertL(1,exp);	aFix.InsertL(2,exp);	exp=aFix.ExpandL(1);	test(aFix[0]==_L("abc AbC"));	test(aFix[1]==_L(""));	test(aFix[2]==_L("abc AbC"));	test(aFix[3]==_L("abc AbC"));		//Test ResizeL and InsertReplL	//Resize the array to containing 20 records,	//copying a record into any new slots.	TBuf<0x10> res(_L("bbbbb"));	aFix.Reset();	aFix.ResizeL(20,res);	for(TInt i=0;i<20;i++)	    {        test(aFix[1]==_L("bbbbb"));	    }	}

void dninst_(int *n,         int *nev,    double *sigmar,             double *sigmai, int *colptr, int *rowind,              double *nzvals, double *dr,  double *di,               double *z,      int *ldz,    int *info,     double *ptol)/*  Arguement list:    n       (int*)    Dimension of the problem.  (INPUT)    nev     (int*)    Number of eigenvalues requested.  (INPUT/OUTPUT)                      This routine is used to compute NEV eigenvalues                      nearest to a shift (sigmar, sigmai).                      On return, it gives the number of converged                      eigenvalues.    sigmar  (double*) Real part of the shift. (INPUT)    sigmai  (double*) Imaginar part of the shift. (INPUT)    colptr  (int*)    dimension n+1. (INPUT)                      Column pointers for the sparse matrix.    rowind  (int*)    dimension colptr[*n]-1. (INPUT)                      Row indices for the sparse matrix.    nzvals  (double*) dimension colptr[*n]-1. (INPUT)                      Nonzero values of the sparse matrix.                      The sparse matrix is represented by                      the above three arrays colptr, rowind, nzvals.    dr      (double*) dimension nev+1.  (OUTPUT)                      Real part of the eigenvalue.    di      (double*) dimension nev+1.  (OUTPUT)                      Imaginar part of the eigenvalue.    z       (double*) dimension ldz by nev+1. (OUTPUT)                      Eigenvector matrix.                      If the j-th eigenvalue is real, the j-th column                      of z contains the corresponding eigenvector.                      If the j-th and j+1st eigenvalues form a complex                      conjuagate pair, then the j-th column of z contains                      the real part of the eigenvector, and the j+1st column                      of z contains the imaginary part of the eigenvector.    ldz      (int*)   The leading dimension of z. (INPUT)    info     (int*)   Error flag to indicate whether the eigenvalues                       calculation is successful. (OUTPUT)                      *info = 0, successful exit                      *info = 1, Maximum number of iteration is reached                                 before all requested eigenvalues                                  have converged.*/{    int    i, j, ibegin, iend, ncv, neqns, token, order=2;    int    lworkl, ldv,  nnz, ione = 1;    double tol=1.0e-10, zero = 0.0;    double *workl, *workd, *resid, *workev, *v, *ax;    double numr, numi, denr, deni;    int    *select, first;    int    ido, ishfts, maxitr, mode, rvec, ierr1, ierr2;    int    iparam[11], ipntr[14];    char   *which="LM", bmat[2], *all="A";#ifdef USE_COMPLEX    doublecomplex *cvals, *cx, *crhs;#endif    neqns = *n;    *info = 0;    tol   = *ptol;    if ( tol < 1.0e-10 ) tol = 1.0e-10;    if ( tol > 1.0e-1  ) tol = 1.0e-1;    if (*n - *nev < 2) {       *info = -1000;       fprintf(stderr, " NEV must be less than N-2!/n");       goto Error_handle;     }    /* set parameters and allocate temp space for ARPACK*/    ncv = max(*nev+20, 2*(*nev));    if (ncv > neqns) ncv = neqns;    /* Convert from 1-based index to 0-based index */    nnz = colptr[neqns]-1;    for (j=0;j<=neqns;j++) colptr[j]--;    for (i=0;i<nnz;i++) rowind[i]--;    /* Subtract shift from the matrix */       if ( *sigmai == 0.0) {       /* real shift */       for (j = 0; j<neqns; j++) {          ibegin = colptr[j];          iend   = colptr[j+1]-1;          for (i=ibegin;i<=iend;i++)              if (j == rowind[i]) nzvals[i] = nzvals[i] - *sigmar;       }//.........这里部分代码省略.........

int main(int argc, char* argv[]){#ifdef _DIST_    CnC::dist_cnc_init< my_context > dc_init;#endif    bool verbose = false;    int max_row = 100;    int max_col = 100;    int max_depth = 10000;    int ai = 1;    if (argc > ai && 0 == strcmp("-v", argv[ai]))    {        verbose = true;        ai++;    }    if (argc == ai+3)    {        max_row = atoi(argv[ai]);        max_col = atoi(argv[ai+1]);        // you can't change a global variable when using distributed memory,        // so this is in the context and must be set before the first tag/item is put        max_depth = atoi(argv[ai+2]);     }    else    {        fprintf(stderr,"Usage: mandel [-v] rows columns max_depth/n");        return -1;    }    complex z(1.0,1.5);    std::cout << mandel(z, max_depth) << std::endl;        double r_origin = -2;    double r_scale = 4.0/max_row;    double c_origin = -2.0;    double c_scale = 4.0/max_col;    int *pixels = new int[max_row*max_col];    // set max-depth in the constructor    my_context c( max_depth );    tbb::tick_count t0 = tbb::tick_count::now();        for (int i = 0; i < max_row; i++)     {        for (int j = 0; j < max_col; j++ )        {            complex z = complex(r_scale*j +r_origin,c_scale*i + c_origin);            c.m_data.put(pair(i,j),z);            c.m_position.put(pair(i,j));        }    }    c.wait();    c.m_pixel.begin(); //in distCnC case, this gathers all items in one go    for (int i = 0; i < max_row; i++)     {        for (int j = 0; j < max_col; j++ )        {            c.m_pixel.get(pair(i,j), pixels[i*max_col + j]);        }    }        tbb::tick_count t1 = tbb::tick_count::now();    printf("Mandel %d %d %d in %g seconds/n", max_row, max_col, max_depth,           (t1-t0).seconds());        int check = 0;    for (int i = 0; i < max_row; i++)     {        for (int j = 0; j < max_col; j++ )        {            if (pixels[i*max_col + j ] == max_depth) check += (i*max_col +j );        }    }    printf("Mandel check %d /n", check);        if (verbose)    {        for (int i = 0; i < max_row; i++)         {            for (int j = 0; j < max_col; j++ )            {                if (pixels[i*max_col + j] == max_depth)                {                    std::cout << " ";                }                else                {                    std::cout << ".";                }            }            std::cout << std::endl;        }    }    return 0;}

void TurretShape::getCameraTransform(F32* pos,MatrixF* mat){   // Returns camera to world space transform   // Handles first person / third person camera position   if (isServerObject() && mShapeInstance)      mShapeInstance->animateNodeSubtrees(true);   if (*pos == 0) {      getRenderEyeTransform(mat);      return;   }   // Get the shape's camera parameters.   F32 min,max;   MatrixF rot;   Point3F offset;   getCameraParameters(&min,&max,&offset,&rot);   // Start with the current eye position   MatrixF eye;   getRenderEyeTransform(&eye);   // Build a transform that points along the eye axis   // but where the Z axis is always up.   {      MatrixF cam(1);      VectorF x,y,z(0,0,1);      eye.getColumn(1, &y);      mCross(y, z, &x);      x.normalize();      mCross(x, y, &z);      z.normalize();      cam.setColumn(0,x);      cam.setColumn(1,y);      cam.setColumn(2,z);      mat->mul(cam,rot);   }   // Camera is positioned straight back along the eye's -Y axis.   // A ray is cast to make sure the camera doesn't go through   // anything solid.   VectorF vp,vec;   vp.x = vp.z = 0;   vp.y = -(max - min) * *pos;   eye.mulV(vp,&vec);   // Use the camera node as the starting position if it exists.   Point3F osp,sp;   if (mDataBlock->cameraNode != -1)    {      mShapeInstance->mNodeTransforms[mDataBlock->cameraNode].getColumn(3,&osp);      getRenderTransform().mulP(osp,&sp);   }   else      eye.getColumn(3,&sp);   // Make sure we don't hit ourself...   disableCollision();   if (isMounted())      getObjectMount()->disableCollision();   // Cast the ray into the container database to see if we're going   // to hit anything.   RayInfo collision;   Point3F ep = sp + vec + offset;   if (mContainer->castRay(sp, ep,         ~(WaterObjectType | GameBaseObjectType | DefaultObjectType | sTriggerMask),         &collision) == true) {      // Shift the collision point back a little to try and      // avoid clipping against the front camera plane.      F32 t = collision.t - (-mDot(vec, collision.normal) / vec.len()) * 0.1;      if (t > 0.0f)         ep = sp + offset + (vec * t);      else         eye.getColumn(3,&ep);   }   mat->setColumn(3,ep);   // Re-enable our collision.   if (isMounted())      getObjectMount()->enableCollision();   enableCollision();   // Apply Camera FX.   mat->mul( gCamFXMgr.getTrans() );}

void SdMsgBaseCanvas::save(QTextStream & st, QString & warning) const{    nl_indent(st);    st << "to ";    dest->save(st, TRUE, warning);    nl_indent(st);#ifdef FORCE_INT_COORD    // note : << float bugged in Qt 3.3.3    st << "yz " << (int) y() << " " << (int) z();#else    QString sy, sz;    st << "yz " << sy.setNum(y()) << " " << sz.setNum(z());#endif    if (msg != 0) {        // not a lost, dest is duration        if (msg->deletedp()) {            warning += QString("<b>") + the_canvas()->browser_diagram()->full_name() +                       "</b> reference the class <b>" +                       ((SdDurationCanvas *) dest)->get_line()->get_obj()->get_class()->full_name() +                       "</b> deleted operation <b>" + msg->definition(TRUE, FALSE) + "</b><br>/n<br>/n";            if (warning[0] == '!') {                st << " msg ";                msg->save(st, TRUE, warning);            }            else {                st << " explicitmsg ";                save_string(msg->get_browser_node()->get_name(), st);            }        }        else {            st << " msg ";            msg->save(st, TRUE, warning);        }    }    else if (! explicit_msg.isEmpty()) {        st << " explicitmsg ";        save_string(explicit_msg, st);    }    else        st << " unspecifiedmsg";    if (stereotype != 0) {        nl_indent(st);        st << "stereotype ";        save_string(stereotype->get_name(), st);        save_xyz(st, stereotype, " xyz");    }    nl_indent(st);    st << "show_full_operations_definition " << stringify(show_full_oper)       << " drawing_language " << stringify(drawing_language)       << " show_context_mode " << stringify(show_context_mode);    if (label != 0) {        if (! args.isEmpty()) {            nl_indent(st);            st << "args ";            save_string(args, st);        }        nl_indent(st);        save_xy(st, label, "label_xy");    }}

int main(int argc, char *argv[]){  int ierr = 0, forierr = 0;  bool debug = false;#ifdef EPETRA_MPI  // Initialize MPI  MPI_Init(&argc,&argv);  int rank; // My process ID  MPI_Comm_rank(MPI_COMM_WORLD, &rank);  Epetra_MpiComm Comm( MPI_COMM_WORLD );#else  int rank = 0;  Epetra_SerialComm Comm;#endif  bool verbose = false;  // Check if we should print results to standard out  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;  int verbose_int = verbose ? 1 : 0;  Comm.Broadcast(&verbose_int, 1, 0);  verbose = verbose_int==1 ? true : false;  //  char tmp;  //  if (rank==0) cout << "Press any key to continue..."<< std::endl;  //  if (rank==0) cin >> tmp;  //  Comm.Barrier();  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting  int MyPID = Comm.MyPID();  int NumProc = Comm.NumProc();  if(verbose && MyPID==0)    cout << Epetra_Version() << std::endl << std::endl;  if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc		    << " is alive."<<endl;  bool verbose1 = verbose;  // Redefine verbose to only print on PE 0  if(verbose && rank!=0) 		verbose = false;  int NumMyEquations = 10000;  int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);  if(MyPID < 3)     NumMyEquations++;  // Construct a Map that puts approximately the same Number of equations on each processor  Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);    // Get update list and number of local equations from newly created Map  int* MyGlobalElements = new int[Map.NumMyElements()];  Map.MyGlobalElements(MyGlobalElements);  // Create an integer vector NumNz that is used to build the Petra Matrix.  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor  int* NumNz = new int[NumMyEquations];  // We are building a tridiagonal matrix where each row has (-1 2 -1)  // So we need 2 off-diagonal terms (except for the first and last equation)  for (int i = 0; i < NumMyEquations; i++)    if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))      NumNz[i] = 1;    else      NumNz[i] = 2;  // Create a Epetra_Matrix  Epetra_CrsMatrix A(Copy, Map, NumNz);  EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);  EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);    // Add  rows one-at-a-time  // Need some vectors to help  // Off diagonal Values will always be -1  double* Values = new double[2];  Values[0] = -1.0; 	Values[1] = -1.0;  int* Indices = new int[2];  double two = 2.0;  int NumEntries;  forierr = 0;  for (int i = 0; i < NumMyEquations; i++) {//.........这里部分代码省略.........

double leekesler::sdep(){    double tr = T/Tcr;    return tr*I() + J() - log(z());}

int main (int argc, char **argv){#if defined(_DIST_)    CnC::dist_cnc_init< HeatEquation_bl_context > AAA;#endif    //    CnC::debug::set_num_threads(1);    tbb::tick_count t0 = tbb::tick_count::now();    if( argc < 3 ) {        std::cerr << "expecting 2 arguments: <file> <blocksize>/n";        exit( 1 );    }    global_data d;    int& block_size = d.block_size;    int& Nx = d.Nx;    int& Nt = d.Nt;    int& N_bl = d.N_bl;    double& Xa = d.Xa;    double& Xb = d.Xb;    double& T0 = d.T0;    double& T1 = d.T1;    double& k = d.k;    double& hx = d.hx;    double& ht = d.ht;    block_size = atoi(argv[2]);    if(block_size < 1) {        std::cerr<<"Bad block size/n";        return 0;    }    {        std::ifstream from(argv[1]);        if( ! from ) {            std::cerr << "couldn't open " << argv[1] << std::endl;            exit( 2 );        }        from >> d.Xa >> d.Xb >> d.Nx >> d.T0 >> d.T1 >> d.Nt >> d.k;        from.close();    }    if( block_size > Nx+1 ) {        block_size = Nx+1;    } else {        Nx = ((Nx+1+block_size-1)/block_size)*block_size-1;    }    HeatEquation_bl_context c;    hx = ( Xb - Xa ) / Nx;    ht = ( T1 - T0 ) / Nt;    //    for ( int i = 0; i < Nx + 1; i++ ) { c.X.put(i,Xa+i*hx); }    //    for ( int i = 0; i < Nt + 1; i++ ) { c.T.put(i,T0+i*ht); }    N_bl = (Nx+1) / block_size;    c.gd.put(0,d);    for ( int i = 0; i < N_bl; i++ ) {        my_array<double> z(block_size);        for(int j = 0; j < block_size; j++ ){            double x = calc_x(block_size*i+j,Xa,Xb,Nx);            z.data[j] = Analitical( x, T0 );        }        c.H.put(Pair(0,i),z);    }    Pair p;    p.It = 1;    for ( int j = 0; j < N_bl; j++ ) {        p.Jx = j;        c.Tag.put( p );    }    // Wait for all steps to finish    c.wait();    tbb::tick_count t1 = tbb::tick_count::now();    std::cout<<"Time taken: "<< (t1-t0).seconds()<<" /n";    if (argc >= 4){        for (int i = 0; i <= Nt; i++)            {            for (int j = 0; j < N_bl; j++){                my_array<double> z;                c.H.get(Pair(i,j),z);                for (int k = j*block_size; k < j*block_size+block_size && k <= Nx; k++){                    printf("%.6lf ",double(z.data[k-j*block_size]));                }            }            puts("");        }    }    Eo(N_bl);    Eo(Nt);    return 0;}

int main(int argc, char *argv[]) {  Teuchos::GlobalMPISession mpiSession(&argc, &argv);  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.  int iprint     = argc - 1;  Teuchos::RCP<std::ostream> outStream;  Teuchos::oblackholestream bhs; // outputs nothing  if (iprint > 0)    outStream = Teuchos::rcp(&std::cout, false);  else    outStream = Teuchos::rcp(&bhs, false);  int errorFlag  = 0;  // *** Test body.  try {    std::string filename = "input.xml";    Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );    parlist->sublist("General").set("Inexact Hessian-Times-A-Vector",true);#if USE_HESSVEC    parlist->sublist("General").set("Inexact Hessian-Times-A-Vector",false);#endif    // Define Status Test    Teuchos::RCP<ROL::StatusTest<RealT> > status = Teuchos::rcp(new ROL::StatusTest<RealT>(*parlist));    *outStream << "/n/n" << ROL::ETestObjectivesToString(ROL::TESTOBJECTIVES_ROSENBROCK) << "/n/n";    // Initial Guess Vector     Teuchos::RCP<std::vector<RealT> > x0_rcp = Teuchos::rcp( new std::vector<RealT> );    ROL::StdVector<RealT> x0(x0_rcp);    // Exact Solution Vector    Teuchos::RCP<std::vector<RealT> > z_rcp = Teuchos::rcp( new std::vector<RealT> );    ROL::StdVector<RealT> z(z_rcp);    // Get Objective Function    Teuchos::RCP<ROL::Objective<RealT> > obj = Teuchos::null;    ROL::getTestObjectives<RealT>(obj,x0,z,ROL::TESTOBJECTIVES_ROSENBROCK);    // Get Dimension of Problem    int dim = Teuchos::rcp_const_cast<std::vector<RealT> >(                (Teuchos::dyn_cast<ROL::StdVector<RealT> >(x0)).getVector())->size();    parlist->sublist("General").sublist("Krylov").set("Iteration Limit", 2*dim);    // Iteration Vector    Teuchos::RCP<std::vector<RealT> > x_rcp = Teuchos::rcp( new std::vector<RealT> (dim, 0.0) );    ROL::StdVector<RealT> x(x_rcp);    x.set(x0);    // Error Vector    Teuchos::RCP<std::vector<RealT> > e_rcp = Teuchos::rcp( new std::vector<RealT> (dim, 0.0) );    ROL::StdVector<RealT> e(e_rcp);    e.zero();    for ( ROL::EDescent desc = ROL::DESCENT_STEEPEST; desc < ROL::DESCENT_LAST; desc++ ) {      parlist->sublist("Step").sublist("Line Search").sublist("Descent Method").set("Type", ROL::EDescentToString(desc));      *outStream << "/n/n" << ROL::EDescentToString(desc) << "/n/n";      for (ROL::ELineSearch ls = ROL::LINESEARCH_BACKTRACKING; ls < ROL::LINESEARCH_USERDEFINED; ls++) {        // Define Step        parlist->sublist("Step").sublist("Line Search").sublist("Line-Search Method").set("Type",ROL::ELineSearchToString(ls));        Teuchos::RCP<ROL::LineSearchStep<RealT> > step = Teuchos::rcp(new ROL::LineSearchStep<RealT>(*parlist));              // Define Algorithm        ROL::Algorithm<RealT> algo(step,status,false);        // Run Algorithm        x.set(x0);        algo.run(x, *obj, true, *outStream);        // Compute Error        e.set(x);        e.axpy(-1.0,z);        *outStream << "/nNorm of Error: " << e.norm() << "/n";        //errorFlag += (int)(e.norm() < std::sqrt(ROL::ROL_EPSILON));       }          }  }  catch (std::logic_error err) {    *outStream << err.what() << "/n";    errorFlag = -1000;  }; // end try  if (errorFlag != 0)    std::cout << "End Result: TEST FAILED/n";  else    std::cout << "End Result: TEST PASSED/n";  return 0;}

// ------------------------------------------------------------// SFCPartitioner implementationvoid SFCPartitioner::_do_partition (MeshBase& mesh,				    const unsigned int n){  libmesh_assert_greater (n, 0);  // Check for an easy return  if (n == 1)    {      this->single_partition (mesh);      return;    }// What to do if the sfcurves library IS NOT present#ifndef LIBMESH_HAVE_SFCURVES  libmesh_here();  libMesh::err << "ERROR: The library has been built without"    << std::endl	        << "Space Filling Curve support.  Using a linear" << std::endl	        << "partitioner instead!" << std::endl;  LinearPartitioner lp;  lp.partition (mesh, n);// What to do if the sfcurves library IS present#else  START_LOG("sfc_partition()", "SFCPartitioner");  const unsigned int n_active_elem = mesh.n_active_elem();  const unsigned int n_elem        = mesh.n_elem();  // the forward_map maps the active element id  // into a contiguous block of indices  std::vector<unsigned int>    forward_map (n_elem, libMesh::invalid_uint);  // the reverse_map maps the contiguous ids back  // to active elements  std::vector<Elem*> reverse_map (n_active_elem, NULL);  int size = static_cast<int>(n_active_elem);  std::vector<double> x      (size);  std::vector<double> y      (size);  std::vector<double> z      (size);  std::vector<int>    table  (size);  // We need to map the active element ids into a  // contiguous range.  {//     active_elem_iterator       elem_it (mesh.elements_begin());//     const active_elem_iterator elem_end(mesh.elements_end());    MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();    const MeshBase::element_iterator elem_end = mesh.active_elements_end();    unsigned int el_num = 0;    for (; elem_it != elem_end; ++elem_it)      {	libmesh_assert_less ((*elem_it)->id(), forward_map.size());	libmesh_assert_less (el_num, reverse_map.size());	forward_map[(*elem_it)->id()] = el_num;	reverse_map[el_num]           = *elem_it;	el_num++;      }    libmesh_assert_equal_to (el_num, n_active_elem);  }  // Get the centroid for each active element  {//     const_active_elem_iterator       elem_it (mesh.const_elements_begin());//     const const_active_elem_iterator elem_end(mesh.const_elements_end());    MeshBase::element_iterator       elem_it  = mesh.active_elements_begin();    const MeshBase::element_iterator elem_end = mesh.active_elements_end();    for (; elem_it != elem_end; ++elem_it)      {	const Elem* elem = *elem_it;	libmesh_assert_less (elem->id(), forward_map.size());	const Point p = elem->centroid();	x[forward_map[elem->id()]] = p(0);	y[forward_map[elem->id()]] = p(1);	z[forward_map[elem->id()]] = p(2);      }  }  // build the space-filling curve  if (_sfc_type == "Hilbert")    Sfc::hilbert (&x[0], &y[0], &z[0], &size, &table[0]);//.........这里部分代码省略.........

int main(int argc, char** argv){  char port[] = "27015";  char ip[] = "10.0.0.67";  int numbytes;  char buf[MAXDATASIZE];  // Get input parameters  for(int a=0; a < argc; a++)    {      if( strcmp( argv[a], "--ip" ) == 0 )	{	  strcpy(ip, argv[a+1]);	}      if( strcmp( argv[a], "--port" ) == 0 )	{	  strcpy(port, argv[a+1]);	}    }  TCP myClient = TCP(port, ip);  ros::init(argc, argv, "joints");  ros::NodeHandle n;  std::string receivedString;  std::string attribute;  std::stringstream ss;  std::istringstream iss;  tf::TransformBroadcaster br;  tf::Transform transform;  ros::Rate loop_rate(30);  if( myClient.Connect() == 0 )    {      // Receive data      numbytes = recv(myClient.s,buf,MAXDATASIZE-1,0);      while (numbytes != 0 && ros::ok())	{	  numbytes = recv(myClient.s,buf,MAXDATASIZE-1,0);	  buf[numbytes]='/0';	  receivedString.assign(buf);	  // debug	  // std::cout << "Ahoy! Received " << numbytes << " bytes. " << receivedString << std::endl;	  std::istringstream iss(receivedString);	  std::vector<double> x (25, 0.0);	  std::vector<double> y (25, 0.0);	  std::vector<double> z (25, 0.0);    std::vector<std::vector<double> > X (6, x);    std::vector<std::vector<double> > Y (6, y);    std::vector<std::vector<double> > Z (6, z);    std::vector<double> qx (25, 0.0);	  std::vector<double> qy (25, 0.0);	  std::vector<double> qz (25, 0.0);	  std::vector<double> qw (25, 1.0);    std::vector<std::vector<double> > QX (6, qx);    std::vector<std::vector<double> > QY (6, qy);    std::vector<std::vector<double> > QZ (6, qz);    std::vector<std::vector<double> > QW (6, qw);	  int j;	  int count = 6;	  int id;	  std::vector<int> foundIdx;    std::string frame;    std::string joint;      /* code */	  while(!iss.eof())	    {      // joint = "joint_";      iss >> attribute;      frame = "person_";	      if( strcmp(attribute.c_str(), "count:") == 0 )		{		  iss >> count;		  // ROS_INFO("count %d",count);		}	      else if( strcmp(attribute.c_str(), "id:") == 0 )		{		  iss >> id;//.........这里部分代码省略.........

int main(int argc, char *argv[]){#ifdef HAVE_MPI  MPI_Init(&argc, &argv);  Epetra_MpiComm Comm(MPI_COMM_WORLD);#else  Epetra_SerialComm Comm;#endif    bool verbose = (Comm.MyPID() == 0);  // set global dimension to 5, could be any number  int NumGlobalElements = 5;  // create a map  Epetra_Map Map(NumGlobalElements,0,Comm);  // local number of rows  int NumMyElements = Map.NumMyElements();  // get update list  int * MyGlobalElements = Map.MyGlobalElements( );  // ============= CONSTRUCTION OF THE MATRIX ===========================  // Create a Epetra_Matrix  Epetra_CrsMatrix A(Copy,Map,3);  // Add  rows one-at-a-time  double *Values = new double[2];  Values[0] = -1.0; Values[1] = -1.0;  int *Indices = new int[2];  double two = 2.0;  int NumEntries;  for( int i=0 ; i<NumMyElements; ++i ) {    if (MyGlobalElements[i]==0) {	Indices[0] = 1;	NumEntries = 1;    } else if (MyGlobalElements[i] == NumGlobalElements-1) {      Indices[0] = NumGlobalElements-2;      NumEntries = 1;    } else {      Indices[0] = MyGlobalElements[i]-1;      Indices[1] = MyGlobalElements[i]+1;      NumEntries = 2;    }    A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);    // Put in the diagonal entry    A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i);  }    // Finish up  A.FillComplete();  // ================ CONSTRUCTION OF VECTORS =======================    // build up two distributed vectors q and z, and compute  // q = A * z  Epetra_Vector q(A.RowMap());  Epetra_Vector z(A.RowMap());  // Fill z with 1's  z.PutScalar( 1.0 );  // ================ USE OF TIME AND FLOPS =========================    Epetra_Flops counter;  A.SetFlopCounter(counter);  Epetra_Time timer(Comm);  A.Multiply(false, z, q); // Compute q = A*z  double elapsed_time = timer.ElapsedTime();  double total_flops =counter.Flops();  if (verbose)    cout << "Total ops: " << total_flops << endl;  double MFLOPs = total_flops/elapsed_time/1000000.0;  if (verbose)    cout << "Total MFLOPs  for mat-vec = " << MFLOPs << endl<< endl;  double dotProduct;  z.SetFlopCounter(counter);  timer.ResetStartTime();  z.Dot(q, &dotProduct);  total_flops =counter.Flops();  if (verbose)    cout << "Total ops: " << total_flops << endl;  elapsed_time = timer.ElapsedTime();  if (elapsed_time != 0.0)    MFLOPs = (total_flops / elapsed_time) / 1000000.0;  else    MFLOPs = 0;  if (verbose)  {    cout << "Total MFLOPs for vec-vec = " << MFLOPs << endl<< endl;    cout << "q dot z = " << dotProduct << endl;  }//.........这里部分代码省略.........

unsigned long CSysSolve::CG_LinSolver(const CSysVector & b, CSysVector & x, CMatrixVectorProduct & mat_vec,                                           CPreconditioner & precond, su2double tol, unsigned long m, bool monitoring) {	int rank = 0;#ifdef HAVE_MPI	MPI_Comm_rank(MPI_COMM_WORLD, &rank);#endif    /*--- Check the subspace size ---*/  if (m < 1) {    if (rank == MASTER_NODE) cerr << "CSysSolve::ConjugateGradient: illegal value for subspace size, m = " << m << endl;#ifndef HAVE_MPI    exit(EXIT_FAILURE);#else    MPI_Barrier(MPI_COMM_WORLD);    MPI_Abort(MPI_COMM_WORLD,1);    MPI_Finalize();#endif  }    CSysVector r(b);  CSysVector A_p(b);    /*--- Calculate the initial residual, compute norm, and check if system is already solved ---*/  mat_vec(x, A_p);    r -= A_p; // recall, r holds b initially  su2double norm_r = r.norm();  su2double norm0 = b.norm();  if ( (norm_r < tol*norm0) || (norm_r < eps) ) {    if (rank == MASTER_NODE) cout << "CSysSolve::ConjugateGradient(): system solved by initial guess." << endl;    return 0;  }    su2double alpha, beta, r_dot_z;  CSysVector z(r);  precond(r, z);  CSysVector p(z);    /*--- Set the norm to the initial initial residual value ---*/  norm0 = norm_r;    /*--- Output header information including initial residual ---*/  int i = 0;  if ((monitoring) && (rank == MASTER_NODE)) {    WriteHeader("CG", tol, norm_r);    WriteHistory(i, norm_r, norm0);  }    /*---  Loop over all search directions ---*/  for (i = 0; i < (int)m; i++) {        /*--- Apply matrix to p to build Krylov subspace ---*/    mat_vec(p, A_p);        /*--- Calculate step-length alpha ---*/    r_dot_z = dotProd(r, z);    alpha = dotProd(A_p, p);    alpha = r_dot_z / alpha;        /*--- Update solution and residual: ---*/    x.Plus_AX(alpha, p);    r.Plus_AX(-alpha, A_p);        /*--- Check if solution has converged, else output the relative residual if necessary ---*/    norm_r = r.norm();    if (norm_r < tol*norm0) break;    if (((monitoring) && (rank == MASTER_NODE)) && ((i+1) % 5 == 0)) WriteHistory(i+1, norm_r, norm0);        precond(r, z);        /*--- Calculate Gram-Schmidt coefficient beta,		 beta = dotProd(r_{i+1}, z_{i+1}) / dotProd(r_{i}, z_{i}) ---*/    beta = 1.0 / r_dot_z;    r_dot_z = dotProd(r, z);    beta *= r_dot_z;        /*--- Gram-Schmidt orthogonalization; p = beta *p + z ---*/    p.Equals_AX_Plus_BY(beta, p, 1.0, z);  }      if ((monitoring) && (rank == MASTER_NODE)) {    cout << "# Conjugate Gradient final (true) residual:" << endl;    cout << "# Iteration = " << i << ": |res|/|res0| = "  << norm_r/norm0 << "./n" << endl;  }  //  /*--- Recalculate final residual (this should be optional) ---*///  mat_vec(x, A_p);//  r = b;//  r -= A_p;//  su2double true_res = r.norm();//  //  if (fabs(true_res - norm_r) > tol*10.0) {//    if (rank == MASTER_NODE) {//      cout << "# WARNING in CSysSolve::ConjugateGradient(): " << endl;//      cout << "# true residual norm and calculated residual norm do not agree." << endl;//      cout << "# true_res - calc_res = " << true_res - norm_r << endl;//.........这里部分代码省略.........

GLdouble GLVector::length()const{    return sqrt(x()*x() + y()*y() + z()*z());}

double leekesler::Pp(){    return 8314.3*z()*Rho*T/Mw;}